Bayesian Analysis of Femtosecond Pump-Probe Photoelectron-Photoion Coincidence Spectra with Fluctuating Laser Intensities

Abstract

This paper employs Bayesian probability theory for analyzing data generated in femtosecond pump-probe photoelectron-photoion coincidence (PEPICO) experiments. These experiments allow investigating ultrafast dynamical processes in photoexcited molecules. Bayesian probability theory is consistently applied to data analysis problems occurring in these types of experiments such as background subtraction and false coincidences. We previously demonstrated that the Bayesian formalism has many advantages, amongst which are compensation of false coincidences, no overestimation of pump-only contributions, significantly increased signal-to-noise ratio, and applicability to any experimental situation and noise statistics. Most importantly, by accounting for false coincidences, our approach allows running experiments at higher ionization rates, resulting in an appreciable reduction of data acquisition times. In addition to our previous paper, we include fluctuating laser intensities, of which the straightforward implementation highlights yet another advantage of the Bayesian formalism. Our method is thoroughly scrutinized by challenging mock data, where we find a minor impact of laser fluctuations on false coincidences, yet a noteworthy influence on background subtraction. We apply our algorithm to data obtained in experiments and discuss the impact of laser fluctuations on the data analysis.

1. Introduction

Coincidence measurements are a widely-used and powerful experimental technique in physics and chemistry. Photoelectron-photoion coincidence (PEPICO) spectroscopy utilizes not only information obtained from the detection of electrons and ions, but also the fact that they stem from the very same ionization event [1,2,3,4,5,6]. Frequently used in photoionization studies of gas phase molecules or clusters, this technique allows for conclusions about the ionization process such as the disentanglement of competing intramolecular relaxation channels [4,7,8,9] or multiple species [10], a depth of insight that cannot be achieved without assigning electrons to the ions from which they originate. Thus, the success of PEPICO is based on these recordings of pairs being unambiguous, energy-resolved for electrons, and mass-resolved for the cations. Yet, the correct pairwise assignment (true coincidence) may be affected by certain experimental conditions: If a laser pulse triggers a number of simultaneous ionization events arising from different neutral molecules, the assignment of correlated electron-ion pairs is impaired and causes so-called false coincidences [11]. The three possible events are: (1) Not exactly one electron and one ion are detected; in this case the event is rejected. (2) One electron and one ion are detected, which can originate from the same molecule (true coincidence) or (3) an electron and ion can originate from different molecules (false coincidence). Let this be illustrated using the example of exactly two ionization events. If two electrons and/or two ions are detected, the measurement would simply be discarded since no unambiguous assignments could be made. Yet, with imperfect detectors, there is a non-negligible probability of the following event: The electron from Molecule 1 is detected, the electron from Molecule 2 is not detected; Cation 1 is not detected; Cation 2 is detected. Hence, the experimentalist sees a false coincidence, where Electron 1 is wrongly assigned to Cation 2. Obviously, false coincidences only arise if both detectors are not perfect, i.e., the detection probabilities are less than unity, and are thus to some extent present in any such experiment. An easy way out is to work with low ionization rates, for the price of either a bad signal-to-noise ratio or time-consuming measurements. In principle, detector noise or ionization events not caused by the laser pulse might also lead to false coincidences, but, for the situation at hand, are sufficiently low to be neglected.

Time-resolved studies are typically carried out as pump-probe experiments [10,12] as depicted in Figure 1. Excitation by a laser pulse, commonly referred to as the pump pulse, triggers dynamical processes in the molecule after which a time-delayed second laser pulse, commonly referred to as the probe pulse, ionizes the molecule. The transient change of photo-electron and -ion signals associated with the excited states, as a function of the time-delay, provides insight into the underlying processes. Unfortunately, pump and/or probe pulses on their own can ionize the molecule as well, leading to signals that are referred to as pump-only and probe-only further on. This background signal is superimposed on the excited state signal, and in many experimental situations, the pump-only and/or the probe-only signals significantly contribute to the pump-probe signal, e.g., if multiphoton transitions are applied for pumping or probing, or if high photon energies are used for probing [13]. In order to extract the excited state transients, the pump-only and/or the probe-only signals are measured separately and usually subtracted from the pump-probe signals, obviously resulting in increased noise if the background and pump-probe spectra energetically overlap with each other. The application of the Bayesian formalism to background subtraction alone was already presented for astrophysical applications [14,15] and for photo-induced X-ray emission spectroscopy (PIXE) [16,17,18,19,20,21]. Yet, here, it has to be considered that the pump-only, the probe-only, and the pump-probe measurements have different ionization rates. Since the statistics of the coincidences depend on the ionization rates, the statistics of sole pump-only and probe-only measurements differ from those in the pump-probe measurement, and simple subtraction turns out not to be an unbiased estimator any longer [22].

Although it is feasible to distinguish between true and false coincidences by pure experimental finesse, such as in cold target recoil ion momentum spectroscopy (COLTRIMS) [23,24], it demands quite a technical and financial effort and is entirely impossible for time-of-flight detection, as used in the presented experiment. Covariance mapping, which is based on the calculation of the covariance for the photoelectron and mass spectra measured with each laser shot [25,26,27], does not guarantee that the reconstructed spectrum is positive, is restricted to Poisson processes, and leads to systematic deviations in other scenarios [26,27]. Further limitations are outlined in [25].

We recently presented a Bayesian approach to PEPICO, which treats both coincidences and background subtraction on the same footing [22,29]. In this work, we extend our theory to the prevalent experimental situation of unstable laser intensity, i.e., ionization rates fluctuating from pulse to pulse. We provide our software, including introductory examples, at https://github.com/fslab-tugraz/PEPICOBayes/. The experiment described in Section 2 is treated with the Bayesian formalism developed in Section 3. It will be tested by some challenging mock data in Section 4 and applied to real experimental data in Section 5.

2. Experiment

The analyzed experiment is of the type depicted in Figure 1 and described in detail in previous publications [7,28]. To apply our method, we choose our specimen and excitation-ionization-scheme according to Figure 2, since we expect the effects described above to be of particular importance in this scenario. Acetone molecules are excited by a three-photon transition to high-lying Rydberg states and ionized in the extraction region of a time-of-flight spectrometer, which measures both the electron kinetic energies and the ion masses [4,7]. The electron and ion flight times are then analyzed by a coincidence algorithm to produce photoelectron spectra corresponding to either an intact parent acetone ion or a fragmented acetyl ion. The concept and the necessity of time-resolved PEPICO (TR-PEPICO) in this context are further elucidated in Animations A1 and A2 in the Supplementary Material.

