The kinetic energy of PAH dication and trication dissociation determined by recoil-frame covariance map imaging

We investigated the dissociation of dications and trications of three polycyclic aromatic hydrocarbons (PAHs), fluorene, phenanthrene, and pyrene. PAHs are a family of molecules ubiquitous in space and involved in much of the chemistry of the interstellar medium. In our experiments, ions are formed by interaction with 30.3 nm extreme ultraviolet (XUV) photons, and their velocity map images are recorded using a PImMS2 multi-mass imaging sensor. Application of recoil-frame covariance analysis allows the total kinetic energy release (TKER) associated with multiple fragmentation channels to be determined to high precision, ranging 1.94–2.60 eV and 2.95–5.29 eV for the dications and trications, respectively. Experimental measurements are supported by Born–Oppenheimer molecular dynamics (BOMD) simulations.


.1 General computational details
Computation of the dissociation mechanisms of the dication and trication of the fluorene was performed using an augmented Born-Oppenheimer molecular dynamics (aBOMD) approach, [1][2][3][4] where the motion of the nuclei is governed by the ground electronic state gradients, and the electronic excitation is represented via an internal excess energy (IEE), external energy reservoir. The IEE energy is assigned to the molecule using a recently proposed semi-heuristic model. [4] Following IEE assignment, the energy is redistributed into the vibrational degrees of freedom with a rate given by the internal conversion (IC) constant, which was computed using a classical electron-nuclei collision model. [4] The potential energy surfaces were computed through the GFN2-xTB method [5] from the XTB semi-empirical package, [6] which has been shown to work well for molecular fragmentation dynamics. [3,4] A detailed description of the computational methods used can be found in the PyRAMD manual, [7] and the source code is also available for download. [8] Fragmentation dynamics were computed for fluorene with a time step of 1 fs covering a total time of 10 ps in each BOMD trajectory. The XUV or XUV and IR pulses were applied simultaneously at the first step of the simulation.

Initial conditions sampling
In PyRAMD, we use an intermediate initial conditions sampling variant: beyond standard classical Maxwell-Boltzmann sampling, but simpler than standard Wigner sampling from the vibrational wavefunctions in the harmonic oscillator approximation. We start from the assumption that for each degree of freedom of each nucleus at temperature T = 0 K, we need to fulfill the Heisenberg inequality for the position and momentum: where ∆x 0 and ∆p 0 are the uncertainties for the position and momentum of the particle at T = 0 K. This can be achieved by sampling each coordinate (x) and momentum (p) from two Gaussian distributions: where x 0 is the initial coordinate value at a given XYZ structure, and we assume that initially each particle is motionless. This is equivalent to sampling the coordinates and momentums from the Wigner function: for the Gaussian wavepacket localized around x 0 : that is: To fulfill the Heisenberg inequality, we take the uncertainties parameterized as the following: and ∆p 0 = hm 2τ where m is the mass of the nuclei, τ is the arbitrarily chosen time, and α ≥ 1 is an arbitrarily chosen scaling coefficient. If α = 1 the inequality turns into equality. These parameterizations correspond to the minimal spreading Gaussian wave packet for the time τ . For a hydrogen atom (m = 1 a.m.u.) at τ = 1 fs and α = 1, this yields ∆x 0 = 0.06Å.
We approximate the effect of temperature on the momentum of the particle by adding extra momentum p ′ (p → p + p ′ ) distributed according to a Maxwell-Boltzmann distribution (p ′ ∼ exp − βp ′2 2m , where β −1 = k B T ). This leads to a momentum distribution: After integration over p ′ , we get the distribution: where ∆p 2 T = ∆p 2 0 + mk B T , i.e. instead of pure Maxwell-Boltzmann sampling, we have a sampling with adjusted temperature T eff = ∆p 2 0 /(mk B ) + T . By setting ∆x 0 = 0 and ∆p 0 = 0, we obtain the classical Maxwell-Boltzmann sampling routine.
In our the simulations, we applied the default PyRAMD values of τ = 10 fs, α = 1, and T = 0 K.

