A 99 , 013409 ( 2019 ) Dynamics of interatomic Coulombic decay in neon dimers by XUV-pump – XUV-probe spectroscopy

We apply the Monte Carlo wave-packet approach to study the interatomic Coulombic decay (ICD) dynamics of neon dimers after removing a 2s electron from one of the Ne atoms by one-photon absorption from an XUV pulse. This method reproduces well both the lifetime for the 2s inner-valence vacancy in Ne2+(2s−1) and the kinetic energy release (KER) spectra for the coincident Ne+ and Ne2+ fragments following triple ionization, i.e., two photoionizations and one ionization via ICD, of Ne2 measured in an XUV-pump–XUV-probe experiment [K. Schnorr et al., Phys. Rev. Lett. 111, 093402 (2013)]. Comparisons between the calculated and measured nuclear KER spectra give physical insights in the considered process. For example, an analysis of the ratios between the lowand high-energy peaks in the nuclear KER spectra for large delays provides an estimate of the photoionization cross sections for removing a 2p electron from the excited states in Ne2+(2s−1). Such comparisons also allow an estimate for the ICD rates for the 2s inner-valence vacancy in the single-site two hole state in Ne22+(2s−12p−1). Finally, the influence of photon statistics of the free electron laser pulses on the nuclear KER spectra is considered.