The addressed excited state lies energetically close to the ionization continuum, resulting in a certain probability of four-photon ionization from the ground state caused by the pump pulse alone (measurement $\alpha$ and Channel 1), leading to a background signal. In a separate pump-probe measurement (measurement $\beta$), the pump process is the same as in the pump-only case. The population is generated in the excited states, which in turn is ionized by a time-delayed probe pulse. Due to the low laser intensity of the probe pulse, ground state molecules are not ionized by the probe pulse alone, i.e., there is no probe-only background. Consequently, the measured pump-probe spectrum consists of pump-only ionization events (Channel 1) and pump-probe ionization events (Channel 2). Cations, produced in both channels, dependent on the ionization path, can either be stable and detected as parent ions or undergo fragmentation into neutral and ionic fragments. Coincidence detection of electrons and ions allows obtaining separate electron spectra for each ion, i.e., parent and fragment. The excited state of the molecule at the moment of ionization is identified by the measured electron kinetic energy in combination with the energy of the ionizing photon and knowledge of the vertical ionization energy of the excited state. In addition to the information of species and electronic state that is ionized, the related ion mass of the PEPICO spectrum provides insight into the fragmentation behavior. For example, the assignment of the photoelectron kinetic energy to an excited electronic state of the unfragmented molecule and coincidence detection of an ion fragment show that the molecule was intact at the moment of ionization and that fragmentation must have occurred afterwards. Moreover, the population in the excited state can decay to energetically lower states quite quickly, e.g., on a femtosecond timescale [4,7,28]. It is due to this decay that the Channel 2 signal can become significantly smaller than the Channel 1 background, in particular for long delay times, causing a poor signal-to-noise ratio.

3. Bayesian Data Analysis

3.1. Preliminary Considerations

We now introduce our notation and develop the Bayesian algorithm for analyzing the data generated in the experiment described in Section 1 and Section 2. We consider the following standard setup consisting of two experiments on the same target: pump-only and pump-probe, denoted by $\alpha$ and $\beta$, respectively. Each experiment consists of ${\mathcal{N}}_p$ measurements. A measurement of the $\alpha$ experiment is performed with exactly one pump pulse, while a measurement of the $\beta$ experiment comprises exactly one pump pulse and one probe pulse. During one measurement, two types of elementary coincidence events are detected, either a molecule is ionized from its ground state (referred to as Channel 1) or from its excited state (Channel 2). The latter is only possible in a pump-probe measurement ($\beta$). We assume that the number $m_j$ of ionization events in channel $j\in \{1,2\}$ is Poisson distributed with some mean ionization rate ${\lambda }_j$. For the sake of readability, we suppress the index j in the following considerations. Furthermore, we assume that in each experiment, characterized by a defined delay time between pump and probe pulse, $\lambda$ is independent of the occupation of the states, which means that we neglect population depletion effects. We additionally presume the laser intensity, and thus the ionization rate $\lambda$, to fluctuate. Mikosch and Patchkovskii [26,27] proposed to describe $\lambda$ with a Gaussian Probability Density Function (PDF). We rather resort to a $\mathsf{\Gamma }$-distribution of the latter instead, since: (1) the assignment of a Gaussian PDF for $\lambda$ is inconsistent with the fact that $\lambda \ge 0$, while the $\mathsf{\Gamma }$-distribution includes this constraint naturally; (2) it turns out to be quite convenient mathematically, while the effect of this choice on the result is deemed negligible, which will become apparent later on. We parametrize the $\mathsf{\Gamma }$-distribution, $p_{\mathsf{\Gamma }}(\lambda |\underline{\lambda },\sigma )$, with their expectation value $\underline{\lambda }$ and variance ${\sigma }^2$. All results match the findings of our previous paper [22] in the limit $\sigma \to 0$.

In one elementary event, the involved molecule can have mass $M_{\mu }$ and the emitted electron energy $E_{\nu }$. For brevity, we will refer to this particular event as $(\mu \nu )$. The ion masses and the electron energies are discretized, $\mu ,\nu \in \mathbb{N}$, due to the finite resolution of the time-of-flight spectrometer. We will also use the symbol $\rho \in \{\alpha ,\beta \}$, if we refer to the measurements/experiments $\alpha$ or $\beta$, the symbol j for Channel 1 or 2, and x for the combination of both sets, i.e., $x\in \{1,2,\alpha ,\beta \}$. Given that an elementary event happens in channel $j\in \{1,2\}$, the probability that it corresponds to $(\mu \nu )$ is denoted by:

$$\begin{array}{cc}\hfill q_{\mu \nu }^{(j)}& =P(M_{\mu },E_{\nu }|\mathrm{one}\mathrm{elementary}\mathrm{event}\mathrm{in}\mathrm{channel}j,\mathcal{I})\end{array}$$

where $\mathcal{I}$ denotes the conditional complex. The probabilities are properly normalized:

$$\begin{array}{c}\hfill \sum _{\mu \nu }{q}_{\mu \nu }^{(j)}=1\forall j\in \{1,2\}.\end{array}$$

(1)

In the pump-only measurement ($\alpha$), all molecules are in their respective ground state; therefore, only Channel 1 is allowed and:

$$\begin{array}{c}\hfill q_{\mu \nu }^{(\alpha )}=q_{\mu \nu }^{(1)}.\end{array}$$

We now introduce the spectrum:

$$\begin{array}{c}\hfill {\tilde{q}}_{\mu \nu }^{(1)}=P(E_{\nu },M_{\mu }|\mathrm{SC},q^{(1)},\pi ,\mathcal{I})\end{array}$$