Ionization and excitation procedure
Upon interaction with the laser pulses, we assume that the molecule M absorbs N ph ≥ 0 photons γ, emits N els ≥ 0 electrons (e − ), resulting in an excited ionic state M N els + * . Such interaction can be described via a pseudo-chemical reaction: The excitation energy of the M N els + * with respect to the ground ionic electronic state M N els + is the previously described internal excess energy (IEE). The integer numbers N ph and N els are sampled on the fly from the Poisson-distribution, where the chance of getting the result n = 0, 1, 2, . . . is given as p n = λ·exp(−λ) n! with parameter λ = ⟨n⟩ being the mean n. The excitation energy obtained by the molecule from the laser is thus given as E exc = N ph · hν, where ν is the laser pulse frequency.
The parameters λ for the N els and N ph in the simulations were the following: • For the XUV pulse: λ ph = ⟨N ph ⟩ = 1, λ els = ⟨N els ⟩ = 2, ν = 40.9 eV.
The IEE assignment algorithm operates as follows: • Two main parameters are generated from the Poisson distribution: excitation energy provided by the laser pulse E exc = N ph · hν, and N els . The parameter IPStretching is provided by the user. This parameter gives a threshold for the molecular excitation, by how much it can exceed the IP without emitting the electrons (i.e. if E exc ≤ IP · (1 + IPStretching) with IP being the ionization potential, and N els = 0, the molecule stays in the neutral state). We used IPStretching=0.2.
• The IPs and HOMO-LUMO gaps of the molecule are calculated to find the maximal ionization that can be reached, applying the restrictions that the number of removed electrons n cannot exceed N els , and the required energy for ionization of n electrons cannot exceed E exc . If N els = 0 but the E exc exceeds the threshold IP · (1 + IPStretching), then one electron will be removed.
• In the final ionic state, the IEE is assigned by the following rules: -If no electrons were removed (n = 0), the excitation will increase the IEE to the value by E exc . However, if E exc is smaller than the HOMO-LUMO gap of the molecule, the additional IEE = 0.
-If the ionization occurs, the IEE will be drawn from the beta-distribution, such as the probability density for getting the IEE = IEE max · x (x ∈ [0; 1]) given by: where B is the beta-function, D is the scaling law for the density of states (taken as D = 3), and IEE max = E exc − EI, where EI (energy of ionization) is the sum of IPs to remove n electrons. [4]

On-the-fly internal conversion rate estimation
The heuristic model for the internal relaxation extends a model introduced in previous work ( [4]): • electrons and nuclei are modelled as two gases, • the electrons are modelled as hot gas and the nuclei are modelled as cold gas, • the excess energy of electrons (IEE) is converted to the nuclei.
The rate of the relaxation, k, is given by the equation for the collision period τ (see e.g. [9]) where m e and m n are the masses of the electrons and the nucleus, κ is a dimensionless parameter, and z is given by equation: with n e being the concentration of the electrons, σ en is the collisional cross section, and T is the temperature of the electronic gas. The internal excess energy (IEE) is taken to be proportional to temperature (IEE ∝ T ). From here, a set of heuristic assumptions are applied. First, the total rate of IC is taken as k IC = n 1 τn , where n enumerates nuclei, and the τ n is the relaxation time for the n-th nucleus. In each τ n , we assume the cross sections, σ, to be proportional to the nuclear mass as σ n = σ 0 m n /m amu , with σ 0 representing shared cross section per unit of mass, and m amu is the atomic mass unit. The concentration of electrons we will take as n e = N e /V mol with N e being the total number of elections, and V mol being the molecular volume. The total volume of the molecule is approximated as: where Z n is the charge of the n-th nucleus in atomic units, R nm is the distance between nuclei, and C nm is the connectivity matrix of the nuclei with: where R n and R m are the covalent radii of the atoms. Physically, this approximation represents the molecular volume as a set of cylinders located along the bonds, the lengths of the cylinder as the bond distance, and the cross section is σ 0 . By taking the ratio of masses to be memn (me+mn) 2 ≈ me mn (since m n ≫ m e ), we get the following expression for the IC rate: where N n is the total number of nuclei in the system. To stabilize this equation if no bonds are found (i.e. all C nm = 0), we include a stabilizing term to the denominator (L 0 ). This yields the final equation: L 0 was varied manually, whilst κ was optimized using a training set of data from literature values [10,11] for the lifetimes of short lived cations obtained using XUV pulses with the following reaction: An acceptable result was found for L 0 = 10 [Bohr] ≈ 5[Å] and κ = 1.2815616845. The error-weighted RMSD for the training set was calculated to be 11.3 fs, and 11.4 fs for the test set. (see Figure 1). More details on the fitting Figure 1: Comparison of the experimental and predicted internal conversion rate lifetimes for the training and test sets of the molecules. procedure can be found in the PyRAMD repository: https://stash.desy.de/projects/PYRAMD/repos/pyramd/ browse/addData/FitDynamicalStuff/NewGasCollisionModel.