I. INTRODUCTION
Modern light sources such as free electron lasers (FEL) can produce ultrashort and ultrastrong light pulses with very energetic photons.The photon energies can be large enough to induce ejection of an inner-valence or even a core shell electron from an atom via one-photon absorption.Following such ionization, the remaining cations are left in highly excited electronic states.These states may be well above the lowest double ionization threshold and can thus undergo an efficient and ultrafast decay to more stable electronic states by a secondary electron emission.The well-known Auger decay [1] is such an example and it takes place in the following way: A core shell electron is first removed by the external radiation, an outer-valence electron then fills the core hole and the energy released by this transition is sufficient to eject another outer-valence electron from the same atom.Thus the Auger process is intra-atomic.The case is different when an inner-valence electron is initially removed from an atom.Even though the outer-valence electron can still refill the inner-valence vacancy, the released energy is insufficient to remove an outer-valence electron from the same atom or the neighboring atom in the case of a molecule.Thus no electronic decay occurs in this case.It would, however, be quite different if the atom has weakly bound neighbors with which it can form dimers, trimers, and so on.In this case, the released energy may be sufficient to kick out an outer-valence electron from a weakly bound neighbor.This is Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
because the lowest double ionization threshold for the weakly bound clusters is well below the lowest double ionization threshold for atoms and molecules.This process, known as interatomic or intermolecular Coulombic decay (ICD), results in two positively charged ions which repel with each other by the Coulomb force between them.ICD is ubiquitous in weakly bound clusters and it has received a lot of interest because of its capacity in reflecting cluster structures and as a means of energy transport.
The study of the ICD process in weakly bound clusters started by the first theoretical prediction in 1997 [2].Since then, experimental evidence for the ICD process has been obtained in van der Waals clusters [3][4][5][6] and hydrogen bonded clusters [7][8][9].Depending on the clusters bond type and size [10], ICD occurs at a timescale from a few femtoseconds to hundreds of femtoseconds, much faster than typical radiative decays.Theoretically, this ultrafast ICD dynamics can be traced by applying wave-packet propagation [11].For van der Waals clusters, like neon clusters, the ICD dynamics are comparable in timescale with the nuclear dynamics and thus the nuclear degree of freedom should be taken into account.The influence of vibrational motion on the ICD dynamics has been studied both theoretically [12][13][14][15] and experimentally [16,17].In the earlier works, photoelectron spectral linewidths of resonances [2,3] and direct calculation of ICD rates [10,[18][19][20][21][22] were used to obtain the ICD lifetimes.Now as an alternative, it is possible to apply time-resolved tools and pump-probe measurements [23][24][25].It is, however, very challenging for theory to perform quantum computations of both the nuclear and electronic dynamics in pump-probe experiments since the process of interest involves multiple ionization events and associated continua.For instance, in the recent experiment [23], coincident Ne 2+ and Ne + kinetic energy release (KER) spectra at different pump-probe delays were measured to obtain the lifetimes of ICD processes in the neon dimer.In that work, the nuclear KER spectra were only qualitatively reproduced by solving Newton's equation of motion for the nuclear motion.Since then, to the best of our knowledge, no results of the nuclear KER spectra from quantum calculations have been achieved.As a result, the aim of this work is to reproduce the measured KER spectra in Ref. [23] by applying a quantum mechanical method.It was shown in previous works that the Monte Carlo wave-packet (MCWP) approach is especially useful in reproducing the experimental nuclear KER spectra following dissociative double ionization of H 2 [26], D 2 [27], and O 2 [28].In addition, this method was applied to predict nuclear KER spectra in the dissociative ionization of small molecules in the long wavelength regime [29], a regime that gains increasing experimental interest [30], and more recently the method was extended to describe breakup in XUV-IR pump-probe spectroscopy [31].In this work, we extend the MCWP approach to accomplish the task of predicting nuclear KER spectra following triple ionization of Ne 2 to resolve the ICD dynamics in Ne 2 .The comparison between the theoretical and experimental results shows the potential of the MCWP approach in resolving the ICD dynamics in weakly bound clusters.
The paper is organized as follows.Section II provides a brief discussion of the physical processes considered.Section III presents a summary of the implementation of the MCWP approach to treat triple ionization of the neon dimer and includes information about the parameters used.In Sec.IV, we present the calculated nuclear KER spectra for the coincident Ne 2+ and Ne + fragments following triple ionization including an ICD process of Ne 2 within the XUV-pump-XUV-probe setup used experimentally [23].The applied laser pulses are either coherent pulses or noisy pulses modeling the photon statistics of FEL sources.We also provide a comparison between the calculated and measured nuclear KER spectra.Section V concludes.

II. PROCESS
We begin with a brief discussion of the processes to be studied.In the experiment [23], the KER spectra for the coincident Ne + and Ne 2+ fragments following triple ionization of the neon dimer were measured as a function of delay between pump and probe XUV pulses.During triple ionization of the neon dimer, four charge states, i.e., Ne .Such triply charged ions finally result in the repulsive Ne + and Ne 2+ pairs.The measured KER spectra for these pairs are dependent on the ICD process in Ne 2 and thus they were used to extract the ICD lifetime in Ne 2 .In brief, the lifetime of the process denoted ICD 1 in Fig. 1 can be estimated by monitoring the signal in the Ne 2+ (2p −2 )-Ne + (2p −1 ) channel as a function of delay.As will be detailed below, this is possible because the transition from the Ne + (2p −1 )-Ne + (2p −1 ) channel to the Ne 2+ (2p −2 )-Ne + (2p −1 ) channel is nearly independent of internuclear distance R, and therefore of delay.The potential energy curves for the electronic states involved in our simulation are taken from Ref. [32].These were also used in Ref. [23] and are plotted in Fig. 1.The initial state is the electronic and vibrational ground state of Ne 2 .The dipole couplings between the ground state and the other bound states in Ne 2 can be neglected because the XUV photon energy is much larger than the energy differences between these states.After removing an inner-valence electron, i.e., the 2s electron, from one of the neon atoms by the XUV pulse, highly excited repulsive 2 2 + g and attractive 2 2 + u states are reached (Fig. 1).The dipole couplings between these two states can also be neglected due to the large XUV photon energy compared to the energy spacing between these two states.Both states are unstable and decay via ICD to the two-site doubly outer-valence-ionized Ne 2 2+ .These two-site doubly ionized states are energetically very close to each other (some of them are even degenerate).To capture the mean effect of all these states, it is reasonable to involve only one effective state in the simulation (its potential energy curve is pointed to by the arrow labeled as ICD 1 in Fig. 1).Further photoionization can take Ne 2 2+ to the two-site triply ionized Ne 2 3+ by absorption of another photon from the remaining part of the pump pulse or from the delayed probe pulse.The potential energy curves for the three series of energetically allowed two-site triply outer-valence-ionized states are nearly parallel at the internuclear distances where the nuclear wave packets evolve [23].Thus, we can capture the nuclear dynamics in Ne 2 3+ by including only one effective state and by using the total photoionization cross section from Ne 2 2+ (2p −2 ) to Ne 2 3+ (2p −3 ).In fact, there is a probability to open another ICD channel (the ICD 2 channel in Fig. 1): Before ICD from the highly excited 2 2 + g and 2 2 + u states takes place, photoionization by the XUV pulse may remove another outer valance electron from these two states, reaching highly excited states in Ne 2 2+ .These states can also decay to two-site triply ionized Ne 2 3+ via ICD.Here only one state in the excited Ne 2 2+ (2s −1 2p −1 ) is considered for simplicity.We will return to a discussion of the parameters needed for the modeling of the experiment once the methodology has been presented.

III. METHOD
In this work, the MCWP approach is applied to study the dynamics of neon dimers outlined in Sec.II (Fig. 1).Detailed discussions of its implementation, equivalence with the master equation on the Lindblad form, and validation in pump-probe settings are provided in earlier works [26,31,33].This method treats ionization as a decay process by adopting a non-Hermitian Hamiltonian.To describe the decay from a given charge state to its adjacent charge state with one electron less, the Hamiltonian for the former charge state is constructed by adding a non-Hermitian term H d to the system Hamiltonian H s , i.e., (1) The system Hamiltonian H s in Eq. ( 1) is constructed by three parts, i.e., , where T N is the nuclear kinetic energy operator, H el is the electronic Hamiltonian, including the nuclei-nuclei interaction operator, satisfying the time-independent Schrödinger equation with χ m ( R, t ) the nuclear wave packet corresponding to the electronic state |φ el R,m .The wave function of each charge state in Eq. ( 2) is a coherent superposition of bound electronic states and the continuum electronic states form the surroundings of the system.The interaction of the system in a given charge state with its surroundings is embraced in the non-Hermitian term in the Hamiltonian H .By inserting Eq. ( 2) into the TDSE and projecting on φ el R,m |, the evolution of the nuclear wave packet χ m ( R, t ) along the potential energy surface E el m ( R) is obtained.We assume the neon dimers to be rotational frozen since the rotational timescale is much larger than the vibrational timescale.As a result, the direction of the three-dimensional internuclear distance vector R is fixed.The evolution of the nuclear wave packet χ m (R, t ) along its potential curve E el m (R) can be expressed as where nonadiabatic coupling terms have been neglected.Here R denotes the Ne-Ne distance along the linear polarization direction and m Ne 2 is the reduced mass of the neon dimer and the sum over k in the last term in Eq. ( 3) describes dipole couplings between |φ el R,m and the other states in the same charge state.Specifically, the matrix V mk (t ) is given by . The split-operator fast Fourier transform method [34] is applied to numerically solve Eq. ( 3) in each charge state.In our simulation, we use the time step t = 4 a.u. and the space step R = 0.0215 a.u..It seems that t is too large compared to the period of the XUV pulse.It is, however, acceptable since the laser-system interaction considered in this work is dependent on the slow-varying intensity instead of on the fast varying instantaneous field.The simulation box extends from R min = 3.32 a.u. to R max = 91.15a.u..For the nuclear wave packets in different charge states, the connections between them are made by quantum jumps, in other words, ionization events, described by mn (R, t ), as we shall see see below.The results were checked for convergence by changing t, R, R min , R max , and sampling density of quantum jumps.
The application of the MCWP method to treat triple ionization of the neon dimer is as follows.The initial nuclear wave packet is the ground vibrational eigenstate of the ground electronic state of the neon dimer.We first solve Eq. ( 3) for the neutral dimer with the above initial condition at t = 0 to obtain the nuclear wave packets in the neutral system.The probability of residing in the neutral system decreases over time due to the presence of the ionization rates mn (R, t ), and the drop of this probability comes from the occurrence of the first ionization.As we will now discuss, the ionization can be described by treating the quantum jumps from Ne 2 to Ne 2 + in either a stochastic or a deterministic manner, see Ref. [33] for more detail.

A. Stochastic sampling
For the stochastic sampling strategy, the basic idea is that whether the jump occurs or not at a given time step depends on the comparison between the drop in the probability, i.e., dP , at that time step and a random number (0 < 1) picked at the same time step.If dP is smaller than for a given time step, no jump occurs and the system keeps evolving in the neutral system.The comparison between dP and continues in each time step until dP is larger than at a certain time t 1 .Then the jump occurs at t 1 and the system evolves in the singly charged system after a renormalization of the projected nuclear wave packet.It is noteworthy that there may be several jump pathways.Suppose there are M states in the neutral system and N states in the singly charged system.Then there would be M × N jump pathways with their respective + .This pathway determination strategy can be easily extended to other cases.Once the jump pathway is known, the renormalized initial nuclear wave packet in Ne 2 + at the jump time t 1 is determined by the nuclear wave packet in Ne 2 at t 1 and the ionization rates at t 1 , i.e., χ n (R, t 1 ) ∝ √ mn (R, t 1 )χ m (R, t 1 ).The evolution of the nuclear wave packets in the singly charged system is similarly obtained by solving Eq. ( 3) from t = t 1 .The probability in Ne 2 + is also decreasing over time due to the second ionization from Ne 2 + to Ne 2 2+ .We evolve the nuclear wave packets in the singly charged system until a second jump takes place at t 2 , i.e., the drop of the probability at t 2 is larger than a random number generated at t 2 .There may be also several second jump pathways and the pathwaydetermination strategy is the same as that for the first jumps.The initial nuclear wave packet in the doubly charged system can be obtained in a similar way as discussed for the singly charged system.When the drop of the probability at a certain time step t 3 is larger than a random number, the third jump occurs and the system will be in a triply charged system.Similarly, we pick up a new random number to determine which jump pathway to take if there are several pathways.Once Ne 2 3+ is reached, i.e., the destination channel in this work, we achieve a successful trajectory, which is uniquely specified by the first, second, and third jump times and their respective jump pathways.By projecting the nuclear wave packet χ fi (R, t 3 ) in Ne 2 3+ just after the jump on the eigenstates χ E (R) of the potential energy curves involved in Ne 2 3+ , one obtains the nuclear KER spectrum of a single trajectory, i.e., , where E represents the eigenenergy of the eigenstate χ E (R) and is the potential energy curve for the considered state in Ne 2 3+ .The nuclear KER spectra from a huge number of trajectories reaching Ne 2 3+ with different t 1 , t 2 , and t 3 are summed together to produce the prediction of the nuclear KER spectrum following triple ionization.Access to the quantum trajectories and the ionization (jump) times allows us to unravel the dynamics.For example, we showed in Refs.[29,31] that a trajectory analysis, which is a unique possibility in the MCWP method, is capable of identifying the dominant trajectories contributing the most to a given structure in the nuclear KER spectra and in this way reveal the nuclear dynamics involved.

B. Deterministic sampling
For the deterministic sampling strategy, the jumps are assumed to take place at every time step from every jump pathways with known probabilities.Let us take the example of a given deterministic trajectory, whose first jump occurs from the mth state in Ne 2 to the nth state in Ne 2 + at t 1 , the second jump from the j th state in Ne 2 + to the kth state in Ne 2 2+ at t 2 and the third jump from the lth state in Ne 2 2+ to the qth state in Ne 2 3+ at t 3 .The probability of this trajectory is obtained by a product of the probabilities for the first, second, and third jumps and the relative probabilities for the first, second, and third jump pathways, i.e., P mnj klq (t 1 , t 2 , t 3 ) = dP (t 1 ) × P 1mn (t 1 ) × dP (t 2 ; t 1 ) × P 2jk (t 2 ; t 1 ) × dP (t 3 ; t 1 , t 2 ) × P 3lq (t 3 ; t 1 , t 2 ).Here the probabilities for the three jumps are indeed the drops of the probabilities at their respective instants of time and the relative probabilities for the pathways corresponding to the second and third jumps, i.e., P 2jk and P 3lq , are obtained similarly as the above-mentioned P 1mn .Thus, in this case, the final nuclear KER spectrum is obtained by a weighted sum of the KER spectra from the deterministic trajectories as in Ref. [31].All trajectories after deterministic sampling are successful trajectories ending in the destination channel.As a result, this sampling strategy greatly reduces the number of quantum trajectories, which finally results in a great reduction of the computational efforts.It is this sampling strategy that ensures that the computational cost for simulating multiple ionization of molecules by the MCWP method reaches a manageable level and we apply this sampling strategy in this work.To reduce the computational efforts, when applying the deterministic sampling strategy, we selectively take the trajectories whose first jumps occur at the intensity extrema of the pump-probe pulses since the first jump probabilities are the largest at the intensity extrema.The second and third quantum jumps are taken to occur at every 40 time steps.

C. Parameters
To model the experiment [23], we solve Eq. (3) for each charge state.Thus, within each charge state, the potential energy curves E el m (R) and the ionization rates mn (R, t ) for the electronic states involved are required.In our model there are no dipole transitions within a given charge state.We show the potential energy curves for the electronic states of interest in Fig. 1.Our choice of the ionization rates is as follows.The R-dependent ICD rates from the cationic 2 2 + u state to Ne 2 2+ (2p −2 ) dissociative state, i.e., ICD 1u are taken from Ref. [35] for R 3.4 Å.The ICD rates from this state at internuclear distances larger than 3.4 Å are approximated by a/R 6 with a = 16.03a.u.obtained by continuation of the numerical rates given at smaller internuclear distance.This procedure is justified since the ICD rates of the neon dimer at large internuclear distances display a 1/R 6 behavior [10,22].For the 2 2 + g state in Ne + 2 , the ICD rates are assumed identical to those from the 2 2 + u state.This is a reasonable approximation since the ICD rates from these two states are very close [13,22,36].There are no data available for the R-dependent ICD rates from the excited state in Ne 2 2+ (2s −1 2p −1 ) to Ne 2 3+ in the literature and we initially approximate them by the ICD rates from the 2 2 + u state in Ne + 2 to Ne 2 2+ .For the photoionization rates p , they are obtained from the photoionization cross sections σ p by using the relation p = σ p I /ω, where I is the intensity of the applied laser pulse with angular frequency ω.The photoionization cross section for the 2s electron in Ne 2 is approximated by the photoionization cross section of the 2s electron in Ne [37], i.e., σ 2s = 0.36 Mb at the considered wavelength.The photoionization cross section of the 2p electron in the two-site doubly ionized Ne 2 2+ is approximated by the photoionization cross section of the 2p electron in Ne(2p −1 ), i.e., σ 2p/2p = 8.08 Mb at the considered wavelength [38].Since there are no available data in the literature for the photoionization cross section of the 2p electron in the highly excited 2 2 + g and 2 2 + u states of Ne + 2 at the considered wavelength, we initially assume a relatively small number, i.e., σ 2p/2s = 0.24 Mb.As we will show later, we can make a more accurate estimate for this photoionization cross section after conducting a comparison between the experimental and theoretical results for the nuclear KER spectra.

IV. RESULTS AND DISCUSSION
In this section we present the nuclear KER spectra of Ne 2 3+ by simulating the triple ionization process of Ne 2 studied previously experimentally [23].The pump and probe pulses are identical except that the probe pulse is delayed by τ compared to the pump pulse.Two sets of XUV pulses are applied in the simulation.One set is composed of coherent pulses and the other one uses noisy pulses.The latter set of pulses represents pulses from the FEL sources in the self-amplified spontaneous emission (SASE) mode [39][40][41].These pulses exhibit statistically varying temporal and spectral shapes from shot to shot.We adopt the partial coherence method [42] to model the noisy FEL pulses.The laser parameters of the applied XUV pulses are chosen to be identical to those stated in the experimental work [23]: The central frequency is ω = 2.14 a.u., the peak intensity I 0 = 10 12 W/cm 2 , and the FWHM of the temporal intensity profile is τ p = 60 fs.
We first present the results for the coherent pulse model and consider the noisy case later.We show in Fig. 2 the KER spectra of coincident Ne 2+ and Ne + pairs as a function of delay between the pulses τ when coherent XUV pump and probe pulses irradiate Ne 2 .We take Gaussian-enveloped coherent pulses and the temporal electric field for the pump pulse is ) where F 0 is the peak field strength.The coherent pulse is show in Fig. 5(a).There is a nice agreement between our spectra [Fig.2(a)] and the experimental spectra [Fig.2(a) in Ref. [23]].The two important features in the measured spectra, i.e., the energy-decreasing branch with maximum at around 5 eV for τ 400 fs and the delay-independent high-energy branch (peaking around 8.5 eV), are reproduced at similar KER regions.The KER of the final Ne 2+ and Ne + fragments is mainly composed of three parts: the kinetic energy obtained in Ne 2 + before the second ionization E k1 , the kinetic energy obtained in Ne 2 2+ before the third ionization E k2 , and the 2/R Coulomb repulsion energy between the two fragments just after the third ionization.The potential energy curves for the 2 2 + g and 2 2 + u states in Ne 2 + (Fig. 1) are very flat around the equilibrium internuclear distance R = 6 a.u. in the potential energy curve for the ground state in Ne 2 and thus E k1 obtained through the nuclear wave packet propagating in Ne 2 + is comparably small.For the triple ionization process when the second ICD channel (ICD 2 in Fig. 1) is involved, the potential energy curve for the effective state in Ne 2 2+ (2s −1 2p −1 ) is also very flat.Thus the nuclear wave packet moves very slowly along this curve and the kinetic energies E k2 obtained for this case are also very small.As a result, the KER values of the coincident Ne 2+ and Ne + fragments from an initial nuclear wave packet in Ne 2 3+ center around 1/R 1/6 9 eV.Things are quite different when the first ICD channel (ICD 1 in Fig. 1) is involved.After ICD via this first channel, the nuclear wave packet moves relatively fast along the potential energy curve for the effective state in Ne 2 2+ (2p −2 ) (this curve is parallel to the 1/R curve at large internuclear distances) to relatively large internuclear distances.This results in relatively large kinetic energies E k2 1/R 1 − 1/R 2 and relative small Coulomb repulsion energies 2/R 2 in Ne 2 3+ , where R 1 (around 6 a.u.) and R 2 are the internuclear positions for the initial nuclear wave packets in Ne 2 2+ (2p −2 ) and Ne 2 3+ , respectively.Thus, the KER values of the Ne 2+ and Ne + fragments from an initial nuclear wave packet in Ne 2 3+ center around by 1/6 + 1/R 2 .Based on the above analysis, the signal in the high-energy branch is mainly from two possibilities: (1) trajectories where ICD from Ne + 2 (the second ionization) and the second photoionization (the third ionization) from Ne 2 2+ occurs very quickly, say, within the same pump or probe pulse, i.e., R 2 6 a.u., and (2) trajectories whose second jumps take place via photoionization from Ne 2 + to Ne 2 2+ and third jumps via the second ICD channel.It is the photoionization cross section from Ne 2 + (2s −1 ) to Ne 2 2+ (2s −1 2p −1 ) that deter-mines whether this possibility gives a small or large contribution.For the case of σ 2s/2p = 0.24 Mb, this possibility is very small and its influence cannot be seen in the calculated nuclear KER spectra.The signal in the delay-dependent branch mostly comes from trajectories where the time difference between the ICD from Ne 2 + , i.e., the ICD 1 channel, and the second photoionization from Ne 2 2+ are comparatively large, meaning that the second photonionization takes place during the probe pulse.The decreasing behavior of the KER values for this branch for an increase of delay comes from the fact that when the delay is larger, the nuclear wave packet evolves a longer time in Ne 2 2+ (2p −2 ) to larger internuclear distances R 2 , which results in smaller KER values converging to around 1/6 = 4.5 eV for the coincident Ne 2+ and Ne + fragments.
In Fig. 2(b), we plot a curve which integrates the signals in the energy-decreasing branch, i.e., for KER values between 3 and 7 eV, as a function of τ .The lifetime of the process denoted ICD 1 in Fig. 1 can be estimated by monitoring the signal in the Ne 2+ (2p −2 )-Ne + (2p −1 ) channel as a function of delay.This is possible because the transition from the Ne + (2p −1 )-Ne + (2p −1 ) channel to the Ne 2+ (2p −2 )-Ne + (2p −1 ) channel is nearly independent of internuclear distance R, and therefore of delay.The lifetime can now be extracted by fitting the signal in the Ne 2+ (2p −2 )-Ne + (2p −1 ) channel, i.e., in the energy-decreasing branch, to an exponential function a − b exp (− τ τ ICD ).This procedure results in the curves in Fig. 2(b).The fitting gives τ ICD = 102 fs, which, within the experimental uncertainty, reproduces the lifetime of 150 ± 50 fs reported experimentally [23].
There are also discrepancies between the theoretical and experimental results.We follow the authors of Ref. [23] and present a sum over spectra in the limit of large delays (τ 300 fs) in Fig. 2(c).Note that for τ 300 fs all the spectra look very similar.We cannot, however, cover delays up to 600 fs as in Ref. [23] since such large delays impose very challenging conditions on the size of our simulation box.Comparing Fig. 2(c) to Fig. 2(c) in Ref. [23], we see that the ratio of the KER signals between the low-energy peak at around 5 eV and the high-energy peak at around 8.5 eV is around 7 in our simulation, while the experimental ratio is around 2. This means that our simulation either overestimates the low-energy peak or underestimates the high-energy peak at around 8.5 eV.