with a subtle distinction of $q^{(1)}$ in the explicit propositions. $q^{(1)}$ is conditioned on $(\mu \nu )$ being a true single coincidence, whereas ${\tilde{q}}_{\mu \nu }^{(1)}$ is conditioned on $(\mu \nu )$ being a single coincidence, either true or false, and hence also on $\pi$. The latter is what is actually observed in a pump-only PEPICO measurement, while the former is what is desired. It will become apparent in the following how this additional condition distorts our statistics. We summarize all unknown parameters in the variable $\pi :=\{{\underline{\lambda }}_1,{\sigma }_1,{\underline{\lambda }}_2,{\sigma }_2,{\xi }_i,{\xi }_e\}$. ${\underline{\lambda }}_1$, ${\sigma }_1$, ${\underline{\lambda }}_2$, and ${\sigma }_2$ describe the fluctuations of ${\lambda }_1$ and ${\lambda }_2$ according to $\mathsf{\Gamma }$-distributions, and ${\xi }_i$ and ${\xi }_e$ are the detection probabilities of ions and electrons, respectively. The probability of detecting an electron with energy $E_{\nu }$ and an ion with mass $M_{\mu }$ in a single coincidence measurement is:

$$\begin{array}{cc}\hfill P(E_{\nu },M_{\mu },\mathrm{SC}|q^{(1)},\pi ,\mathcal{I})& =\underset{0}{\overset{\infty }{\int }}d{\lambda }_{1}P(E_{\nu },M_{\mu },\mathrm{SC}|q^{(1)},{\lambda }_1,\mathcal{I})p_{\mathsf{\Gamma }}({\lambda }_1|{\underline{\lambda }}_1,{\sigma }_1)\hfill \\ & ={\xi }_e{\xi }_i\left(q_{\mu \nu }^{(1)}\langle {\lambda }_1^1\rangle +{\overline{\xi }}_e{\overline{\xi }}_i{q}_{\mu .}^{(1)}{q}_{.\nu }^{(1)}\langle {\lambda }_1^2\rangle \right),\hfill \end{array}$$

(2)

with the marginals $q_{.\nu }^{(j)}={\sum }_{\mu }{q}_{\mu \nu }^{(j)}$ and $q_{\mu .}^{(j)}={\sum }_{\nu }{q}_{\mu \nu }^{(j)}$ and the detection failure probabilities ${\overline{\xi }}_e=1-{\xi }_e$ and ${\overline{\xi }}_i=1-{\xi }_i$. In the second line, we use that $P(E_{\nu },M_{\mu },\mathrm{SC}|q,{\lambda }_1,\mathcal{I})$ was already derived in Appendix 1 of [22]. The appearing integrals are of the type:

$$\begin{array}{c}\hfill \langle {\lambda }_j^n\rangle =\underset{0}{\overset{\infty }{\int }}d{\lambda }_j{\lambda }_j^n{e}^{-{\lambda }_j(1-{\overline{\xi }}_e{\overline{\xi }}_i)}{p}_{\mathsf{\Gamma }}({\lambda }_j|{\underline{\lambda }}_j,{\sigma }_j)\end{array}$$

(3)

and solved in Appendix A. Note that the PDF describing the laser fluctuations enters the whole algorithm only within integrals of the above type. Thus, the theory can easily be adapted to different descriptions of the fluctuations by merely exchanging those integrals. The spectrum including false coincidences, ${\tilde{q}}_{\mu \nu }^{(1)}$, is then given by:

$$\begin{array}{cc}\hfill {\tilde{q}}_{\mu \nu }^{(1)}& =P(E_{\nu },M_{\mu }|\mathrm{SC},q^{(1)},\pi ,\mathcal{I})=\frac{P(E_{\nu },M_{\mu },\mathrm{SC}|q^{(1)},\pi ,\mathcal{I})}{{\displaystyle \sum _{\mu \nu }}P(E_{\nu },M_{\mu },\mathrm{SC}|q^{(1)},\pi ,\mathcal{I})}=\frac{q_{\mu \nu }^{(1)}+{\kappa }_1{q}_{\mu .}^{(1)}{q}_{.\nu }^{(1)}}{1+{\kappa }_1}.\hfill \end{array}$$

(4)

${\kappa }_1$ turns out to be ${\kappa }_1={\overline{\xi }}_e\overline{{\xi }_i}\frac{\langle {\lambda }_1^2\rangle }{\langle {\lambda }_1^1\rangle }$. The false coincidences are represented by the term ${\kappa }_1{q}_{\mu .}^{(1)}{q}_{.\nu }^{(1)}$, and therefore ${\kappa }_1$ is a measure of the amount of false coincidences.

We now define ${\tilde{q}}^{(\beta )}$ analogously to ${\tilde{q}}^{(1)}$ and derive its dependence on $q^{(1)}$ and $q^{(2)}$, ${\underline{\lambda }}_1$ and ${\underline{\lambda }}_2$, and ${\sigma }_1$ and ${\sigma }_2$, abbreviated in the appendix as q, $\underline{\lambda }$, and $\sigma$. Therefore, we have to reassign the single coincidence event to the electron stemming from channel $i\in \{1,2\}$ and the ion from channel $j\in \{1,2\}$, respectively.

$$\begin{array}{cc}\hfill {\tilde{q}}_{\mu \nu }^{(\beta )}& =\frac{{\displaystyle \sum _{ij}}P(E_{\nu },M_{\mu },{\mathrm{SC}}_{ij}|q^{},\pi ,\mathcal{I})}{{\displaystyle \underset{ij}{{\displaystyle \sum _{\mu \nu }}}}P(E_{\nu },M_{\mu },{\mathrm{SC}}_{ij}|q^{},\pi ,\mathcal{I})}.\hfill \end{array}$$

Distinguishing the case where the electron and ion are measured in coincidence and stem from the same ($i=j$) or a different ($i\ne j$, or $\neg i=j$) channel, we obtain in Appendix B:

$$\begin{array}{cc}\hfill P(E_{\nu },M_{\mu },{\mathrm{SC}}_{ii}|q,\pi ,\mathcal{I})& ={\xi }_e{\xi }_i\left(q_{\mu \nu }^{(i)}\langle {\lambda }_i^1\rangle +{\overline{\xi }}_e{\overline{\xi }}_i{q}_{\mu .}^{(i)}{q}_{.\nu }^{(i)}\langle {\lambda }_i^2\rangle \right)\langle {\lambda }_{\neg i}^0\rangle \hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill P(E_{\nu },M_{\mu },{\mathrm{SC}}_{i,\neg i}|q,\pi ,\mathcal{I})& ={\xi }_e{\xi }_i{\overline{\xi }}_e{\overline{\xi }}_i{q}_{.\nu }^{(i)}\langle {\lambda }_i^1\rangle q_{\mu .}^{(\neg i)}\langle {\lambda }_{\neg i}^1\rangle .\hfill \end{array}$$

