Approaching black-box calculations of pump-probe fragmentation dynamics of polyatomic molecules


 A general framework for the simulation of ultrafast pump-probe time resolved experiments based on Born-Oppenheimer molecular dynamics (BOMD) is presented. Interaction of the molecular species with a laser is treated by a simple maximum entropy distribution of the excited state occupancies. The latter decay of the electronic excitation into the vibrations is based on an on-the-fly estimation of the rate of the internal conversion, while the energy is distributed in a thermostat-like fashion. The approach was tested by reproducing the results of previous femtosecond studies on ethylene, naphthalene and new results for phenanthrene.


Introduction
The pump-probe technique is a key approach to understand ultrafast chemistry. It is used for elucidating reaction mechanisms and tracing transition states along with short-lived species. [1,2,3,4,5,6,7,8,9,10] In this technique one pulse (the pump) excites some inter-/intra-molecular dynamics, and a second delayed pulse (the probe) changes or visualizes the transient dynamics of the molecule, allowing for a scan of the intermediate states. Depending on the pump, the probe, and the observed signal, there is a large variety of different experimental methods to study the induced ultrafast dynamics in molecules: they differ in how the induced pump-probe delay signal is observed. For example, the response can be measured with ion/electron spectroscopy [3] by measuring the yields of different electrons and fragments from the molecule; absorption spectroscopy that measures the amount of absorbed photons of certain kind; [4,11,12,13] ultrafast diffraction/microscopy, which investigates the changes in scattering patterns, [8,9,10] etc. An understanding of these types of experiments can be obtained through the framework of explicitly simulating the underlying processes. [14,15,16] In this paper, we focus on a theoretical description of pump-probe experiments used to investigate the fragmentation dynamics of molecules. The most straightforward way is to simulate femtosecond chemical responses based on the quantum mechanical treatment of the nuclear motions in several electronic states (beyond the Born-Oppenheimer approximation). Examples of numerical approaches of this sort are wavepacket dynamics, [17] multiconfiguration time-dependent Hartree (MCTDH) [18,19,20,21,17] and ab initio multiple spawning (AIMS) [22,23,24] methods. However, performing quantum dynamics of this sort is computationally very demanding, since their computational time scales exponentially with the number of degrees of freedom. [17] Thus, these types of simulations are typically applied to full simulations of few-atom molecular systems, [17,22,23,24] or they do not account for all of the degrees of freedom, only a few selected ones. [19,20] By neglecting the quantum nature of nuclei, it becomes possible to reduce the cost of simulations. [25,26,27,28] This leads to molecular dynamic (MD) approaches. For motions in multiple electronic states, non-adiabatic MD methods must be applied, for instance trajectory surface hopping dynamics. [29,30,26,27,28] These types of simulations can be applied to polyatomic molecules in a straightforward manner: in fact, molecular dynamics based simulations have already been successfully used for simulating the pump-probe yields of various molecular fragments in polyatomic systems. [28,31] However, non-adiabatic MD methods are more computationally demanding than the Born-Oppenheimer MD (BOMD) method, which is performed only on the ground state potential energy surface.

