Ultrafast formation dynamics of D3+ from the light-driven bimolecular reaction of the D2–D2 dimer

The light-driven formation of trihydrogen cation has been attracting considerable attention because of its important role as an initiator of chemical reactions in interstellar clouds. To understand the formation dynamics, most previous studies focused on creating H3+ or D3+ from unimolecular reactions of various organic molecules. Here we observe and characterize the ultrafast formation dynamics of D3+ from a bimolecular reaction, using pump–probe experiments that employ ultrashort laser pulses to probe its formation from a D2–D2 dimer. Our molecular dynamics simulations provide an intuitive representation of the reaction dynamics, which agree well with the experimental observation. We also show that the emission direction of D3+ can be controlled using a tailored two-colour femtosecond laser field. The underlying control mechanism is in line with what is known from the light control of electron localization in the bond breaking of single molecules. H3+ and D3+ serve as initiators of many chemical reactions in interstellar clouds. Now the ultrafast formation dynamics of D3+ from a light-driven bimolecular reaction starting from D2–D2 dimers have been measured. It has also been shown that the emission direction of D3+ can be controlled by driving the reaction with a more complex two-colour laser pulse.

The light-driven formation of trihydrogen cation has been attracting considerable attention because of its important role as an initiator of chemical reactions in interstellar clouds. To understand the formation dynamics, most previous studies focused on creating H 3 + or D 3 + from unimolecular reactions of various organic molecules. Here we observe and characterize the ultrafast formation dynamics of D 3 + from a bimolecular reaction, using pump-probe experiments that employ ultrashort laser pulses to probe its formation from a D 2 -D 2 dimer. Our molecular dynamics simulations provide an intuitive representation of the reaction dynamics, which agree well with the experimental observation. We also show that the emission direction of D 3 + can be controlled using a tailored two-colour femtosecond laser field. The underlying control mechanism is in line with what is known from the light control of electron localization in the bond breaking of single molecules.
When exposed to a strong laser field, molecular potential surfaces suffer serious deformation, leading to chemical bond cleavage and formation. This gives rise to new reaction products, even in unimolecular reactions. The light-driven formation dynamics of the trihydrogen cation H 3 + (the simplest but most important triatomic ion) and its isotope D 3 + have attracted much attention in the past two decades 1-14 due to the important role of H 3 + as the initiator for most chemical reactions in interstellar clouds 15-18 , including those leading to water and organic molecules. For almost one century, it has generally been believed that the formation of this trihydrogen cation in interstellar molecular clouds is triggered by cosmic ray ionization of hydrogen molecules 19 (that is, H 2 + + H 2 → H 3 + + H). This reaction is regarded as the start of ion-neutral reaction chains in interstellar clouds. Previous experimental studies of this bimolecular reaction were mostly aimed at obtaining the kinetic energy and angular distribution of the products, the reaction rate coefficients and the state-selective cross-sections [20][21][22] . However, the time information of trihydrogen cation formation dynamics was obscured in the beam and spectroscopy measurements [23][24][25][26][27][28] .
Recently, the ultrafast formation dynamics of the trihydrogen cation from a laser-induced unimolecular reaction in organic molecules have attracted much attention 7-12 since an unambiguous time trigger of the reaction can be identified. Rich formation dynamics were revealed, including H 2 molecule roaming 12,13 and hydrogen atom migration 28,29 . Although these works provide abundant information to aid our understanding of the formation dynamics of trihydrogen cations from the fragmentation of large organic molecules, a more concerning reaction in the interstellar clouds would be to produce H 3 + or D 3 + from simpler elements or molecules rather than the opposite.
Article https://doi.org/10.1038/s41557-023-01230-0 Interestingly, branch B showed continuously decreasing energy with an increase in the delay, indicating a joint contribution from the pump and probe pulses. Formation of the (D 3 + , D) pair was first triggered by the pump pulse and started to dissociate. The probe pulse then ionized the neutral D atom while it was separating from D 3 + . The kinetic energy release accumulated after the probe pulse when the charged D + and D 3 + ions repelled each other. The kinetic energy released decreased with the delay as the distance between D and D 3 + increased. For the time-independent branch A, double ionization of the dimer occurred within the single-pulse duration of 35 fs, resulting in the formation of (D 3 + , D + ) fragments. In contrast, the low-energy branch B showed a time-dependent yield that increased with delay. This meant that the probe pulse not only ionized the neutral D atom but also interfered with the formation of D 3 + , which may have led to three-or complete four-body fragmentation. To extract the formation time of D 3 + from the single-ionization pathway, we plotted the D 3 + yield as a function of the time delay in Fig. 2c. The increase in yield clearly indicates that more D 3 + was formed at larger delays, when the D 3 + was stably formed to the extent that the probe pulse could not interrupt the formation process any more. By fitting with an exponential function of y = A 0 exp(−t/τ) + y 0 , where A 0 is the amplitude, τ is the time constant and y 0 is the offset caused by the high-energy pathway, we obtained an effective formation time of τ = 139 ± 7 fs, as shown in Fig. 2c.
Very recently, trihydrogen cation formation from water molecules attached to nanoparticle surfaces was observed 14 , which serves as a catalytic environment to boost trihydrogen cation formation. However, there are still many relevant questions to be answered. Can this reaction happen without the catalysis of nanoparticles? How fast is the bimolecular reaction to produce trihydrogen cations compared with the unimolecular reaction in organic molecules? And most importantly, can we manipulate the trihydrogen cation formation dynamics in the bimolecular reaction?
In this Article, we investigate the ultrafast formation dynamics of D 3 + from a bimolecular reaction of a gas phase D 2 -D 2 dimer driven by ultrashort laser pulses. We find that D 3 + can be created from a single D 2 -D 2 dimer without any interaction with the surrounding environment (for example, nanoparticles 14 or analagous nanoscale dust in the interstellar clouds). By scanning the time delay between the pump and probe pulses, we identify two formation pathways corresponding to the double-and single-ionization processes, respectively. The fast double-ionization pathway, denoted as D 2 -D 2 + nħω → D 2 + + D + D + → D 3 + + D + , is induced by a single pulse with a duration of ~35 fs, where n is the number of absorbed photons of angular frequency ω. In contrast, the slow single-ionization pathway, denoted as D 2 -D 2 + nħω → D 2 + + D 2 → D 3 + + D, takes ~139 fs. By performing a molecular dynamics simulation, we track the reaction dynamics in time. Furthermore, by constructing a tailored two-colour femtosecond laser field, we realize control of the emission direction of the outgoing D 3 + from the bimolecular reaction by finely tuning the relative phase between the two colours 30-33 .