A. Effects of values of photoionization cross sections and ICD rates for the ICD 2 channel
Before we go directly to the discrepancy between the theoretical and experimental ratios for the low-and highenergy peaks, we present a discussion on the delay-dependent ratios in the following.As we mentioned before, the signals in the high-energy branch in Fig. 2(a) are mainly from the trajectories where the ICD from Ne 2 + and the photoionization from Ne 2 2+ occur within the same pump pulse or within (near and within for ICD) the same probe pulse.An increase of the pump-probe delay would decrease the contribution from the second case.This is because, for a larger τ , a smaller fraction of Ne 2 + remains to undergo ICD to Ne 2 2+ when the probe pulse arrives.In the very extreme case τ → ∞, the population g and 2 2 + u states to Ne 2 2+ at τ = 400 fs when the ICD rates are identical to the rates used in Fig. 2.
in Ne 2 + has totally decayed before the arrival of the probe pulse and thus in this case the signal in the high-energy branch comes only from trajectories where the ICD and the second photoionization takes place within the pump pulse.Therefore, the ratios of the signals between the low-and high-energy peaks in the nuclear KER spectra for very large delays can also be used to reflect the ICD lifetime in Ne 2 .Generally speaking, a smaller ICD lifetime would lead to smaller ratios of the signals between the low-and high-energy peaks for very large delays.Now we turn to resolving the discrepancy of the ratios of the low-and high-energy peaks between theory and experiment.As we mentioned before, another ICD channel is open during the triple ionization of Ne 2 (ICD 2 in Fig. 1).The trajectories via this channel are specified by the first photoionization from Ne 2 to Ne 2 + , the second photoionization from Ne 2 + to Ne 2 2+ and the following ICD from Ne 2 2+ to Ne 2 3+ .It is clear from Fig. 1 that the potential energy curve for the quasibound excited state converging to Ne 2+ (2s −1 2p −1 )-Ne possesses a shallow equilibrium minimum at R ≈ 4. Thus the nuclear wave packets along this curve would move much more slowly at smaller internuclear distances than the nuclear wave packets propagating along the potential energy curve for the Ne 2 2+ state converging to Ne + (2p −1 )-Ne + (2p −1 ).As a result, the trajectories via the second ICD channel possess larger Coulomb repulsion energies at the instants of the third ionization and they would contribute to the high-energy part of the coincident Ne 2+ and Ne + signals.To make these trajectories more important, we can increase the photoionization cross sections σ 2s/2p from the excited 2 2 + g and 2 2 + u states to Ne 2 2+ .We present in Fig. 3 the nuclear KER spectra for increasing σ 2s/2p from 0.24 to 8 Mb at τ = 400 fs when the ICD rates are same as in Fig. 2. Just as expected, the ratio of the signals between the low-and high-energy peaks becomes smaller for a larger σ 2s/2p .When σ 2s/2p is about 8 Mb, the ratio is close to the measured result of 2. The discrepancy in the ratios between the theoretical and experimental low-and high-energy peaks in the nuclear KER spectra for large delays could indicate that the photoionization cross section from the 2 2 + g and 2 2 + u states is around 8 Mb at the considered photon energy.
It has been shown in Ref. [43] that a comparison between theoretical and experimental nuclear KER spectra allows an estimate for the lifetimes of the ICD processes of interest.We can clearly see from Fig. 3 that when the ratio of the signals between the low-and high-energy peaks is close to the experimental ratio of 2, the signal at around 9 eV is mainly from triple ionization of Ne 2 involving the ICD 2 channel.The contribution from the trajectories via the ICD 2 channel is about four times larger than that via the ICD 1 channel.Thus we can use the signal at around 9 eV in the nuclear KER spectra to estimate the ICD rates for the ICD 2 channel.In the following we will investigate the influence of the ICD rates for the ICD 2 channel, i.e., ICD 2 , on the nuclear KER spectra following triple ionization of Ne 2 when the photoionization cross section from the 2 2 + g and 2 2 + u states in Ne 2 + to Ne 2 2+ are taken as σ 2s/2p = 8 Mb.In Fig. 4, we show the nuclear KER spectra for the coincident Ne + and Ne 2+ fragments at τ = 400 fs when the rates for the ICD 2 channel are taken as 1, 2, 3, and 4 × ICD 1u , respectively.Also shown in in Fig. 4 is the measured spectrum extracted from Ref. [23].We can see that an increase of the ICD rates for the ICD 2 channel results in a slight decrease of the width of the peak at around 9 eV in the spectrum.When the rates for the ICD 2 channel are taken as 2 × ICD 1u or larger, the spectra agree well with the experimental result.This means that the ICD rates for the second ICD channel are expected not to be smaller than 2 × ICD 1u .In the following, the nuclear KER spectra are obtained by taking 2 × ICD 1u as the ICD rates for the ICD 2 channel.