The general structure of the simulation
The idea behind the presented simulation for mass spectra in pump-probe experiments was adapted from the EIMS simulations within the QCEIMS method. [32,33,34,35] However, unlike in EIMS, in this pump-probe scheme there are two subsequent excitations of the system instead of one. Therefore, the algorithm for such modeling has the following structure. 1. Generation of initial conditions (coordinates and velocities of the atoms) that represent the initial ensemble of the molecules. In principle, this can be done in the same way as in QCEIMS [35] by imitating a hot molecular vapor in the mass-spectrometer chamber with BOMD simulations. In this case, the vibration of the initial molecule in the electronic ground state is calculated  with canonical ensemble ( ) MD at elevated temperatures (for example, at ∼ 500 K [35]), followed by choosing a set of random time frames within the trajectory as the initial conditions. In principle, this procedure samples the full conformational space of the molecule. In contrast to EIMS, which deals with any vaporizable substance, the molecular species in pump-probe experiments are often chosen to have a limited conformational space; therefore, performing the full electronic ground state sampling is a redundant operation. In addition, the temperatures of the molecular ensemble in these experiments are often cooled down to a few Kelvins using, for instance, the adiabatic expansion technique, which makes the classical MD at finite temperatures an unrealistic model.
In the examples we considered in this work (ethylene, naphthalene, and phenanthrene), there is only one conformer, and most of the vibrations have frequencies > 800 cm −1 , which correspond to Debye temperatures D > 1200 K ( D = ℎ B , where ℎ and B are the Planck's and Boltzmann's constant, respectively). This means that most of the motions of our molecules at the experimental temperatures are in the vibrational ground state, and the classical Maxwell-Boltzmann distribution is not a suitable representation of the molecular ensemble compared with the quantum distribution at = 0 K. Based on these assumptions, we have generated initial conditions from the Wigner distribution [36] obtained from a harmonic oscillator approximation of the molecule at an equilibrium geometry. This ensemble of initial conditions obtained from the ground state represents the molecules at = 0 K.
2. BOMD simulations are performed for each of the initial conditions. The overall time duration for each trajectory is PP + free , where PP is the pump-probe delay and free is the additional time for the molecule to fragment, around 2-6 ps in our cases. The electronic excitation is modeled via an additional degree of freedom, namely the internal excess energy (IEE ≥ 0), as introduced in QCEIMS. This energy is constantly drained into the intramolecular vibrations via the first order differential equation IEE = − IC IEE, where IC is an internal conversion constant representing the electronic excitation decay via non-adiabatic processes, such as the passing of conical intersections and avoided crossings. 3. The interaction with the pump and probe pulses brings about changes in the state of the molecule. The electronic excitations that occur without ionization can be emulated as an increase of the IEE by a pre-set value or as a change in the electronic state. The ionization by an XUV pulse, in addition to previous changes, increases the charge state of the molecule. The first pulse comes at the first step of the trajectory, and the second acts after time PP has passed. 4. Fragmentation of the molecule during BOMD is detected based on the analysis of a molecular graph constructed from interatomic distances by a standard criterion of bonds as ≤ · ( + ), where is the interatomic distance between atoms and , is the atomic radius, taken in our case from Ref. [37], and is a scale factor between 1 and 2 (in our case = 1. , where = A ( +) , B ( +) and e is the number of electrons. 5. Finally, the number of the fragments formed from different initial conditions is counted to yield the mass spectrum for the particular pump-probe delay PP between the lasers.
This sequence of actions has to be taken at different pump-probe delays, PP , yielding a two-dimensional mass spectrum MS ( / , PP ), which makes the simulation scheme more demanding in computational resources when compared to EIMS. To obtain a better consistency with the experimental data, the calculated fragmentation yields MS ( / , ) can be convoluted with the instrumental func-  ) along the coordinate PP , where cc is the cross-correlation time consisting of the widths of the pump and the probe pulses.