Results and discussion
As illustrated in Fig. 1a, two sets of pump-probe experiments were performed to investigate the formation dynamics of D 3 + . We only give a brief description here (see Methods for details). For the sake of simplicity, we started with a single-colour (790 nm) pump-probe experiment to clock the formation dynamics. The pump and probe pulses were polarized along the z and y axes respectively, such that the products from the pump and probe could be distinguished. The pump pulse initiated the fast and slow D 3 + formation pathways, while the probe pulse disrupted the formation process 11,12 and ionized the neutral D atom in the slow pathway to make it detectable. For the two-colour pump-probe experiment (pump: 790 and 395 nm; probe: 790 nm), the probe pulse was fixed at a time delay of ~1.2 ps with respect to the phase-controlled two-colour pump pulse. By scanning the relative phase between the two-colour pump pulses, we controlled the emission direction of the formed D 3 + . Figure 2a shows the measured coincidence map of the two-body fragmentation channels, where the horizontal and vertical axes are the time of flight of the first and second ions, respectively. Different two-body fragmentation channels induced by Coulomb explosion could clearly be identified along the off-diagonal axis. The overwhelming signal was from the double ionization of single D 2 molecules (that is, D 2 + nħω → D + + D + ), termed as (D + , D + ) hereafter. Besides the prominent (D + , D + ) channel, two fragmentation channels of (D 2 + , D 2 + ) and (D 3 + , D + ) appeared only when the molecular jet was pre-cooled to produce D 2 -D 2 dimers.

Timing the formation of D 3 +
To understand the formation dynamics of the (D 3 + , D + ) channel, we first performed a single-colour pump-probe experiment. Figure 2b displays the kinetic energy spectrum of the (D 3 + , D + ) channel as a function of the time delay between the pump and probe pulses. Two branches, denoted as A and B hereafter, could clearly be distinguished. Branch A, with a kinetic energy peaking at ~4.0 eV, showed a flat distribution that was delay independent. This high-energy channel was created from direct double ionization by the pump pulse. It was also present when only a single pulse was applied (blue curve in the left panel of Fig. 2b). Potential energy (eV)