B. Effects of chaotic SASE FEL pulses
We are now ready to investigate the interaction of Ne 2 with chaotic SASE FEL pulses.The partial coherence method [42] to model these noisy pulses is implemented as follows: From the average spectrum, ideally measured experimentally I (ω), we obtain the average electric field spectrum F (ω) = √ I (ω).Then we introduce random noise to F (ω) by adding a random spectral phase φ(ω) between −π and π , i.e., F 0 (ω) = F (ω) exp[−iφ(ω)].The corresponding temporal electric field F 0 (t ) is obtained by a Fourier transform of F 0 (ω).This field is infinitely long in the time domain due to the random spectra phase.To generate pulses with the same duration as mea- sured experimentally, a Gaussian temporal-amplitude filtering function G(t ) = exp{−4 ln 2[(t − t c )/t p ] 2 } with a duration equal to the known duration of the FEL pulse (t p = 60 fs) is multiplied by the temporal field F 0 (t ).Thus the final temporal electric field is F (t ) = F 0 (t )G(t ) for a particular realization of random phases φ(ω).The average spectrum for the FEL pulses used experimentally [23] was not available so we assume a Gaussian shape with a FWHM of ω = 0.14 a.u., i.e., I (ω) ∝ exp{−4 ln 2[(ω − ω 0 )/ ω] 2 } with ω 0 = 2.14 a.u. the center angular frequency of the FEL pulse.As an example, we show four random pump pulses obtained by the partial coherence method in Figs.5(b) to 5(e).For comparison, the Gaussian-enveloped coherent pump pulse is plotted in Fig. 5(a).
We finally present in Fig. 6 the nuclear KER spectra for Ne 2+ and Ne + pairs when the XUV pulses are noisy pulses from the SASE FEL sources.For the presented results we used 51 random pulses (one for each delay).For these results we used 2 × ICD 1u as the ICD rates for the second ICD channel and 8 Mb as the photoionization cross sections from the 2 2 + g and 2 2 + u states in Ne 2 + to Ne 2 2+ .The other parameters are those used for Fig. 2. Because of the random pulses from shot to shot, the nuclear KER spectra are not smooth with respect to the pump-probe delay.In spite of the loss of smoothness, the two main structures: the high-energy branch and the energy-decreasing branches are reproduced.Different from the coherent case in Fig. 2, we observe a ratio of the signals between the low-and high-energy peaks close to 2 in Fig. 6 smoothness makes it difficult to extract the ICD lifetime from Fig. 6(b).We can average over many nuclear KER spectra, which are individually obtained in the same way as Fig. 6, to reduce the random noise and the average spectrum can be used to extract the ICD lifetime.We expect to extract a similar ICD lifetime of around 100 fs from the averaged spectrum since the ICD process is independent of the external fields.The computational efforts are greatly increased since the number of first jumps taking place at intensity extrema are considerably increased when the light is no longer coherent.