Internal excess energy after ionization
QCEIMS assigns the internal excess energy (IEE) using the Poisson-type distribution with two adjustable parameters [35] that describe the process of electron impact ionization However, a straightforward application of the proposed distribution to the laser excitation described by reaction is not obvious. Therefore, we have derived an oversimplified analytical model to estimate the IEE after a generalized ionization reaction of the molecule (see Fig. 2) The processes given by Eqs. 1 and 2 are particular examples of such reaction. The IEE is basically the energy difference between the electronic energy of the final state of M ( +) and its ground state electronic energy. Here we assume that the ionizing radiation acts with an energy of r, which is equal to r = ℎ in the case of ionization by photons with frequency and to r = e 2 2 in the case of electronic impact ionization ( e is the mass of the electron and is its velocity). The formula for the photoelectric effect [38] in this case is where IP + IEE corresponds to the work function of the molecule (IP is the ionization potential, and the IEE is an electronic excitation) and T e, is the kinetic energy of the ′ electrons in the right-hand side of ionization reaction given by Eq. 3.
If we assume that there is no resonant excitation, then the system can end up in any of the excited states ( > 0) as well as in the ground state ( = 0, IEE 0 = 0) with equal probability. Therefore, we can determine the distribution of occupancies of the molecule's excited states by applying the maximum entropy principle, which is known to be useful for the interpretation of fragmentation patterns. [39,40] The first contribution to the total entropy is the entropy of the ion M ( +) , where = 0 + 1 + . . . + + . . . is the total number of ionized molecules, is the number of ions in the -th electronic state, and B is the Boltzmann constant. [41] The second contribution to entropy is the entropy of the leaving electrons ( e). The number of states for a gas with f degrees of freedom in the vicinity of energy Te is gas ∝ T The requirements r ≥ IP + IEE and IEE being smaller than the ionization potential of M ( +) should be fulfilled for the application of this distribution. The f of the electrons is equal to 3 in the case of photoionization ( ′ = ) and to 3( + 1) in the case of electron impact ( ′ = + 1), due to the fact that the ionizing electron also leaves the molecule. The distribution given by Eq. 6 in contrast to the one in QCEIMS does not have any free parameters, and it only requires the energies of the excited states of the ion M ( +) .
A better understanding of the resulting distribution can be found if we consider a simple toy model. Let us assume that electronically excited states of the molecule can be described via the density of states function DoS(IEE) ∝ IEE ( −1) , and the maximum excitation energy IEEmax = r − IP is smaller than the ionization potential of ion M ( +) . Then the discrete distribution obtained from Eq. 6 can be reduced to a continuum probability density for obtaining an IEE, which then becomes a beta-distribution:

Internal excess energy after excitation
In this section we present the model for a non-ionizing electronic excitation of the molecule by a laser, which can be described by the reaction M ( +) + · ℎ → (M ( +) ) * . For simplicity we assume that the excitation is caused by a laser pulse that is off-resonant with both vibrational and electronic degrees of freedom: the energy of the photon should be larger than the vibrational frequencies and smaller than the energy gap between ground and first excited electronic states, which is the case of 800-nm laser excitations. To simplify even further, we will model the laser pulse as a two level system of many independent particles: |photon⟩ with the energy of a single photon (ℎ ) and |no photon⟩ with energy 0. The molecular entropy is given by the same expression given by Eq. 5. Ignoring the possible coherencies between the photons, we can estimate the entropy of a laser field as field where is the number of photons in state |photon⟩ after interaction, ini is the initial number of photons, and thus ( ini − ) = is the amount of photons that were absorbed by the molecule. The overall entropy should be maximized with two conditions: energy is conserved, thus ini ℎ + 0 This leads to the expression for the occupancies of the molecular electronic levels: where = / ini ∈ (0, 1) is the transmittance (portion of the photons that have not interacted with the molecule). This distribution is a Boltzmann distribution with an effective temperature of B eff = ℎ / ln . However, since 1 − 1 can take any value from 0 to +∞, its logarithm can range from −∞ to +∞. If < 0.5 (i.e. most of the photons were absorbed), then eff > 0, and the distribution is Boltzmann-like. In the case of > 0.5, eff < 0 leads to the non-thermal distribution with inverted Boltzmann populations. At = 0.5 the distribution becomes uniform over all the possible electronic states. Illustration of these features of Eq. 7 is given in Fig. 3. Since there is no obvious way to obtain , we leave it as a free parameter of the model that can be used to fine tune the simulations.