Fig. 1 | Schematic of the experiments. a,
A supersonic gas jet of D 2 -D 2 molecular dimers was exposed to pump-probe pulses in an ultrahigh vacuum chamber of a cold target recoil ion momentum spectrometer. Two sets of pump-probe schemes were implemented. Scheme (1) was a single-colour pump-probe measurement with a continuously scanned time delay t, whereas scheme (2) was a phase-scanning two-colour experiment with a fixed pump-probe time delay T. Inset, schematic of the two-colour field for the relative phase ϕ L = 0 or π and the corresponding phase-dependent directional dissociative ionization of one D 2 molecule in the D 2 -D 2 dimer. TOF, time of flight. b, Potential energy curves of 1sσ g and 2pσ u of D 2 + as a function of internuclear distance R. The dominating net-2ω FW and 1ω SH dissociation pathways are illustrated with vertical red and blue arrows, representing 790 and 395 nm photons, respectively.
Article https://doi.org/10.1038/s41557-023-01230-0 In the following, we analyse the reaction geometry of branches A and B. Figures 2d and 2e display the two-dimensional momentum distribution of the (D 3 + , D + ) channel in the y-z plane at different time delays, where the outer and inner rings correspond to branches A and B, respectively. The outer ring showed a pronounced anisotropy along the polarization of the pump pulse at large delays, while the inner ring was almost isotropic, as shown in Fig. 2f. The difference in angular distribution indicates that the two pathways undergo different reaction dynamics. For the fast pathway, the two electrons were simultaneously removed from two D 2 molecules constitute the dimer. If the two D 2 + molecular ions are at the electronic ground state, this will lead to the two-body Coulomb explosion of (D 2 + , D 2 + ). Only if one of the two D 2 + molecular ions is excited to the dissociative state can D 3 + be formed. This can happen when the neutral D atom from the dissociation of the excited D 2 + is flying towards the other D 2 + (that is, D 2 -D 2 + nħω → D 2 + + D + D + → D 3 + + D + ). This scenario is supported by the strong polarization dependency of the outer ring, since dissociative ionization is strongly favoured when the molecular axis of the dissociating D 2 + is parallel to the laser polarization 34 .
For the slow pathway, only one electron is removed from one of the D 2 molecules within the dimer by the pump pulse. We observed that the angular distribution of branch B was almost isotropic. These two pieces of evidence indicate that the ionized D 2 + is in the ground state and launches the effective ion-neutral exothermic reaction of D 2 -D 2 + nħω → D 2 + + D 2 → D 3 + + D. Ionization to the non-dissociative ground state is slightly enhanced when the molecular axis is parallel to the laser polarization 35 . This polarization dependence is finally  https://doi.org/10.1038/s41557-023-01230-0 averaged out due to the fast rotation of both the D 2 + molecular ion and the D 2 molecule, of which the rotational revival is ~550 fs [36][37][38] . To this end, we have revealed that the two formation pathways of D 3 + have different formation times and reaction dynamics.