(6)

In summary, we have:

$$\begin{array}{c}\hfill {\tilde{q}}_{\mu \nu }^{(\beta )}=\frac{q_{\mu \nu }^{(2)}+{\kappa }_2\left(q_{\mu .}^{(2)}{q}_{.\nu }^{(2)}+{\alpha }_{\nu }{q}_{\mu .}^{(2)}+{\beta }_{\mu }{q}_{.\nu }^{(2)}\right)+{\gamma }_{\mu \nu }}{Z_2}\end{array}$$

(7)

with the parameters:

$$\begin{array}{cccc}\hfill {\alpha }_{\nu }& =(\mathsf{\Omega }-1)q_{.\nu }^{(1)}\hspace{1em}\hspace{1em}\hfill & \hfill {\gamma }_{\mu \nu }& =\frac{\langle {\lambda }_1^1\rangle \langle {\lambda }_2^0\rangle }{\langle {\lambda }_1^0\rangle \langle {\lambda }_2^1\rangle }(1+{\kappa }_1){\tilde{q}}_{\mu \nu }^{(1)}\hfill \\ \hfill {\beta }_{\mu }& =(\mathsf{\Omega }-1)q_{\mu .}^{(1)}\hspace{1em}\hspace{1em}\hfill & \hfill Z_2& =1+{\kappa }_2(2\mathsf{\Omega }-1)+{\gamma }_{..}\hfill \end{array}$$

(8)

where ${\gamma }_{..}={\sum }_{\mu \nu }{\gamma }_{\mu \nu }$ and:

$$\begin{array}{cccc}\hfill {\kappa }_2& ={\overline{\xi }}_e{\overline{\xi }}_i\frac{\langle {\lambda }_2^2\rangle }{\langle {\lambda }_2^1\rangle }\hspace{1em}\hspace{1em}\hfill & \hfill \mathsf{\Omega }& =1+\frac{\langle {\lambda }_1^1\rangle \langle {\lambda }_2^1\rangle }{\langle {\lambda }_1^0\rangle \langle {\lambda }_2^2\rangle }.\hfill \end{array}$$

Up to now, we have defined all variables and dependencies we need for the derivation of the desired posterior distribution presented in the following sections.

3.2. The Posterior PDF

We use Bayes’ theorem to determine the posterior probability of the parameters we want to estimate, $\pi :=\{{\underline{\lambda }}_1,{\sigma }_1,{\underline{\lambda }}_2,{\sigma }_2,{\xi }_i,{\xi }_e\}$, and the spectra $q^{(1)}=\{q_{\mu \nu }^{(1)}\}$ and $q^{(2)}=\{q_{\mu \nu }^{(2)}\}$.

$$\begin{array}{c}\hfill \underset{\mathrm{Prior}}{\underbrace{p(q^{(1)},q^{(2)},\pi |\mathcal{I})}}\underset{\mathrm{Likelihood}}{\underbrace{P(D_1,D_2|q^{(1)},q^{(2)},\pi ,\mathcal{I})}}=\underset{\mathrm{Evidence}}{\underbrace{P(D_1,D_2|\mathcal{I})}}\underset{\mathrm{Posterior}}{\underbrace{p(q^{(1)},q^{(2)},\pi |D_1,D_2\mathcal{I})}},\end{array}$$

(9)

where capital P shall denote discrete and lower case p continuous distributions. $D_1$ and $D_2$ are two datasets, and $D_1$ contains the count rates $n^{(\alpha )}=\{n_{\mu \nu }^{(\alpha )}\}$ and $n^{(\beta )}=\{n_{\mu \nu }^{(\beta )}\}$. $\{n_{\mu \nu }^{(\rho )}\}$ counts how often the pair $(E_{\nu },M_{\mu })$ was detected as a single coincidence event during the experiment $\rho$. Dataset $D_2$ contains $N_{N_e,N_i}^{(\alpha )}$ and $N_{N_e,N_i}^{(\beta )}$, which counts how many measurements lead to the detection of $N_e$ electrons and $N_i$ ions during the experiments $\alpha$ and $\beta$, respectively. In this case, it is expedient to use all detected events, not just single coincidences.

3.3. The Prior PDF

To define the prior PDF in Equation (9), we proceed with:

$$\begin{array}{c}\hfill p(q^{(1)},q^{(2)},\pi |\mathcal{I})=p(q^{(1)},q^{(2)}|\pi ,\mathcal{I})p(\pi |\mathcal{I}).\end{array}$$

For the parameters’ prior, $p(\pi |\mathcal{I})$, we use Jeffreys’ prior for scale variables for ${\underline{\lambda }}_1$ and ${\underline{\lambda }}_2$ and a flat prior for ${\xi }_i$, ${\xi }_e$, ${\sigma }_1$, and ${\sigma }_2$. The parameter ranges are ${\underline{\lambda }}_1,{\underline{\lambda }}_2,{\sigma }_1,{\sigma }_2\in {\mathbb{R}}^{+}$, and $0<{\xi }_i,{\xi }_e<1$. The spectra are normalized according to Equation (1) and positive, $q_{\mu \nu }^{(1)},q_{\mu \nu }^{(2)}>0$. Our choice for the prior of the spectra,

$$\begin{array}{c}\hfill p(q^{(1)},q^{(2)}|\pi ,\mathcal{I})=p(q^{(2)}|q^{(\alpha )},\pi ,\mathcal{I})p(q^{(1)}|\pi ,\mathcal{I}),\end{array}$$

(10)

are Dirichlet priors [30] in the spectra ${\tilde{q}}_{\mu \nu }^{(\alpha )}$ and ${\tilde{q}}_{\mu \nu }^{(\beta )}$,

$$\begin{array}{c}\hfill p({\tilde{q}}^{(\beta )}|q^{(\alpha )},\pi ,\mathcal{I})=\frac{1}{B(\{c_{\mu \nu }^{(\beta )}\})}\prod _{\mu \nu }{\tilde{q}}_{\mu \nu }^{(\beta )}{}^{c_{\mu \nu }^{(\beta )}-1}\delta (\tilde{S}-1),\end{array}$$