Internal conversion rate
QCEIMS performs an estimation of the internal conversion (draining of the electronic excitation into the molecular vibrations) rate using an energy-gap law with two parameters. [35] Such treatment of internal conversion assumes that the molecule has multiple channels allowing fast electronic de-excitation. In this work we have applied another model for the internal conversion rates based on the representation of the internal conversion as a hot electronic gas heating the cold nuclei. Therefore, the internal conversion in this case is the excess energy from electronic excitation transfered by collisions. For describing this process we will consider two colliding classical particles: one of them (electron) has a mass of e and energy e 2 0 2 = IEE/3, the second one (nucleus) has a mass and is motionless before the collision. We will assume the collision to be elastic, and therefore, after the application of energy and momentum conservation laws, we will obtain the velocity of the electron after the collision to be 1 = 0 · e− e+ . Thus the kinetic energy loss of the electron is If the molecule contains multiple nuclei with masses 1 , 2 , . . ., and the electrons collide with the -th nuclei with frequency , the decay of the IEE can be described with a first-order differential equation: where IC is the internal conversion rate. An approximation for the collision frequency of the electrons with the nucleus can be made by the so-called electron plasma (or Langmuir) frequency: [43,44]

Uncertainties for yields in different fragmentation channels
The ab initio MD simulations are computationally and memory demanding, which leads to low statistics of the results. Therefore, in the prediction of ratios for different fragmentation channels an important issue that arises is that of calculating the confidence intervals. In our simulations we have a molecule M in certain excitation conditions. The question is: what is the probability that fragment A will be produced in a fragmentation reaction M → A + . . .? A simple answer can be obtained using a binomial distribution assumption where we have two possible outcomes for this channel: 1) "A was produced" and 2) "A was not produced". If we obtained trajectories and A appeared in of them, then the naive estimation of the probability of fragmentation yielding A corresponds to = with an uncertainty where is the upper 1 − /2 quantile of the standard normal distribution [45,46] with = 1 for an uncertainty corresponding to 1 . However, this scheme gives adequate estimations only in the limit of large ; therefore, for small data sets we have used a modified version of this method: [45,46]

General considerations
The simulations were performed using a home-written molecular dynamics code in Python (see SI) that relies on the XTB software [47,48] as the source of potential energy surfaces at the GFN2-xTB [48] level of theory. This semi-empirical theory has been shown to perform well in the simulation of EIMS spectra. [49] Mass spectra at each pump-probe delay for ethylene and phenanthrene were computed using 500 trajectories, whilst for naphthalene 400 trajectories were used. The integration step for BOMD in all simulations was 0.5 fs, and the free evolution time free after the pumping and probing was 2 ps for ethylene, 6 ps for napthalene, and 2.5 ps for phenanthrene. The ionization potentials and excited state levels for ethylene and naphthalene were calculated using DFT/sTD-DFT [50] with the Orca 4 program suite [51] at the PBE0/def2-TZVPP [52,53] level of theory with the RIJCOSX approximation [54,55] utilizing the def2/J auxiliary basis sets. [56] The molecular geometry was optimized at the ground charge state, followed by single point energy calculations and sTD-DFT excited state calculations of different charge states at their frozen geometry. Analogous computations for phenanthrene were done using the B97-D functional [57] and canonical TD-DFT calculations. The ensemble of initial conditions was calculated using a Python script from the DFTBaby [58] package based on GFN2-xTB harmonic force fields. All the calculations were enabled by Maxwell computational resources operated at Deutsches Elektronen-Synchrotron (DESY), Hamburg, Germany.