Molecular dynamics simulation
To obtain a detailed picture of the reaction dynamics of D 3 + formation, we performed a molecular dynamics simulation. Previous studies have shown that there are several structural configurations of the neutral D 2 -D 2 dimer 39 . For the sake of simplicity, we first investigated reactions starting in the T-shape configuration, where two molecular axes are perpendicular to each other with the lowest energy at equilibrium. Figure 3 shows a typical reaction trajectory starting close to the equilibrium with zero initial velocity and an artificially frozen T-shape geometry. By defining two reaction coordinates of R 1 and R 2 , the potential energy surface of the T-type cationic (D 2 -D 2 ) + dimer is shown in Fig. 3 (see Methods for computational details). The initial geometry of the neutral dimer is indicated by step 1, with R 1 = 7.635 atomic units (a.u.) and R 2 = 1.390 a.u., respectively. Upon Frank-Condon excitation, the cationic state initially around the equilibrium geometry of the neutral state is populated. Under the force field of the cationic state, the nuclear dynamics evolve along the trajectory, as shown by the orange curve. Clearly, the D 2 + ion (one of the D 2 molecules whose electron is removed) starts to vibrate around the new equilibrium internuclear distance of R 2 = 2.0 a.u. At the same time, the D 2 + ion moves to approach the neighbouring D 2 molecule. The D 3 + ion is formed when the neutral D atom departs away at a longer time delay, denoted as step 3.
Next, we performed a more realistic on-the-fly nonadiabatic molecular dynamics simulation with a swarm of 300 nuclear trajectories (see Methods). During the formation of D 3 + , the D + ion approaches the neighbouring D 2 molecule. The other D atom departs simultaneously. Figure 4a-e shows the approaching and departing distances as a function of time for five structural configurations of T, L, H, X and Z shapes, respectively. We judged that D 3 + was formed when the approaching distance was smaller than 4 a.u. and the departing distance was larger than 5 a.u. Statistically, 64, 79, 21, 60 and 0.7% trajectories converged to form D 3 + in T, L, X, Z and H shapes, respectively. Obviously, the H shape could hardly form D 3 + , which mostly dissociated into three-body fragments of (D 2 , D, D + ). By collecting the trajectories that led to D 3 + for the T, L, X and Z shapes, we obtained formation times of 125, 245, 200 and 100 fs, respectively. Our simulation results clearly show that the formation time of D 3 + strongly relies on its initial structural configuration. The average formation time after considering the contribution from all of the above shapes was 150 fs, as shown in Fig. 4f. It is noticeable that, from the molecular dynamics simulation, the structure of (D 2 -D 2 ) + rearranges quickly due to rotation during the evolution to form (D 3 + , D) (see Supplementary Videos 1-3); thus, the probe process will not show strong polarization dependence. As a result, the experimentally measured formation time of 139 fs results from the average contribution from all possible configurations, well in agreement with the time scale predicted by our simulation.