with $\tilde{S}={\sum }_{\mu \nu }{\tilde{q}}_{\mu \nu }^{(\rho )}$ and the normalization $B(\{c_{\mu \nu }^{(\rho )}\})$ being the multivariate beta function. We can always choose the prior to be uninformative (flat) by setting all $c_{\mu \nu }^{(\rho )}=1$. The choice of these Dirichlet priors is without loss of generality and justified by the likelihood being multinomial in these variables, as we will see later. Previously, with Equation (7), we found ${\tilde{q}}_{\mu \nu }^{(\beta )}$ to be a function of $q^{(\alpha )}$ and $q^{(2)}$, i.e., ${\tilde{q}}_{\mu \nu }^{(\beta )}={\tilde{Q}}^{(\beta )}(q^{(\alpha )},q^{(2)})$. With the usual transformation rules, the first term in Equation (10) becomes:

$$\begin{array}{c}\hfill p(q^{(2)}|q^{(\alpha )},\pi ,\mathcal{I})=p({\tilde{q}}^{(\beta )}|q^{(\alpha )},\pi ,\mathcal{I})\times \left|\frac{d{\tilde{Q}}^{(\beta )}(q^{(\alpha )},q^{(2)})}{d{q}^{(2)}}\right|.\end{array}$$

(11)

The appearing transformation of the Dirac distribution,

$$\begin{array}{c}\hfill \delta (\tilde{S}-1)=\frac{Z_2}{1+2{\kappa }_2\mathsf{\Omega }}\delta (S-1),\end{array}$$

(12)

with $S={\sum }_{\mu \nu }{q}_{\mu \nu }^{(2)}$, is derived in Appendix C, and the derivation of the Jacobian determinant,

$$\begin{array}{c}\hfill \left|\frac{d{\tilde{Q}}^{(\beta )}(q^{(\alpha )},q^{(2)})}{d{q}^{(2)}}\right|=\frac{{(1+{\kappa }_2\mathsf{\Omega })}^{{\mathcal{N}}_{\mu }+{\mathcal{N}}_{\nu }-2}(1+2{\kappa }_2\mathsf{\Omega })}{Z_2^{{\mathcal{N}}_{\mu }{\mathcal{N}}_{\nu }}},\end{array}$$

(13)

is shown in Appendix D. ${\mathcal{N}}_{\mu }$ and ${\mathcal{N}}_{\nu }$ denote the total number of bins for ions $(\mu )$ and electrons $(\nu )$, respectively. Putting everything together, Equation (11) becomes:

$$\begin{array}{c}\hfill p(q^{(2)}|q^{(\alpha )},\pi ,\mathcal{I})=\frac{{(1+{\kappa }_2\mathsf{\Omega })}^{{\mathcal{N}}_{\mu }+{\mathcal{N}}_{\nu }-2}\delta \left({\sum }_{\mu \nu }{q}_{\mu \nu }^{(2)}-1\right)}{Z_2^{{\mathcal{N}}_{\mu }{\mathcal{N}}_{\nu }-1}B\left(\left\{c_{\mu \nu }^{(\beta )}\right\}\right)}\times \prod _{\mu \nu }{\left({\tilde{Q}}_{\mu \nu }^{(\beta )}\left(q^{(\alpha )},q^{(2)}\right)\right)}^{c_{\mu \nu }^{(\beta )}-1}.\end{array}$$

The second term in Equation (10), $p(q^{(1)}|\pi ,\mathcal{I})$, is obtained by setting $\mathsf{\Omega }=1$ and replacing experiment $\beta$ by $\alpha$ and Channel 2 by 1; hence, the prior is fully determined.

3.4. The Likelihood

In this paper, we only use the single coincidence events $D_1$ for estimating the spectra $q^{(1)}$ and $q^{(2)}$ given the parameters $\pi$. This is justified by the fact that especially these events include relevant information about the spectrum, because the case of detecting more than one electron-ion pair does not allow one to link an electron to the ion it originates from, and the Bayesian approach would be different. The dataset $D_2$ will be used to determine the unknown parameters $\pi$. Therefore, the likelihood splits into:

$$\begin{array}{cc}\hfill P(D_1,D_2|q^{(1)},q^{(2)},\pi ,\mathcal{I})& =P(D_1|q^{(1)},q^{(2)},\pi ,\mathcal{I})P(D_2|\pi ,\mathcal{I}).\hfill \end{array}$$

(14)

The first term separates into:

$$\begin{array}{cc}\hfill P(D_1|q^{(1)},q^{(2)},\pi ,\mathcal{I})& =P(n^{(\alpha )}|q^{(1)},\pi ,\mathcal{I})P(n^{(\beta )}|q^{(1)},q^{(2)},\pi ,\mathcal{I}).\hfill \end{array}$$

Marginalizing over ${\tilde{q}}_{\mu \nu }^{(\alpha )}$ and ${\tilde{q}}_{\mu \nu }^{(\beta )}$ introduces the multinomial distribution:

$$\begin{array}{c}\hfill P(n^{(\beta )}|{\tilde{q}}^{(\beta )},\pi ,\mathcal{I})=\mathcal{M}(\{n^{(\beta )}\}|\{{\tilde{q}}^{(\beta )}\},\pi ,\mathcal{I})\end{array}$$

according to Appendix 2 of [22]. We now turn to the second term in Equation (14),

$$\begin{array}{cc}\hfill P(D_2|\pi ,\mathcal{I})& =P(\{N_{N_e,N_i}^{(\alpha )}\}|\pi ,\mathcal{I})P(\{N_{N_e,N_i}^{(\beta )}\}|\pi ,\mathcal{I}).\hfill \end{array}$$

With ${\mathcal{N}}_P$ being the number of laser pulses; the terms are multinomial distributions, as well,

$$\begin{array}{c}\hfill P(\{N_{N_e,N_i}^{(\rho )}\}|\pi ,\mathcal{I})=\mathcal{M}(\{N_{N_e,N_i}^{(\rho )}\}|\{{\tilde{P}}_{N_e,N_i}^{(\rho )}\},{\mathcal{N}}_P^{(\rho )}).\end{array}$$

Since the PDF is the same for both experiments, $\alpha$ and $\beta$, we suppress the index $\rho$ in the following considerations. The probability of detecting $N_e$ electrons and $N_i$ ions is:

$$\begin{array}{c}\hfill {\tilde{P}}_{N_e,N_i}:=P(N_e,N_i|\underline{\lambda },\sigma ,{\xi }_e,{\xi }_i,\mathcal{I})=\underset{0}{\overset{\infty }{\int }}d\lambda \underset{=:P_{N_e,N_i}}{\underbrace{P(N_e,N_i|\lambda ,{\xi }_e,{\xi }_i,\mathcal{I})}}{p}_{\mathsf{\Gamma }}(\lambda |\underline{\lambda },\sigma ).\end{array}$$