Ethylene
To initially test this approach we chose ethylene (C 2 H 4 ). For this system, we reproduced the data collected in Ref. [24] in which C 2 H 4 was investigated with a 160 nm UV laser (pump) and an XUV probe pulse with photon energy of 17-23 eV. That study involved simulating pump-probe yields using the AIMS method based on SA-2-CAS(2/2)-MSPT2/6-31G* (state-averaged complete active space self-consistent field wavefunction with an active space of two electrons and two orbitals using the multi-state variant of second-order perturbation theory) quantum chemistry calculations. [22] In our simulations we assumed the following conditions. -If the probe (XUV) pulse hits the molecule first, it ionizes the molecule from the ground (singlet) state and produces an average IEE of 1.4 eV. Since the pump pulse is not in resonance with possible electronic transitions anymore, it does not contribute to the IEE. -If the pump (UV) hits first, it performs a resonant electronic excitation → * promoting an electron from HOMO to LUMO and increasing the IEE by ℎ − ( (C 2 H 4 * ) − (C 2 H 4 )) = 7.7eV − 5.8eV = 1.9 eV, where the excitation energy of 5.8 eV was obtained using sTD-DFT. The second pulse (probe) can access more states starting from a higher lying C 2 H 4 * state; therefore, the leftover IEE after ionization will be 8.8 eV.
Upon first glance the resulting pump-probe ion yields from the BOMD simulations (Fig. 4, panel b) are similar to both the experimental values (Fig. 4, panel a) and the ones computed with AIMS (Fig. 4, panel c). First of all, the two non- trivial fragmentation products (CH 3 + and H 2 + ) appear in our MD simulations.
The parent C 2 H 4 + ion (green line) shows the same "S"-shaped decay as in the experiment, CH 2 + (red) and CH 3 + (blue) show the similar transient behavior in all three cases with peaks around PP = 50 fs and PP = 100 fs, respectively, and H + and H 2 + are increasing by a similar "S"-shaped curve.
However, there are some discrepancies that need to be explained. The yield of the parent in our simulation goes to zero, whereas in the experiment it only decreases to 90%. This difference is caused by the fact that we model only a single channel with a photochemical yield of 100%, which is not what happens in the experiment. The other discrepancies, i.e., lower maximum yield of the CH 2 + (2% in our case vs. 10% in AIMS) and the absence of the delayed formation of H 2 + with respect to H + observed both in the experiment and AIMS simulations might be caused by the quality of the quantum chemical methods producing the PES. The absence of nuclear quantum effects [59] in our BOMD simulations might also contribute to the observed differences. Nevertheless, taking into account the exceptionally small computational demands of our simulations in comparison to AIMS, the agreement of our simulations both with the experimental data and more sophisticated theory (see Introduction) is promising and gives a good first (and fast) approximation for molecular movies.