V. CONCLUSION
In this work, we simulate the process of triple ionization, including an ICD process, of the neon dimer investigated in an XUV-pump-XUV-probe experiment by applying the MCWP approach.Characteristic structures of the nuclear KER spectra for the coincident Ne 2+ and Ne + fragments measured at different pump-probe delays are captured by our simulation.The measured lifetime for the 2s inner-valence vacancy in the neon dimer was also reproduced by our calculation.In addition, we observe a discrepancy in the theoretical and experimental ratios of the signals between the low-and high-energy peaks in the nuclear KER spectra for large delays between the pump and the probe pulses τ 400 fs.An analysis of this discrepancy revealed that another ICD channel could also be involved, i.e., ICD following double photoionization of Ne 2 .The decay rate for this ICD channel was estimated thorough a comparison between the theoretical and experimental spectra.We expect that the nuclear KER spectra following triple ionization of Ne 2 would agree better with the experimental spectra in Ref. [23] if more accurate data for the photoionization cross sections and possibly for the ICD rates were available.To summarize, the MCWP simulations in combination with the experimental results allow a critical assessment of the accuracy of the input data, e.g., the ICD rates and photoionization cross sections as well as an identification of reaction pathways.
describes the lasersystem interaction with μ the dipole operator and F (t ) the electric field of the applied laser pulse.The non-Hermitian term H d = − i 2 mn C † mn C mn is responsible for the decay from the |φ el R,m states in a given charge state to the |φ el R,n states in the charge state with one electron less.C mn = d R mn ( R, t )|φ el R,n φ el R,m | ⊗ | R R| is the quantum jump operator specifying a transition from |φ el R,m to |φ el R,n .Here | R is the position eigenket of the nuclear coordinate and mn ( R, t ) denotes the decay or ionization rates from |φ el R,m to |φ el R,n .The wave function in each charge state satisfies the time-dependent Schrödinger equation (TDSE) of the Hamiltonian H , i.e., i| ˙ (t ) = H | (t ) , where the dot on top of the ket on the left-hand side denotes the time derivative.We make use of the ansatz