(15)

$P_{N_e,N_i}$ is already derived in Appendix 4 in [22] and has the general form:

$$\begin{array}{c}\hfill P_{N_e,N_i}=\sum _{n=\max (N_e,N_i)}^{N_e+N_i}{a}_n({\xi }_1,{\xi }_i){\lambda }^n{e}^{-\lambda (1-{\overline{\xi }}_e{\overline{\xi }}_i)}.\end{array}$$

Inserting in Equation (15) produces:

$$\begin{array}{cc}\hfill {\tilde{P}}_{N_e,N_i}& =\sum _{n=\max (N_e,N_i)}^{N_e+N_i}{a}_n({\xi }_e,{\xi }_i)\langle {\lambda }^n\rangle \hfill \end{array}$$

A list of ${\tilde{P}}_{N_e,N_i}$ for $max(N_e,N_i)\le 3$ is given in Appendix E.

3.5. Remarks on the Posterior Sampling

In the previous sections, we fully determined the posterior $p(q^{(1)},q^{(2)},\pi |D_1,D_2,\mathcal{I})$. Ultimately, we are interested in the spectrum $q^{(2)}$, or rather say its probability $p(q^{(2)}|D_1,D_2,\mathcal{I})$, which can easily be achieved by integrating out $q^{(1)}$ and $\pi$ in the posterior. For performing this integration and, more generally, computing expectation values, variances, and covariances of the posterior, we resort to numerical techniques, posterior sampling specifically. A suitable technique for sampling from this PDF is Markov Chain Monte Carlo (MCMC), which is based on a Markov chain that has the desired distribution as its equilibrium distribution. The technique is standard in Bayesian probability theory (see [30], especially Chapter 30 and the references therein). To generate the Markov chain, we used the Metropolis–Hastings algorithm with local updates in $q^{(1)}$, $q^{(2)}$, and $\pi$. We chose the step size for the update of each parameter separately to achieve an acceptance probability of ca. 50%. We discarded the first 20% of the Markov chain as means of thermalization and tested the convergence to the equilibrium distribution with several initial states. Correlations were checked with binning.

4. Mock Data Analysis

In this section, we demonstrate the performance of our algorithm including $\lambda$-fluctuations. There are two disturbing influences in the reconstruction of the spectra: false coincidences and pump-only background. We study these influences separately. First, we investigate the influence of the $\lambda$-fluctuations on the false coincidences. In the second part, we scrutinize our algorithm in the presence of a challenging pump-only background and $\lambda$-fluctuations.

4.1. False Coincidences

In the simulation, we have two different ion masses $\mu \in \{\mathrm{p},\mathrm{f}\}$, called parent (p) and fragment (f). We use the test spectra related to the ones reported in [22,26]. The electronic spectrum of the parent has a step-like form, and the electronic spectrum of the fragment consists of Gaussian peaks; see Figure 3. These challenging test spectra exhibit a strongly-varying parent/fragment-ratio as a function of the electron energy. Since we are firstly interested in the effect of $\lambda$-fluctuations on false coincidences, we suppress the background in our simulations (${\underline{\lambda }}_1\to 0$ and ${\sigma }_1\to 0$). Before studying the simulation results in detail, we investigate the influence of $\lambda$-fluctuations on false coincidences by analyzing Equation (7). Since we have only signal contributions ($q^{(2)}$), Equation (7) reduces to Equation (4) with the superscripts $(2)$ and 2 instead of $(1)$ and 1, respectively. The equation describes the connection between the measured spectrum ${\tilde{q}}^{(2)}$, including false coincidences, and the true spectrum $q^{(2)}$. ${\kappa }_2$ gives the weight between false and true coincidences,

$$\begin{array}{c}\hfill {\kappa }_2={\overline{\xi }}_e\overline{{\xi }_i}\frac{\langle {\lambda }_2^2\rangle }{\langle {\lambda }_2^1\rangle }={\overline{\xi }}_e\overline{{\xi }_i}{\underline{\lambda }}_2\left\{\frac{1+\frac{{\sigma }_2^2}{{\underline{\lambda }}_2^2}}{1+(1-{\overline{\xi }}_e\overline{{\xi }_i})\frac{{\sigma }_2^2}{{\underline{\lambda }}_2}}\right\}.\end{array}$$

(16)

This equation shows that the influence of $\lambda$-fluctuations on false coincidences is negligible in two cases: first, if the relative fluctuations ${\sigma }_2^2/{\underline{\lambda }}_2^2$ are small; and second, if the relation $(1-{\overline{\xi }}_e\overline{{\xi }_i}){\underline{\lambda }}_2\approx 1$ holds, the influence is negligible regardless of high relative fluctuations.

In our mock data analysis, we choose the parameters ${\xi }_i={\xi }_e=0.5$ and ${\mathcal{N}}_p={10}^7$, ${\sigma }_2=0.5$, and ${\underline{\lambda }}_2=\{1.5,0.5\}$. Now, we turn to the parameter estimates summarized in

Table 1. Estimated parameters ${\underline{\lambda }}_2$, ${\sigma }_2$, ${\xi }_i$, and ${\xi }_e$. In the lines showing the results of the algorithm presented in [22], ${\lambda }_2$ is shown instead of ${\underline{\lambda }}_2$. Each value denotes the mean and standard deviation of the parameter’s distribution.

	λ${}_2$	${\mathit{\sigma }}_2$	${\mathit{\xi }}_{\mathit{i}}$	${\mathit{\xi }}_{\mathit{e}}$
Parameters (Figure 3a,b)	$1.5$	$0.5$	$0.5$	$0.5$
Algorithm in [22]	$1.4177\pm 0.0008$	-	$0.5235\pm 0.0003$	$0.5235\pm 0.0003$
Algorithm with $\lambda$-fluctuations	$1.499\pm 0.001$	$0.501\pm 0.002$	$0.4998\pm 0.0003$	$0.5000\pm 0.0003$
Parameters (Figure 3c,d)	$0.5$	$0.5$	$0.5$	$0.5$
Algorithm in [22]	$0.4365\pm 0.0003$	-	$0.5618\pm 0.0004$	$0.5621\pm 0.0004$
Algorithm with $\lambda$-fluctuations	$0.4992\pm 0.0005$	$0.499\pm 0.001$	$0.5000\pm 0.0004$	$0.5004\pm 0.0004$