Naphthalene
In order to further test our model, we concentrated on simulating two-color pump probe experiments of naphthalene carried out using an XUV pump pulse with photon energies of 20.4, 23.6, 26.7, and 29.8 eV, generated using an HHG source and an IR probe pulse of 800 nm ( [60]). We assumed the lowest XUV photon energy (20.4 eV), which gave an average IEE of C 10 H 8 + after ionization to be 8.7 eV. The probe (IR) laser was heating the system as C 10 H 8 + ℎ IR → C 10 H 8 * → C 10 H 8hot (multiphotonic electronic excitation followed by a decay of this excitation into vibrations); therefore, a model from section 2.3 was applied. The parameter was taken to be 0.5, i.e., we assumed that the IR laser uniformly shuffled the populations of the excited states of naphthalene. In the case when the probe laser was acting before the pump, the average IEE from the distribution given by Eq. 7 (for C 10 H 8 ) was 6.5 eV, and in the case of probing after pumping, the IEE was   The overall yields of the charged fragments decrease with a decreasing of their mass in both theory and experiment. This effect is basically determined by the ratios of the ionization potentials (IPs) of fragments C n1 H x1 and C n2 H x2 , which determine the dissociation energies along two possible ionic breakup processes: Ionization potentials tend to increase with the decrease of the system size, [61] therefore the smaller the pair fragment C n2 H x2 gets, the less energetically favored the second channel becomes. This leads to small fragments being formed due to a long chain of events compared to large fragments that appear at small fragmentation generations (by generation we mean the number of unimolecular dissociation events that have been undertaken to form the desired fragment). This can be directly observed from the simulations (see Fig. 6). However, the second channel, leading to shorter carbon chains, also remains possible with non neglible yield. In the current software implementation the algorithm always chooses the lowest energy dissociation pathway. This most likely causes the deviations in the trend of the ion intensities in our simulations.
In the experiment, in addition to a "S"-shaped pump-probe behavior, there is an additional transient peak, which does not appear in our simulations. This is due to the absence of any of the transient channels in our simulations. An example of such channels can be where (C 10 H 8 + ) * diss denotes the dissociative excited state of the naphthalene cation, and X + = CnHm + is one of the observed fragments, or i.e. ionization of short-lived excited fragments, that are close to the ionization threshold. However, we can reproduce some of this transient behavior directly from our simulation. The dication (C 10 H 2+ 8 ) formation given in the panel (b) of Fig. 5 is mostly given by a pathway as proposed in Ref. [19] The internal conversion time ic of (C 10 H 8 1+ ) * estimated during simulation is 47 fs, which is relatively close to the values obtained from experiment (27 ± 2 fs, 37 ± 3 fs, 37 ± 2 fs, 43 ± 5 fs at XUV photon energies of 20.4, 23.6, 26.7 and 29.8 eV). [60] The ion yield of C 10 H 8 2+ produced by the scheme given in Eq. 9 is described via the analytical expression [62] ( where cc is the cross-correlation parameter from the instrumental function of the form exp(−( / cc) 2 ). This expression is plotted in the panel (a) of Fig. 5 with the experimentally determined cc = 56 fs. [62] Overall, the experimental and theoretical behaviors are similar. Thus, while the current model for interaction of the molecule with the IR laser does not include ionization, this posterior inclusion of transient channels looks promising.

Phenanthrene
The last test case is that of a phenanthrene (C 14 H 10 ) pump-probe study [B.
Manschwetus et al., manuscript in preparation] conducted at the CAMP end station of the Free Electron Laser At Hamburg (FLASH). [63] In this case, a twocolor experiment was done using XUV radiation with photon energy of 41 eV from FLASH and IR radiation produced by a 810 nm Ti:Sa laser. The time 0 , when the XUV and IR lasers were hitting the molecule simultaneously, and cross-correlation time cc = 130 fs were determined from electronic velocity map images (VMI) based on the helium fluorescence line and the first helium sideband. For the MD modeling we have chosen the following excited dication formation channels: where "*" and "**" denote single and double pulse excited molecular species. The MD simulations for phenanthrene were done with the IEE after XUV excitation being 14 eV, and the IR pulse before/after XUV added 7/13 eV into IEE, respectively. Positive pump-probe delays in both the theory and experiment represent the XUV acting before the IR.
First of all, it is interesting to compare the yields of the dication itself. The theoretical curve is "S'-shaped, similar to the case of naphthalene. The data obtained can be compared with the stepwise pump/probe vs. probe/pump behavior contribution in the experimental data. This switch can be described by the expression [62] ( ) = 1 + erf As it can be seen in Fig. 7, our theoretical curve and this particular experimental contribution do match. However, the experimental curve of the phenanthrene dication (C 14 H 10 2+ ) has also a prominent transient peak, which is produced by a channel similar to the one discussed in the case of naphthalene: [19] and also by a similar channel in which the IR laser acts as a pump: Both of these channels are absent in our simulations. Nevertheless, we can introduce them in the same way as we did for naphthalene, namely by calculating the internal conversion times of (C 14 H 10 1+ ) * and (C 14 H 10 ) * . According to GFN2-xTB-based estimations they are 36.7 and 35.7 fs. We include them as [64] Fig. 7: Comparison of the experimental and theoretical pump-probe yields for various channels in the case of phenanthrene (C 14 H 10 ). The "MD" curve represents yields obtained directly from simulation, "MD + transients" is the "MD" curve augmented with an analytical curve given by Eq. 12, "Experimental step function" corresponds to Eq. 11 fitted to the experimental data, "Experiment" represents the experimental data.
where (C14H10 1+ ) * and (C14H10) * are given by Eq. 10 and is a fitting constant. Adding this expression to the C 14 H 10 2+ yield from MD simulation and adjusting the value, we can get an almost perfect agreement with the experiment (curve "MD + transient" in Fig. 7). From the covariance analysis of ionic channels [65,66] we know that a prominent decay channel for the dication is a Coulomb explosion C 14 H 2+ 10 → CnH + x + C 14−n H + 10−x with = 2, 3, 4. Therefore it is interesting to compare the pump-probe yields of these channels in the experiment with our theoretical predictions. The high-mass fragments in the experimental data will be significantly contaminated by the products of the monocation, therefore we focus on smaller fragments (see Fig. 7). As in the case of the dication, the stepwise behavior is reproduced correctly, while the transient increase is not captured. The origin of this increase in the case of C 2 H + x , C 3 H + x and C 4 H + x is not as obvious as in the dication, so that the transient behavior is not treated a posteriori. Nevertheless, some important features of the pump-probe experiment of phenanthrene could be reproduced by our approach.

Drawbacks of the current implementation
The current implementation of the proposed simulation method uses several assumptions and simplifications that were discussed in the previous sections. Here we will address them in more detail and also point out possible improvements to be implemented in a stepwise manner in the near future.
-The values of the IEE produced by laser pulses are pre-assigned before the simulation is started. These values are the expectation values from distributions given by Eqs. 6 and 7. This limitation is crucial from two perspectives. If the molecule is dissociated/isomerized after the first pulse, it will have a different excited state manifold, thus the IEE distribution after interaction with the second pulse will be different. Therefore, a pre-assigned average value will not match with the molecular system in the simulation. The second problem is the influence of the second, third and so forth moments of IEE distributions, which can have a drastic influence on ionic channel ratios. This assumption can be corrected by introducing on-the-fly calculations of the ionization potentials and excited state manifolds at specific trajectory moments. -If several fragmentation channels AB q+ → A (q−n)+ + B (n)+ with different n are possible, the current software always chooses the one with the lowest dissociation energy. However, the least populated channel can also have a significant probability of happening, therefore in our simulation we have biased the results with respect to dissociation probabilities. To improve this behavior we can apply a concept of statistical charges as used in QCEIMS. [35] -The current simulation version relies on the pre-assigned results of the molecule -laser interactions, i.e. the number of electrons removed by the XUV laser is defined beforehand, while the IR cannot ionize the molecule. A more favorable situation would be if the ionization state would be determined on-the-fly during the simulation. It would open the possibility to directly simulate the transient behavior similar to that in Eq. 9. However, this extension would again require on-the-fly calculation of ionization potentials and excited states. Also a generalization of Eqs. 6 and 7 has to be derived to include multiple ionization channels. -The last major drawback of our implementation is the inclusion of the dissociative excited states in the simulation. This issue cannot be resolved within BOMD simulations; however, switching from BOMD to the Car-Parinello approach (CPMD) [67] would probably be a good solution and even decrease the MD cost.
Despite these drawbacks we have demonstrated the applicability of the current methodology to simulate pump-probe experiments involving XUV excitation.

Conclusions and perspectives
We have developed a general framework for the calculations of fragmentation yields of ions in pump-probe experiments. The model is based on the ideas from QCEIMS, [35] but we have developed new few-parametric models for calculating the internal excess energy (IEE) in the molecule after ionization and multiphotonic laser heating. In addition, a model for on-the-fly estimation of the internal conversion rate constants describing the relaxation of the electronic excitation into vibrational degrees of freedom was presented. The methodology was applied to three test cases: ethylene, naphthalene, and phenanthrene. In the first case a good agreement was found with both the experiment and high-level quantum dynamics simulations. In the case of naphthalene and phenanthrene we have reproduced the simple overall behavior of switching between pump before probe and probe before pump cases and also introduced some of the transient dynamics using internal conversion model.
In the end we have described the main drawbacks of the current implementation and potential ways for their improvements, such as on-the-fly random generation of IEE rather then usage of pre-calculated average values.