Controlling the formation of D 3 +
After determining the intermediate products of the two pathways, we proceeded to control the bimolecular reaction by manipulating the emission direction of D 3 + using a two-colour field as the pump pulse, as illustrated in Fig. 1a. Here the probe pulse was fixed at a time delay of 1.2 ps and only served to ionize the neutral D atom in pathway B,   since the formation of D 3 + was completed and would not be disrupted at such a long delay. Figure 5a shows the momentum p z of D 3 + from the (D 3 + , D + ) channel as a function of the relative phase between the two colours. It is very clear that for high-energy branch A, with |p z | > 20 a.u., the D 3 + emitting upwards (p z > 0) or downwards (p z < 0) showed strong phase-dependent modulation. At phase ϕ L = 0 more D 3 + ions were emitted downwards, whereas at phase ϕ L = π the D 3 + ions were preferentially emitted upwards. In contrast, low-energy branch B showed a phase-independent distribution. This is easy to understand since the free rotation of both D 2 + molecular ions and D 2 molecules will smear out the modulation even if the phase information is tagged upon ionization. To quantify the directional emission of D 3 + , we define the asymmetry parameter A(|p z |, ϕ L ) = [N(p z+ , ϕ L ) − N(p z− , ϕ L )]/[N(p z+ , ϕ L ) + N(p z− , ϕ L )], where N(p z+ , ϕ L ) (or N(p z− , ϕ L )) is the measured D 3 + yield for p z > 0 (or p z < 0) at phase ϕ L . Figure 5b shows the momentum phase-dependent asymmetry spectrum of D 3 + calculated from Fig. 5a, where the phase modulation is more visible. The integrated phase-dependent asymmetry parameters for the high-energy part (23 a.u. < |p z | < 27 a.u.) are displayed in Fig. 5c, where the solid curve is the fitting sine function.
To this end, we have shown that the bimolecular reaction to produce D 3 + from a D 2 -D 2 dimer in the gas phase can be well controlled.
To reveal the underlying physics of the phase-dependent asymmetric emission of D 3 + in branch A, we compared the same result for D + from dissociative single ionization of the D 2 monomer with that from the (D 3 + , D + ) channel, as shown in Fig. 5d-i. As discussed above, the high-energy pathway to create D 3 + proceeds via dissociative single ionization of one of the D 2 molecules accompanied by single ionization of the other one; thus, the phase dependence of D + (or D 3 + ) from the (D 3 + , D + ) channel should be in (or out of) phase with that for single D 2 molecules. Figure 5d-f shows the phase dependence of the emission direction of D + from the (D 3 + , D + ) channel, which is opposite to that of D 3 + due to the momentum conservation. Figure 5g shows the measured phase-dependent momentum p z of D + from dissociative single ionization of the D 2 monomer (that is, D 2 + nħω → D + + D), denoted as (D + , D). Generally, several dissociation pathways can be identified, as illustrated in Fig. 1b. As shown by the one-dimensional momentum spectrum integrated over the laser phase in Fig. 5g, the here-observed (D + , D) is dominated by the net-2ω FW and 1ω SH pathways ending with protons of momentum |p z | > 10 a.u., where ω FW and ω SH represent the photon energies of the fundamental wave of 790 nm and its second harmonic of 395 nm, respectively. The asymmetric emission of D + requires superposition of nuclear wave packets with the same kinetic energies but opposite parities. These wave packets are created by dissociation along two potential surfaces of 1sσ g (net-2ω FW pathway) and 2pσ u (1ω SH pathway). Figure 5h shows the asymmetry parameter A(|p z |, ϕ L ) for the emitted protons from the interference between the net-2ω FW and 1ω SH pathways. The interferences of theses pathways were comprehensively discussed in previous studies [30][31][32][33] . Figure 5i shows the one-dimensional phase-dependent modulation of the asymmetry of the proton emission with |p z | > 10 a.u., which agrees well with that of D + from the (D 3 + , D + ) channel. This confirms that the formation of (D 3 + , D + ) in branch A (Fig. 2b) results from the dissociative single ionization of one of the D 2 molecules in the dimer. Thus, we can control the formation dynamics of D 3 + by controlling the directional bond breaking of the dissociation process.
Control over D 3 + formation can be well reproduced using the molecular dynamics simulations. This can be achieved by populating onto different cationic states of (D 2 -D 2 ) 2+ . We employed Mulliken charge analysis to determine the responsible dicationic state of (D 2 -D 2 ) 2+ that leads to the desired product. It was found that populating onto the ground state of (D 2 -D 2 ) 2+ leads to the formation of D 3 + , while populating onto the excited state leads to fragmentation of the dimer. Control over the D 3 + formation is thus achieved by transient charge localization facilitated by populating the two states 40 . When the laser field is in phase, it leads to the formation of D 3 + with efficiencies of 66 and 69% for the T and L shapes, respectively. When the laser field is out of phase, the formation of D 3 + drastically decreases to minimal levels of 0.3 and 0.6% for the T and L shapes, respectively. We note that control for the H and X types cannot be achieved since no matter which direction a D 2 + cation dissociates in, the other D 2 + cation is at no parity, leading to identical yields. Intuitively, D 3 + can only be produced when the dissociating D atom is emitting towards the neighbouring D 2 + , which is preferable for the T and L shapes. Otherwise, it will lead to the three-body breakup into (D 2 + , D + , D).

Conclusion
In conclusion, we show that the D 3 + ion can be produced from a bimolecular reaction of a D 2 -D 2 dimer in the gas phase. The formation time varies for different reaction pathways. Compared with the fast double-ionization pathway occurring within 35fs, the slow single-ionization pathway takes ~139 fs. Furthermore, we achieved control over the formation dynamics of the D 3 + ion by manipulating its emission direction using a tailored two-colour laser field. Our results are interesting not only for the field of photochemistry, but also astronomy, to understand the most fundamental light-induced bimolecular process in interstellar clouds, which was crucial for the start of the Universe.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41557-023-01230-0.