. The algorithm reported in [22], which does not include $\lambda$-fluctuations, provides estimates for ${\lambda }_2$, ${\xi }_i$, and ${\xi }_e$ with about $5\%$ (for ${\underline{\lambda }}_2=1.5$) and $13\%$ (for ${\underline{\lambda }}_2=0.5$) relative deviation from the real values. The differences are caused by the assumption of having a constant ${\lambda }_2$ without any statistical noise, which is simply wrong, since we generated our mock data by drawing ${\lambda }_2$ from the $\mathsf{\Gamma }$-distribution with the parameters ${\underline{\lambda }}_2$ and ${\sigma }_2$. As expected, ${\lambda }_2$ estimated by the algorithm ignoring $\lambda$-fluctuations lies within ${\underline{\lambda }}_2\pm {\sigma }_2$. The algorithm reported in this paper includes the $\lambda$-fluctuations and is therefore able to give reliable estimates for all parameters.

We will now show that even if the parameter estimation of the algorithm not including $\lambda$-fluctuations is slightly off, there is only a small difference in the reconstructed spectra. While at ${\underline{\lambda }}_2=1.5$, $\lambda$-fluctuations do not affect the reconstruction (Figure 3a,b), at ${\underline{\lambda }}_2=0.5$, the algorithm including $\lambda$-fluctuations (green line) is slightly closer to the true spectrum than the algorithm not including $\lambda$-fluctuations (blue line); see Figure 3c,d. In this parameter regime, neglecting $\lambda$-fluctuations leads to ${\kappa }_2=0.125$, while Equation (16) produces ${\kappa }_2=0.182$. Therefore, the algorithm not including $\lambda$-fluctuations underestimates the amount of false coincidences.

Altogether, the influence of $\lambda$-fluctuations on false coincidences in the reconstructed spectra is small in the experimentally-relevant parameter regime, which is in accordance with Mikosch et al. [26]. Still, including $\lambda$-fluctuations is important in the parameter estimation.

4.2. Background Subtraction

To analyze the influence of $\lambda$-fluctuations on the background subtraction, we use the same spectra as in Section 4.1. Now, the step-like spectrum is chosen as the signal and the spectrum consisting of Gaussian peaks as the background. We restrict ourselves to one ion mass to exclude the influence of false coincidences stemming from two ion masses. If we neglect false coincidences, the prefactor $\frac{\langle {\lambda }_1^1\rangle \langle {\lambda }_2^0\rangle }{\langle {\lambda }_1^0\rangle \langle {\lambda }_2^1\rangle }$ in ${\gamma }_{\mu \nu }$ defined in Equation (8) is:

$$\begin{array}{c}\hfill \frac{\langle {\lambda }_1^1\rangle \langle {\lambda }_2^0\rangle }{\langle {\lambda }_1^0\rangle \langle {\lambda }_2^1\rangle }=\frac{{\underline{\lambda }}_1}{{\underline{\lambda }}_2}\left\{\frac{1+(1-{\overline{\xi }}_e\overline{{\xi }_i})\frac{{\sigma }_2^2}{{\underline{\lambda }}_2}}{1+(1-{\overline{\xi }}_e\overline{{\xi }_i})\frac{{\sigma }_1^2}{{\underline{\lambda }}_1}}\right\}.\end{array}$$

This factor represents the ratio between the background (1) and signal (2). The equation shows that the corrections cancel if ${\underline{\lambda }}_1={\underline{\lambda }}_2$ and ${\sigma }_1={\sigma }_2$; see Figure 4a,b. If ${\sigma }_1^2/{\underline{\lambda }}_1<{\sigma }_2^2/{\underline{\lambda }}_2$ (Figure 4c), the algorithm without $\lambda$-fluctuations (blue line) underestimates the background, leading to peaks in every step. In the opposite case, ${\sigma }_1^2/{\underline{\lambda }}_1>{\sigma }_2^2/{\underline{\lambda }}_2$ (Figure 4d), the algorithm without $\lambda$-fluctuations (blue line) overestimates the background, leading to notches in every step. The algorithm including $\lambda$-fluctuations (green line) is able to reconstruct the signal correctly in all scenarios.

Similar to the mock data analysis in the previous section, not including $\lambda$-fluctuations leads to deviations in the parameter estimation in all of the simulated test cases; see

Table 2. Estimated parameters ${\underline{\lambda }}_1$, ${\underline{\lambda }}_2$, ${\sigma }_1$, ${\sigma }_2$, ${\xi }_i$, and ${\xi }_e$. The parameter regimes denoted by the identifications (a–d) are according to Figure 4. For each parameter set, the first line denotes the true value, while Line 2 (3) contains the parameter estimation performed with the algorithm without (with) $\lambda$-fluctuations, respectively. In the lines showing the results of the algorithm ignoring $\lambda$-fluctuations, ${\lambda }_j$ is shown instead of ${\underline{\lambda }}_j$. Each value denotes the mean and standard deviation of the parameter’s distribution.

	λ${}_1$	λ${}_2$	${\mathit{\sigma }}_1$	${\mathit{\sigma }}_2$	${\mathit{\xi }}_{\mathit{i}}$	${\mathit{\xi }}_{\mathit{e}}$
(a)	$0.5$	$0.5$	$0.1$	$0.1$	$0.5$	$0.5$
	$0.4970\pm 0.0003$	$0.4972\pm 0.0005$	-	-	$0.5032\pm 0.0002$	$0.5029\pm 0.0002$
	$0.5000\pm 0.0003$	$0.5007\pm 0.0005$	$0.098\pm 0.002$	$0.101\pm 0.005$	$0.5003\pm 0.0003$	$0.5000\pm 0.0003$
(b)	$0.5$	$0.5$	$0.5$	$0.5$	$0.5$	$0.5$
	$0.4367\pm 0.0003$	$0.4317\pm 0.0004$	-	-	$0.5618\pm 0.0002$	$0.5620\pm 0.0002$
	$0.5002\pm 0.0004$	$0.5009\pm 0.0006$	$0.5026\pm 0.0008$	$0.496\pm 0.002$	$0.4996\pm 0.0003$	$0.4998\pm 0.0003$