FIG. 2 .
FIG. 2. (a) Calculated nuclear KER spectra for coincident Ne 2+ and Ne + fragments following triple ionization of Ne 2 as a function of delay τ between coherent XUV pump and probe pulses (see text for laser parameters).(b) The solid line shows yields of coincident Ne 2+ and Ne + pairs integrated over KER values between 3 and 7 eV in (a).The dashed line is an exponential fit of the solid line to extract the lifetime of the ICD in Ne 2 .(c) Sum of KER spectra for delays between 300 and 420 fs in (a) (see text).

σ
FIG.3.Nuclear KER spectra for different photoionization cross sections σ 2s/2p from the 2 2 + g and 2 2 + u states to Ne 2 2+ at τ = 400 fs when the ICD rates are identical to the rates used in Fig.2.

FIG. 6 . 2 2+are 8
FIG. 6.(a) Calculated nuclear KER spectra for coincident Ne 2+ and Ne + fragments following triple ionization of Ne 2 as a function of τ when the applied laser pulses are modeled by the partial coherence method (see text for laser parameters).(b) Ne 2+ and Ne + yields integrated over all KERs for the delays in (a).(c) Sum of nuclear KER spectra for delays between 300 and 420 fs.The ICD rates for the second ICD channel are taken as 2 × ICD 1u and the photoionization cross sections from the 2 2 + g and 2 2 + u states in Ne 2 + to Ne 2 2+ are 8 Mb.The other parameters are those used to obtain the results presented in Fig. 2.
. .., M and n = 1, 2, ..., N. Which pathway to follow in the jump is determined by a comparison between the P 1mn s and another random number η (0 η < 1).We here illustrate a pathway determination strategy by considering the simple case M = N = 2. (1) If 0 η < P 111 is satisfied then the jump takes place from the m = 1 state in Ne 2 to the n = 1 state in Ne 2 + .(2) If P 111 η < P 111 + P 112 , then the jump takes place from the m = 1 state in Ne 2 to the n = 2 state in Ne 2 + .(3) If P 111 + P 112 η < P 111 + P 112 + P 121 then the jump takes place from the m = 2 state in Ne 2 to the n = 1 in Ne 2 + .(4) If P 111 + P 112 + P 121 η < 1 then the jump takes place from the m = 2 state in Ne 2 to n = 2 in Ne 2