Experimental design
The experiments were performed with ultrashort femtosecond laser pulses (35 fs; 790 nm; 10 kHz) produced from a multipass amplifier Ti:sapphire laser system. The pulses were separated into pump and probe arms via a beam splitter with an intensity ratio of 7:3 and later recombined in a Mach-Zehnder-type interferometer configuration. As depicted in Fig. 1a, in both sets of pump-probe experiments, the pump and probe pulses were linearly polarized along the z and y axes, respectively. In the single-colour pump-probe experiment, the time delay between pump and probe pulse was continuously controlled by a motorized delay stage, which scanned from 60-2,060 fs with a step of 20 fs. The peak intensities of the pump and probe pulses were estimated to be ~2.0 × 10 14 W cm −2 . In the two-colour pump-probe experiment, the time delay between the pump and probe pulses was fixed at ~1.2 ps. The pump pulse was a phase-controlled linearly polarized two-colour laser pulse generated in a collinear scheme by feeding a fundamental wave pulse (790 nm; y polarized) into a 150-μm-thick β-barium borate crystal to produce the second harmonic pulse (395 nm; z polarized). The group delay between the fundamental wave and second harmonic pulse was compensated for by inserting birefringent α-barium borate crystals. The polarization of the fundamental wave pulse was rotated using a dual-wavelength wave plate to be parallel to that of the second harmonic pulse. The electric field of the linearly polarized two-colour field is expressed by E(t) = E ω (t)cos(ωt) + E 2ω (t)cos(2ωt + ϕ L ), where E ω (t) and E 2ω (t) are the envelope of the fundamental wave and second harmonic field and ϕ L is the relative phase between the fundamental wave and second harmonic components. The relative phase ϕ L was finely tuned by controlling the inserted thickness of a pair of fused silica wedges via a delay stage. The absolute value of ϕ L was calibrated by analysing the phase-dependent emission direction of C + from CO molecules, which is more likely ionized when the laser field points from C to O 41 . The peak intensities of the fundamental wave, second harmonic and probe pulse in the interaction region were estimated to be ~8.5 × 10 13 , 3.9 × 10 13 and 1.5 × 10 14 W cm −2 , respectively.
The pump and probe laser pulses were focused onto a supersonic molecular beam using a concave mirror with a focal length of 75 mm inside the ultrahigh vacuum chamber of a cold target recoil ion momentum spectrometer 42,43 . The gas of D 2 molecules was expanded through a 5-μm-diameter nozzle pre-cooled to 85 K with a driving pressure of 2.5 bar, generating a molecular beam composed of D 2 monomers and dimers. The nozzle was cooled using a closed-cycle cryostat and the temperature of the nozzle was controlled using a heating resistor and a temperature controller, giving rise to a stable temperature with a precision of 0.1 K. The photoionization-created ionic fragments were guided by a weak homogeneous electric field (~6 V cm −1 ) to a time-and position-sensitive microchannel plate detector at the end of the spectrometer. We constructed the three-dimensional momenta and kinetic energies of the ions from their times of flight and positions of impact. The Coulomb explosion channels of the multiply ionized dimer were unambiguously identified based on the momentum conservation of the measured ionic fragments.

Molecular dynamics simulation
Computational simulations were carried out to investigate the dynamics during the chemical reactions, which provided the time-resolved information of the bond formation and bond cleavage. In agreement with the literature, we found a few stable configurations of the neutral D 2 dimer. To obtain a list of these stable configurations, we randomly placed the four D atoms in space and relaxed the geometry using the OpenMolcas package 44 . We found stable configurations including the T, L, H, Z and X types, in line with recent experiments 39 . The potential energy surfaces were calculated via the complete active space self-consistent field (CASSCF) method at the CAS(4,4), CAS (3,4) and CAS(2,4) levels with the aug-cc-pVTZ basis for the neutral, cationic and dicationic dimer states, respectively. We used the equilibrium geometry of the neutral dimer as the starting configuration. Upon Frank-Condon excitation, the nuclear trajectory would evolve along the potential energy surface of the respective cationic or dicationic state, which was initially at the equilibrium geometry of the neutral state. In addition, for the sake of simplicity and for the purpose of illustration, Fig. 3 was obtained with the constraint that the dimer always remained at initial symmetry at all times during the reaction. In reality, this of course cannot be true. For more realistic simulations, we thus lifted such constraints. With the typical reaction trajectory in mind, we subsequently carried out an on-the-fly nonadiabatic molecular dynamics simulation with a swarm of 300 nuclear trajectories using the surface hopping including arbitrary couplings (SHARC) package 45,46 . The initial conditions of the trajectories were sampled with the Wigner distribution 47,48 , which allowed setting of the initial coordinate around the equilibrium geometry with an initial distribution of velocity according to the oscillation frequency around the equilibrium geometry of the neutral dimer. Nuclear coordinates were then evolved along the CASSCF-calculated energy profile with a time step of 0.5 fs during 1 ps SHARC simulations. To solve the over-coherence problem, we applied the local diabatization algorithm and energy-based decoherence correction during propagation of the wave function 49,50 . Then, we collected the trajectories that led to the formation of D 3 + and analysed their corresponding formation times.

Data availability
The data that support the plots within this article are available from the Zenodo database at https://doi.org/10.5281/zenodo.7807812.