(c)	$0.5$	$0.5$	$0.1$	$0.5$	$0.5$	$0.5$
	$0.4840\pm 0.0003$	$0.4608\pm 0.0005$	-	-	$0.5201\pm 0.0002$	$0.5201\pm 0.0002$
	$0.4991\pm 0.0003$	$0.5013\pm 0.0006$	$0.098\pm 0.002$	$0.500\pm 0.002$	$0.5001\pm 0.0002$	$0.5001\pm 0.0003$
(d)	$0.5$	$0.5$	$0.5$	$0.1$	$0.5$	$0.5$
	$0.4452\pm 0.0003$	$0.4675\pm 0.0004$	-	-	$0.5441\pm 0.0002$	$0.5440\pm 0.0002$
	$0.4998\pm 0.0004$	$0.5005\pm 0.0006$	$0.5002\pm 0.0009$	$0.102\pm 0.004$	$0.4998\pm 0.0003$	$0.4997\pm 0.0003$

. As expected, ${\lambda }_j$ estimated by the algorithm ignoring $\lambda$-fluctuations lies within ${\underline{\lambda }}_j\pm {\sigma }_j$. Obviously, including $\lambda$-fluctuations in the parameter estimation becomes important especially for high ${\sigma }_1$ or ${\sigma }_2$. In summary, the influence of $\lambda$-fluctuations on the background subtraction is important if the magnitude of the relative $\lambda$-fluctuations of one channel is large and does not compare to the other.

5. Application to Experimental Data

With our setup for time-resolved PEPICO studies on gas phase molecules (see Figure 2), we were not able to produce stable $\lambda$-fluctuations with a sufficiently large $\sigma$ to demonstrate a systematic influence on the reconstructed spectra. The output of our commercial laser system (Coherent Vitara oscillator and Legend Elite Duo amplifier) turned out to be too stable for this purpose. We tried to simulate a less stable setup by operating the Optical Parametric Amplifier (OPA) at intensities below the saturation threshold, which induces $\lambda$-fluctuations. In this setting, however, the pulse energy also undergoes temporal drifts on time scales of the measurement, which are different from the statistical fluctuations and not accounted for in our analysis as a setup would not be operated in this regime.

Nevertheless, by using the algorithm including the $\lambda$-fluctuations, we find that the fluctuations were in the order of $\sigma =\mathcal{O}(0.001)$ in the experimentally-relevant parameter regime. The parameters estimated by the different algorithms are depicted in

Table 3. Estimated parameters ${\underline{\lambda }}_1$, ${\underline{\lambda }}_2$, ${\sigma }_1$, ${\sigma }_2$, ${\xi }_i$, and ${\xi }_e$. Line 1 (2) contains the parameter estimations performed with the algorithm without (with) $\lambda$-fluctuations, respectively. In the line showing the results of the algorithm ignoring $\lambda$-fluctuations, ${\lambda }_j$ is shown instead of ${\underline{\lambda }}_j$. Each value denotes the mean and standard deviation of the parameter’s distribution.

λ${}_1$	λ${}_2$	${\mathit{\sigma }}_1$	${\mathit{\sigma }}_2$	${\mathit{\xi }}_{\mathit{i}}$	${\mathit{\xi }}_{\mathit{e}}$
$0.3328\pm 0.0009$	$0.132\pm 0.001$	-	-	$0.3247\pm 0.0008$	$0.542\pm 0.001$
$0.3328\pm 0.0009$	$0.1335\pm 0.0003$	$0.0004\pm 0.0010$	$0.0004\pm 0.0012$	$0.3245\pm 0.0008$	$0.541\pm 0.001$

. Finally, we note that less stable systems are expected to suffer from significant $\lambda$-fluctuations, which have to be accounted for in the reconstruction process of the spectrum. This will in particular be the case if multiple nonlinear optical processes like in the OPA are used in the pump path of a pump-probe setup. The fluctuations get even more prominent if multiphoton processes are used for the excitation process, because these processes depend nonlinearly on the pulse power.

6. Conclusions

We used Bayesian probability theory to analyze data obtained from pump-probe photoionization experiments with photoelectron-photoion coincidence detection. We extended the algorithm developed previously in [22] by including fluctuations of the laser intensity as a random variable $\lambda$. Based on challenging mock data, we have demonstrated the reliability of the developed algorithm. In accordance with Mikosch et al. [26], the influence of $\lambda$-fluctuations on false coincidences was small in a certain parameter regime. We derived the condition $(1-{\overline{\xi }}_e\overline{{\xi }_i})\underline{\lambda }\approx 1$ for negligible influence of $\lambda$-fluctuations on false coincidences, even at high relative fluctuations. ${\overline{\xi }}_e$ and ${\overline{\xi }}_i$ denote the complement of the detection probabilities of electrons (${\xi }_e$) and ions (${\xi }_i$), respectively. In the case of the pump-only background (Channel 1) and the signal (Channel 2) contained in the pump-probe measurements, neglecting $\lambda$-fluctuations underestimates the background signal in the case ${\sigma }_1^2/{\underline{\lambda }}_1<{\sigma }_2^2/{\underline{\lambda }}_2$. The other way around, the background signal is overestimated. In both scenarios, including $\lambda$-fluctuations is important for estimating the parameters ${\underline{\lambda }}_1$, ${\underline{\lambda }}_2$, ${\sigma }_1$, ${\sigma }_2$, ${\xi }_i$, and ${\xi }_e$ correctly. In our application to the experimental data, we find that the relative laser fluctuations in the experimental setup are fortunately too small to see effects on the reconstructed spectra. With the developed algorithm, we were able to determine the experimental $\lambda$-fluctuations to be in the order of $\sigma =\mathcal{O}(0.001)$. Compared to conventional subtraction of the pump-only spectrum from the pump-probe spectrum, the Bayesian approach provides several important advantages: (i) It results in a significant increase of the signal-to-noise ratio. (ii) It does not overestimate the pump-only contribution and never leads to negative spectra because the relative weight of the pump-only contribution is self-consistently determined. (iii) Spectral signatures based on false coincidences are eliminated, allowing for higher signal rates. (iv) It includes consistently all prior knowledge, such as positivity, and (v) a confidence interval is obtained for the estimated spectrum. (vi) It is applicable to any experimental situation and noise statistics, as demonstrated in this paper for the case of $\lambda$-fluctuations.