Calculation of the Kinetic Energy Release
The fragmentation procedure used in PyRAMD is described in Ref. [4]. Here, we outline the procedure for calculating the kinetic energy release (KER). Let us assume that the molecule AB dissociates into fragments A and B (AB → A + B). We calculate the relative velocity of the fragments in the direction that connects their centers of mass , where E is the electronic energy of the particle in the parentheses. The resulting kinetic energy release (KER) along the dissociation direction is thus given as KER = KE − DE (see Figure  2), and if KER > 0, then the fragmentation is possible. Upon dissociation, the momenta of the fragments are recorded and used to compute the KER. The KER energy will be distributed between fragments A and B in accordance with the momentum and energy conservation laws. The momentum conservation gives us p A + p B = 0, where "p" is the momentum of the fragment along the n AB axis. Kinetic energy conservation = KER, and therefore the KER velocities assigned to the fragments will be  Theoretical values were calculated for the dication produced by the XUV pulse only, for the dication with both the XUV and IR pulses, and for the trication obtained with both XUV and IR following the loss of a H + or H + 2 ion, i.e. a (1,1,1) channel.
We performed four types of simulation of the dissociation dynamics for fluorene (C 13 H 10 ): produced by only the XUV pulse (10927 trajectories), produced by XUV and IR pulses (1835 trajectories); • C 13 H 10 3+ produced by the XUV and IR pulses (5773 trajectories), • C 13 T 10 3+ produced by the XUV and IR pulses (13783 trajectories).
The main dissociation dynamics are driven be absorption of an XUV photon, with the IR pulse adding only a small level of additional energy compared to the XUV only simulations. The simulation of the dication with and without IR pulse shows that the kinetic energy release (KER) is only slightly changed with the addition of the IR pulse ( Figure  3). The dissociation dynamics of the dication was found to mostly proceed via (1,1) channels with a breakup of the carbon-carbon bonds. The KER in these channels seems to be rather constant at around 1 eV. This is lower than the experimental observation, which is attributed to three main reasons: • In BOMD, the dissociation occurs only on the ground state potential, whereas in reality the electronically excited states will also produce the (1,1) channel ions, and thus have a higher KER.
• We use the semi-empirical approach to compute the highly energetic conformations of the molecules, which can give some systematic errors in the potential.
• Dissociation potential energy surfaces of the highly-charged cationic states can have a dissociation barrier, which might be missed in simulations due to earlier detection of the dissociation (to avoid the multi-reference character of the bond breaking). Therefore, the results might underestimate KER values.
In case of fluorene interacting with XUV pulse to produce the dissocation from the dication, the distribution and the average kinetic energies (KE) of the individual C n H + x fragments are plotted in Figure 4. It can be seen that the higher mass fragments generally have lower KEs than lower mass fragments, as expected.
In the case of the fluorene trication (C 13 H 10 3+ ), upon interaction with the XUV pulse, the only fragmentations observed within 10 ps were the loss of one or two hydrogen ions. Therefore, to enhance the carbon-carbon fragmentation, the IR pulse was added to the simulation. Carbon backbone fragmentation of the trication mostly progressed through the (1,1,1) channel: first, the C 13 H 10 3+ lost either H + or H + 2 , resulting in dicationic C 13 H 9 2+ and C 13 H 8 2+ , respectively, and the carbon backbone fragmentation continued via a (1,1) channel ( Figure 3). To quench this process, the protons were replaced by heavier tritium atoms. Of the 13783 C 13 T 10 3+ trajectories, the following three (1,2) channels were observerd: • C 13 H 3+ 10 → C 7 H + 6 + C 6 H 2+ 4 (KER = 3.6 eV), • C 13 T 3+ 10 → C 7 T 2+ 6 + C 6 T + 4 (KER = 3.1 eV), • C 13 T 3+ 10 → C 11 T 2+ 8 + C 2 T + 2 (KER = 2.1 eV). The first 400 fs of these aBOMD trajectories are also given as movies in the electronic supplementary information. Similar to the (1,1) channels, the (1,2) channels are approximately 1-2 eV lower than the experimental values, which is attributed to the same reasons described above.