Laser–particle beam collisions have proven an effective tool in the investigation of strong-field quantum electrodynamics (SFQED). Early experiments utilized collisions between a linear accelerator (linac)-accelerated electron beam and a high-power laser^[1^] to observe non-linear Compton scattering (NLCS)^[2^] and the non-linear Breit–Wheeler (BW) process^[3^]. More recently, so-called all-optical experiments (in which one laser pulse drives a wakefield accelerator^[4^], producing a relativistic electron beam that collides with a second laser pulse) aimed to probe radiation reaction^[5^,6^], the recoil experienced by a charge accelerated in an external field upon emitting a photon. A number of initiatives^[7^,8^] aim to perform high-precision studies of NLCS and non-linear BW pair creation.

The advent of laser facilities capable of attaining ${10}^{21}-{10}^{23}$ W cm⁻² intensities^[9^–12^] will enable the exploration of SFQED in regimes where quantum effects are predicted to be substantial. Efforts are underway to measure vacuum birefringence in collisions between a high-power laser and a brilliant X-ray source^[13^,14^], while future campaigns propose to use high-power lasers to observe exotic phenomena such as photon–photon splitting^[15^] and the self-focusing and self-compression of light in vacuum^[16^].

Many SFQED experiments propose to use collisions between high-power lasers and particles beams, or between high-power lasers, and may thus experience difficulties in interpreting data due to shot-to-shot variation in collision parameters. A recent study^[17^] of the effect of varying collision parameters on model differentiation for radiation reaction studies using electron beam–laser collisions found $20\%\hbox{--} 30\%$ changes in the mean and $5\%\hbox{--} 10\%$ changes in the width of the final electron spectrum for typical variations in longitudinal alignment between the electron beam and laser pulse, while variations in the electron beam duration and chirp induced changes of up to $20\%$ in the final mean electron energy and width. This motivates the development of novel data analysis tools that account for uncertainties due to unknown or measured particle beam, laser and collision parameters self-consistently.

To address this challenge, we have developed a Bayesian inference framework that facilitates parameter inference and model comparison for all-optical NLCS experiments aiming to probe radiation reaction. This framework infers values of collision parameters that directly affect experimental observables but were not measured on-shot, and incorporates knowledge of collision parameters from prior measurements or simulations (in the form of prior distributions). This procedure combines multiple diagnostics into a single, self-consistent analysis and enables a quantitative comparison of three radiation reaction models; the classical, quantum-continuous and quantum-stochastic models outlined in Section 2. While this framework applies to all-optical radiation reaction experiments, some of the techniques used in this work have wider relevance for beam–beam or laser–beam collider experiments.

We identify a number of challenges associated with the implementation of this analysis and propose strategies to address them. Increasing the number of inference (or free) parameters rapidly increases the computational cost of the inference procedure beyond the point where the inference is tractable. For example, in Section 5.3, convergence is achieved for the slowest converging parameter after 8000 steps for the one-dimensional test cases and 31,000 steps for the three-dimensional test cases. On average, the one-dimensional inference procedures conducted for mono-energetic electron spectra required 5 CPUs, 60 GB per CPU and 8 hours of runtime, while the three-dimensional inference procedures typically required 40 CPUs, 60 GB per CPU and 480 hours of runtime.

In addition, an excessive number of free parameters may result in over-fitting. Therefore, a number of collision parameters are assigned fixed values. We assess the impact of fluctuating collision parameters on our experimental observables given their expected shot-to-shot variation. We then demonstrate that degeneracies between free and fixed parameters allow free parameters to replicate the effect on the experimental observables of a fixed parameter having an experimental value that differs from its value in the model. These two considerations inform the selection of free and fixed parameters.

When using this approach, inference parameters should be treated as effective parameters that replicate the collision conditions, rather than physical parameters that accurately represent electron beam and laser properties, or their spatio-temporal alignment. If the experimental value of a fixed parameter deviates from its value in the model by more than a given amount, we find accurate model differentiation is no longer possible. We identify the threshold at which this occurs for transverse misalignments between the electron beam and colliding laser and propose to mitigate this issue by applying the Bayesian analysis to shots with the highest photon yields, for which transverse misalignments are likely to be small.

Finally, we assess the accuracy of model selection and parameter inference using the Bayesian framework via a variety of test cases with electron, laser and collision parameters representative of a recent experiment. We find that, given the experimental uncertainties and the broad, uniform priors we opted for, single-shot model differentiation is infeasible. However, model selection may be achieved by combining model evidences over multiple shots.

Two parameters, namely the electron quantum parameter, $\eta ={E}_\mathrm{RF}/{E}_\mathrm{s}$ , and the dimensionless intensity parameter, ${a}_0=\frac{\mid {E}_\mathrm{L}\mid e}{\omega_\mathrm{L}{m}_\mathrm{e}c}$ , govern the quantum and non-linear character of electron–photon interactions, respectively^[18^]. The electron charge and mass are denoted as $e$ and ${m}_\mathrm{e}$ , respectively, ${\omega}_\mathrm{L}$ and ${E}_\mathrm{L}$ are the frequency and strength of the external electric field in the lab frame, respectively, ${E}_\mathrm{RF}$ is the electric field strength in the electron rest frame and ${E}_\mathrm{s}=1.3\times {10}^{18}\;\mathrm{V}\;{\mathrm{m}}^{-1}$ is the Schwinger field. A classical theory of radiation reaction is expected to be valid when $\alpha {a}_0\eta \simeq 1$ and $\eta \ll 1$ ^[19^,20^], where $\alpha$ is the fine structure constant, while quantum effects become dominant for $\alpha {a}_0\simeq 1$ and $\eta \gtrsim 1$ ^[19^,20^]. Multi-photon and relativistic effects manifest for ${a}_0>1$ .

Classical radiation reaction is expected to be well-described by the Landau–Liftshitz model^[21^], which treats radiation emission as a continuous process and does not impose an upper bound on the frequency of radiation emitted by an electron. For this reason, classical radiation reaction over-predicts electron energy loss compared to quantum models^[22^].

For ${a}_0>1$ , quantum models of electron–photon scattering become non-perturbative. An alternative approach is employed (the Furry picture) in which electrons are ‘dressed’ by the background field^[23^]. In order to treat arbitrary electromagnetic fields, photon emission is assumed to occur over timescales much smaller than the timeframe of electromagnetic field variation, allowing the electric and magnetic fields to be treated as locally constant and orthogonal (locally constant field approximation, or LCFA)^[24^]. Between photon emissions, electrons propagate classically^[25^–27^]. The quantum model of radiation reaction, here termed the quantum-stochastic model, prohibits the emission of photons with energies exceeding the electron energy and treats photon emission as a stochastic process, which gives rise to spectral broadening^[22^,28^,29^].

We also consider a quantum-continuous model, which incorporates first-order quantum effects in a classical framework. This model treats radiation emission as continuous but applies a correction factor, the Gaunt factor^[30^], to the radiation reaction force term to recover the same rate of change of average electron momentum predicted by the quantum-continuous model^[22^,28^]. Both the classical and quantum-continuous models predict spectral narrowing^[22^,28^].

Throughout this paper, the subscripts $\mathrm{qs}$ , $\mathrm{qc}$ and $\mathrm{cl}$ refer to the quantum-stochastic, quantum-continuous and classical models, respectively.

Bayesian inference is a statistical technique that allows the unknown variables, $\rho$ , which parameterize a given model, ${M}^x$ , to be inferred. Central to this technique is the approach of using a forward model (which predicts the experimental observables for a set of inputs) to obtain the probability that the model is accurate given the observed data, $D$ , that is, $P\left({M}^x|D\right)$ , the posterior probability. Crucially, the forward models, and hence their posterior probability distributions, are functions of the inference parameters. The posterior probability is calculated using Bayes’ theorem:

$$\begin{array}{r}P\left(M^x|D\right)=\frac{P\left(D|M^x\right)P\left(M^x\right)}{P(D)},\end{array}$$  (1)

$$\begin{array}{r}{r}_{\mathrm{a},\mathrm{b}}=\frac{\int P(D|{\rho }_{\mathrm{a}},M^{\mathrm{a}})P\left({\rho }_{\mathrm{a}}|M^{\mathrm{a}}\right)\mathrm{d}{\rho }_{\mathrm{a}}}{\int P(D|{\rho }_{\mathrm{b}},M^{\mathrm{b}})P\left({\rho }_{\mathrm{b}}|M^{\mathrm{b}}\right)\mathrm{d}{\rho }_{\mathrm{b}}},\end{array}$$  (2)

The Bayes factor may be challenging to compute as it requires integrals over $f$ -dimensional space, where $f$ is the number of fit parameters. For this reason, $r$ is often approximated numerically. We used leave-one-out cross-validation with Pareto-smoothed importance sampling^[32^], available in the Python package ArviZ^[33^], to compute Bayes factors.

A Markov chain Monte Carlo (MCMC)^[34^] algorithm, implemented using the emcee library in Python^[35^], was used to perform the inference.

The Bayesian inference procedure is summarized in Figure 1. It was not possible to measure the pre-collision electron spectrum for successful collision shots. Furthermore, the pre-collision spectrum for a given shot may not be well-represented by a summary statistic such as the mean electron spectrum for shots with no colliding laser (null spectra) due to substantial shot-to-shot variation in the electron spectrum. For this reason, null spectra were used to train a neural network, which was used to predict pre-collision electron spectra for successful collision shots. The construction, training and testing of this neural network are discussed by Streeter et al.^[36^]. The neural network consisted of an encoder followed by a translator stage and subsequently a decoder. The former compressed information from four diagnostics (the plasma density, laser energy and pointing and the recombination light emitted by the plasma) into the minimum number of parameters that allowed the key features of the inputs to be re-constructed. The decoder performed a similar function but in reverse, reconstructing a full pre-collision electron spectrum from a small number of inputs. The translator section provided a mapping between the outputs and inputs of the encoder and decoder, respectively. Once trained, the encoder, translator and decoder were used to predict an ensemble of pre-collision electron spectra. This was used to estimate the uncertainty due to the availability of training data for each collision shot analysed using the Bayesian inference procedure. The near-median and standard deviation of the predicted distribution (the former is defined as the spectrum closest in shape to the median spectrum of the distribution) were used to approximate the pre-collision electron spectrum and its uncertainty.

In the following section, curled variables, such as $\mathcal{N}$ , denote post-collision variables and pre-collision variables are italicized. The superscripts ^– and ^′ demarcate properties of the gamma spectrum and measured (as opposed to predicted) observables, respectively.

The pre-collision electron spectra were complex and varied significantly from shot to shot, necessitating an analysis procedure capable of treating arbitrary electron spectra. To this end, a routine was developed that decomposed pre-collision electron spectra into a sum over Gaussian sub-spectra, with mean and standard deviation Lorentz factors, ${\left\langle \gamma \right\rangle}_k$ and ${\Theta}_k$ , respectively, where the subscript $k$ iterates over the number of sub-spectra, ${n}_{\left\langle \gamma \right\rangle }$ . The sub-spectra were fed into the MCMC algorithm, which sampled three inference parameters from their prior distributions, namely the laser ${a}_0$ , the longitudinal offset of the collision from the laser focus, ${Z}_\mathrm{d}$ , and the electron bunch duration, ${\tau}_\mathrm{e}$ , illustrated in Figure 2. Given ${\tau}_\mathrm{e}$ , the electron beam phase-space distribution was obtained by assuming the longitudinal electron beam distribution was Gaussian. The spatio-temporal distribution of the laser intensity at the collision was derived using ${a}_0$ and ${Z}_\mathrm{d}$ , assuming the laser is well-described by the fundamental transverse electromagnetic mode solution to the paraxial Helmholtz equation. The laser phase is described by the Gouy term and the transverse offset of the collision from the laser focus is assumed to be zero. The laser ${a}_0$ and electron beam phase-space distribution comprised the inputs for the forward model, that is, the parameterized model of radiation reaction.

Various collision parameters (see Tables 2–4) affect the energy loss of the electron beam and hence the measured post-collision electron and photon spectra. The interpolation procedure used in the forward models scales poorly with the number of inference parameters, and thus the maximum number of inferred parameters is largely limited by the computational expense associated with an increasing number of interpolation dimensions. For this reason, we selected three parameters to infer and fixed the remaining parameters. The choice of inference parameters was contingent upon a number of factors, including the expected impact of the collision parameters on the post-collision observables, data sub-selection and degeneracy. (That is, two or more parameters that cause similar changes in the post-collision observables. If one such parameter is fixed in the forward model but takes a different value experimentally, the inference procedure can change the value of the free, degenerate parameter to accurately recover the post-collision observables.) As a consequence of the use of degeneracy, the inference parameters should be treated as effective parameters that collectively represent the collision distribution of $\eta$ to the first order, rather than physical characteristics of the electron beam and laser pulse.

Each forward model consisted of five four-dimensional interpolation tables produced using a Monte Carlo code written in C++ (see Appendix B). Three of these tables parameterize the post-collision electron spectrum for a Gaussian pre-collision electron spectrum or sub-spectrum as a Weibull distribution (as indicated by simulations), providing its location, $\mu$ , scale, $\lambda$ , and shape, $\kappa$ :

$$\begin{array}{r}\frac{\mathrm{d}\mathcal{N}}{\mathrm{d}\gamma }=\frac{\kappa }{\lambda }{\left(\frac{\gamma -\mu }{\lambda }\right)}^{\kappa -1}{e}^{-{\left(\frac{\gamma -\mu }{\lambda }\right)}^{\kappa }}.\end{array}$$  (3)

The two remaining interpolation tables returned the photon number, $A$ , and the normalized critical energy, ${\overline{\epsilon}}_\mathrm{c}$ , of the photon spectrum:

$$\begin{array}{r}\frac{\mathrm{d}{\bar{\mathcal{F}}}_{k,l}}{\mathrm{d}\overline{ϵ}}=A{\left(\frac{\overline{ϵ}}{{\overline{ϵ}}_{\mathrm{c}}}\right)}^{-\frac{2}{3}}{e}^{-\frac{\overline{ϵ}}{{\overline{ϵ}}_{\mathrm{c}}}},\end{array}$$  (4)

Pre-collision electron spectra input into the Bayesian inference procedure were decomposed into ${n}_{\left\langle \gamma \right\rangle }$ Gaussian sub-spectra that re-created the full spectra when summed, as illustrated in Figure 3(a):

$$\begin{array}{r}\frac{\mathrm{d}N}{\mathrm{d}\gamma }=\sum _{k=1}^{n_{⟨\gamma ⟩}}{h}_k{e}^{-\frac{{(\gamma -{⟨\gamma ⟩}_k)}^2}{2{\mathrm{\Theta }}_k^2}},\end{array}$$  (5)

Given the longitudinal displacement of the collision from the laser focus, ${Z}_\mathrm{d}$ , and the electron bunch duration, ${\tau}_\mathrm{e}$ , sampled by the Bayesian inference procedure, the phase space of the pre-collision electron spectrum, obtained under the assumption that each sub-spectrum has a Gaussian longitudinal distribution as indicated in Figure 3(b), is given by the following:

$$\begin{array}{r}\frac{{\mathrm{d}}^{2}N}{\mathrm{d}\gamma \mathrm{d}z}=\sum _{k=1}^{n_{⟨\gamma ⟩}}\frac{h_k}{\sqrt{2\pi c^2{\tau }_{\mathrm{e}}^2}}{e}^{-\frac{{(\gamma -{⟨\gamma ⟩}_k)}^2}{2{\mathrm{\Theta }}_k^2}}{e}^{-\frac{{(Z-Z_{\mathrm{d}})}^2}{2c^2{\tau }_{\mathrm{e}}^2}},\end{array}$$  (6)

$$\begin{array}{r}{g}_{k,l}=\frac{h_k}{\sqrt{2\pi c{\tau }_{\mathrm{e}}}}{e}^{-\frac{{(\gamma -{⟨\gamma ⟩}_k)}^2}{2{\mathrm{\Theta }}_k^2}}{e}^{-\frac{{(z_{k,l}-Z_{\mathrm{d}})}^2}{2c^2{\tau }_{\mathrm{e}}^2}}.\end{array}$$  (7)

Together, ${a}_0$ , ${g}_{k,l}$ , ${\Theta}_k$ and ${\left\langle \gamma \right\rangle}_k$ comprise the inputs for the interpolation tables, which constitute the forward model.

The forward model for the post-collision electron spectrum consists of three interpolation tables that each generate one output for the four inputs that parameterize each femto-bunch. The three outputs, namely $\mu$ , $\lambda$ and $\kappa$ , are fed into Equation (3) to obtain each post-collision femto-spectrum, $\frac{\mathrm{d}{\mathrm{\mathcal{F}}}_{k,l}}{\mathrm{d}\gamma }$ . This procedure is illustrated in Figure 4.

As the interpolation tables returned normalized electron spectra, the post-collision electron femto-bunches are weighted by the number of electrons in the corresponding pre-collision femto-bunch and then summed over to obtain the full post-collision electron spectrum, $\frac{\mathrm{d}\mathcal{N}}{\mathrm{d}\gamma }$ :

$$\begin{array}{r}\frac{\mathrm{d}\mathcal{N}}{\mathrm{d}\gamma }=\frac{1}{D_{\mathrm{e}}}\sum _{k=1}^{n_{⟨\gamma ⟩}}\sum _{l=1}^{n_z}{g}_{k,l}{\mathrm{\Theta }}_k\frac{\mathrm{d}{\mathcal{F}}_{k,l}}{\mathrm{d}\gamma },\end{array}$$  (8)

The calculation of the uncertainties on the measured and predicted post-collision electron and photon spectra, and the derivation of the posterior probability, up to a constant, are discussed in this section. The neural network that predicted the pre-collision electron spectra returned ${n}_\mathrm{P}=100$ distinct spectra, $\frac{{\mathrm{d}{N}}_j}{\mathrm{d}\gamma }$ , where $j$ enumerates each prediction. At each iteration of the Bayesian inference procedure, this pre-collision distribution is used to obtain the corresponding post-collision electron spectra for a set of collision parameters, as discussed in Section 3.2. The near-median of the distribution of post-collision spectra, $\frac{\mathrm{d}{\mathcal{N}}_{\mathrm{nm}}}{\mathrm{d}\gamma }$ , approximates the post-collision spectrum. The uncertainty on the predicted post-collision electron spectrum, $\zeta \left(\gamma \right)$ , is the root mean squared (rms) deviation of the distribution of predicted post-collision spectra from the near-median spectrum. The uncertainty in the measured post-collision electron spectrum due to the uncertainties in the magnet, lanex screens and gas jet positions, ${\zeta}^{\prime}\left(\gamma \right)$ , is as follows:

$$\begin{array}{r}{\zeta }^{\prime }\left(\gamma \right)=\frac{\frac{\mathrm{d}{\mathcal{N}}_{\mathrm{nm}}}{\mathrm{d}\gamma }(\gamma +{\sigma }_{\gamma })+\frac{\mathrm{d}{\mathcal{N}}_{\mathrm{nm}}}{\mathrm{d}\gamma }(\gamma -{\sigma }_{\gamma })}{2},\end{array}$$  (9)

$$\begin{array}{rl}& ⟨\log \left(P\left(\frac{\mathrm{d}{\mathcal{N}}^{\prime }}{\mathrm{d}\gamma }\right|M^x\right))⟩=\int (\phantom{\frac{{(\frac{\mathrm{d}{\mathcal{N}}^{\prime }}{\mathrm{d}\gamma }-\frac{\mathrm{d}{\mathcal{N}}_{\mathrm{nm}}}{\mathrm{d}\gamma })}^2}{2({\zeta }^2\left(\gamma \right)+{{\zeta }^{\prime }}^2\left(\gamma \right))}}-\frac{1}{2}\log \left(2\pi ({\zeta }^2\left(\gamma \right)+{{\zeta }^{\prime }}^2\left(\gamma \right))\right)\\ & -\frac{{(\frac{\mathrm{d}{\mathcal{N}}^{\prime }}{\mathrm{d}\gamma }-\frac{\mathrm{d}{\mathcal{N}}_{\mathrm{nm}}}{\mathrm{d}\gamma })}^2}{2({\zeta }^2\left(\gamma \right)+{{\zeta }^{\prime }}^2\left(\gamma \right))})\mathrm{d}\gamma .\end{array}$$  (10)

The uncertainty, ${\overline{\zeta}}_\mathrm{pre}$ , in the predicted post-collision gamma spectrum, $\frac{\mathrm{d}\overline{\mathcal{N}}}{\mathrm{d}\overline{\varepsilon }}$ , due to the uncertainty in the predicted pre-collision electron spectrum is given by the standard deviation of the distribution of predicted post-collision gamma spectra.

A separate Bayesian inference routine was used to fit Equation (4) given the measured gamma spectrometer signal. This procedure yielded ${n}_{\mathcal{Y}}$ -valued distributions of ${\overline{\epsilon}}_\mathrm{c}$ and $A$ that yielded the best fit, as discussed in Appendix A. These values of ${\overline{\epsilon}}_\mathrm{c}$ and $A$ were used to generate a distribution of ${n}_{\mathcal{Y}}$ gamma spectra, $\frac{\mathrm{d}{\mathcal{N}}_v^{\prime }}{\mathrm{d}\overline{\epsilon}}$ , where $v$ , which denotes the vth spectrum in the distribution, has values $0\le v\le {n}_{\mathcal{Y}}$ . The mean, $\left\langle \frac{\mathrm{d}{\overline{\mathcal{N}}}^{\prime }}{\mathrm{d}\overline{\epsilon}}\right\rangle$ , and standard deviation, ${\overline{\zeta}}^{\prime}\left(\overline{\epsilon}\right)$ , of the distribution indicate the most probable gamma spectrum and its corresponding uncertainty. The log likelihood, $\left\langle \log \left(p\left(\left.\frac{\mathrm{d}{\overline{\mathcal{N}}}^{\prime }}{\mathrm{d}\overline{\epsilon}}\right|{M}^x\right)\right)\right\rangle$ , for the gamma spectrum is then as follows:

$$\begin{array}{rl}& ⟨\log \left(p\left(\frac{\mathrm{d}{\bar{\mathcal{N}}}^{\prime }}{\mathrm{d}\overline{ϵ}}\right|M^x\right))⟩\\ & =\frac{1}{n_{\mathrm{P}}}\sum _{j=0}^{n_{\mathrm{P}}}\frac{1}{\mathrm{\Delta }\overline{ϵ}}\int (-0.5\log \left(2\pi {\bar{\zeta }}^2\left(\overline{ϵ}\right)\right)-\frac{{(\frac{\mathrm{d}{\bar{\mathcal{N}}}_j}{\mathrm{d}\overline{ϵ}}-\frac{\mathrm{d}{\bar{\mathcal{N}}}^{\prime }}{\mathrm{d}\overline{ϵ}})}^2}{2{\bar{\zeta }}^2\left(\overline{ϵ}\right)})\mathrm{d}\overline{ϵ},\end{array}$$  (11)

$$\begin{array}{rl}& \log \left(P\left(D_{\mathrm{T}}|M^x\right)\right)\\ & =⟨\log \left(P\left(\frac{\mathrm{d}\bar{\mathcal{N}}}{\mathrm{d}\overline{ϵ}}|M^x\right)\right)⟩+⟨\log \left(P\left(\frac{\mathrm{d}{\mathcal{N}}^{\prime }}{\mathrm{d}\gamma }|M^x\right)\right)⟩,\end{array}$$  (12)

The prior distribution, $P\left({a}_0,{Z}_\mathrm{d},{\tau}_\mathrm{e}\right)$ , is given by the sum over the log priors for ${a}_0$ , $P\left({a}_0\right)$ and ${z}_{k,l}$ , $P\left({z}_{k,l}\right)$ :

$$\begin{array}{r}\log \left(P(a_0,Z_{\mathrm{d}},{\tau }_{\mathrm{e}})\right)=\log \left(P\left(a_0\right)\right)+\log \left(P\left(z_{k,l}\right)\right),\end{array}$$  (13)

$$\begin{array}{rl}& \log \left(P\left(a_0\right)\right)=00.1\le a_0\le 35,\\ & \log \left(P\left(a_0\right)\right)=-\mathrm{\infty }{a}_0<0.1,\\ & \log \left(P\left(a_0\right)\right)=-\mathrm{\infty }{a}_0>35,\end{array}$$  (14)

$$\begin{array}{rl}& \log \left(P\left(z_{k,l}\right)\right)=0-4800/{\omega }_{\mathrm{L}}\le z_{k,l}\le 1440/{\omega }_{\mathrm{L}},\\ & \log \left(P\left(z_{k,l}\right)\right)=-\mathrm{\infty }{z}_{k,l}\le -4800/{\omega }_{\mathrm{L}},\\ & \log \left(P\left(z_{k,l}\right)\right)=-\mathrm{\infty }{z}_{k,l}\ge 1440/{\omega }_{\mathrm{L}}.\end{array}$$  (15)

The quantity optimized by the MCMC, $\log \Big(p\left({M}^x|{D}_\mathrm{T}\Big)\right)$ , was calculated using the following:

$$\begin{array}{r}\log \left(p\left(M^x|D_{\mathrm{T}})\right)\right)=\log \left(P\left(D_{\mathrm{T}}|M^x\right)\right)+\log \left(P(a_0,Z_{\mathrm{d}},{\tau }_{\mathrm{e}})\right).\end{array}$$  (16)

Here, $p\left({M}^x|{D}_\mathrm{T}\right)$ is proportional to the posterior probability, $P\left({M}^x|{D}_\mathrm{T}\right)$ , up to a normalization constant, $P\left({D}_\mathrm{T}\right)$ ; $P\left({D}_\mathrm{T}\right)$ is calculated by integrating the posterior probability over the parameter space. If this integral has no analytic solution and is multi-dimensional, as is the case in this work, it can be extremely challenging to compute $P\left({D}_\mathrm{T}\right)$ accurately. However, as $P\left({D}_\mathrm{T}\right)$ is constant, it does not affect the results of the inference procedure. Therefore, $p\left({M}^x|{D}_\mathrm{T}\right)$ may be used to perform the inference instead without loss of accuracy. The MCMC optimized the log posterior as this enabled the inference procedure to treat extremely small probabilities. The $\log \Big(P\left({M}^x|{D}_\mathrm{T}\Big)\right)$ obtained for a given set of ${a}_0$ , ${Z}_\mathrm{d}$ and ${\tau}_\mathrm{e}$ determined the subsequent region of the parameter space to be explored by the MCMC. This procedure was continued iteratively until the MCMC converged. The degree of convergence and the point at which convergence was reached were calculated as described by Sokal^[38^] and the emcee module documentation^[39^], respectively.

The full forward models were benchmarked using the Monte Carlo code QEDCASCADE^[40^,41^], as illustrated in Figures 5(a) and 5(b).

The interpolation tables that comprised the forward models were produced using a Monte Carlo code written in C++. Appendix B details the computational implementation of each model of radiation reaction.

A number of difficulties arise when implementing Bayesian inference in a context such as this, where many parameters affect the observables, are poorly constrained and their effects on the observables are correlated. If the latter statement holds, a change in one parameter may be partially compensated for by a change in another parameter. If the prior is insufficiently restrictive, many regions of the parameter space may exist that optimize the posterior: a unique solution may not exist. This problem is known as degeneracy. There is also a risk of over-fitting; if the number of inference parameters (i.e., the degrees of freedom of the forward model) is increased, the inference procedure will return progressively larger hyper-volumes of the f-dimensional parameter space which optimize the posterior probability. Without sufficiently constraining priors, additional diagnostics or smaller uncertainties, increasing the number of free parameters may not increase the quantity of meaningful information that can be extracted from the data. Furthermore, if an excessive number of poorly constrained free parameters are used, these parameters may compensate for inaccuracies in the models, allowing any model of radiation reaction to be made compatible with the data. This may be avoided by applying strong priors; however, such priors require measurements of the unknown parameters or must otherwise be physically motivated, and many parameters in collider experiments lack either constraint.

We employ three approaches to address these challenges, as discussed below.

Parameter selection The collision parameters whose expected variations have the highest and lowest impact on the post-collision electron and gamma spectra are identified. The latter are assigned fixed values in the forward models.

Degeneracy Degeneracies between the parameters that have the greatest impact on the observables are identified, if present. These degeneracies are used to further reduce the number of free parameters; two degenerate parameters may be replaced by one, so long as it is able to accurately reproduce the underlying physics of the collision. This inference parameter is then treated as an effective parameter. As such, its inferred value will not reflect a physical property of the electron beam or laser, but rather will reproduce, to first order, the distribution of $\eta$ at the collision. This is because $\eta$ indicates the importance of the correction to the electron energy due to photon emission, which in turn indicates the expected impact of radiation reaction on the post-collision observables. From the definition of $\eta$ given in Section 2, the expression $\eta =2{\tilde{a}}_0\gamma \frac{\lambda_\mathrm{C}}{\lambda_\mathrm{L}}$ may be derived, where ${\lambda}_\mathrm{C}=\frac{h}{m_\mathrm{e}c}$ is the Compton wavelength. Thus, the distribution of the inferred effective laser ${a}_0$ , ${\tilde{a}}_0\left({a}_0,{Z}_\mathrm{d},{\tau}_\mathrm{e}\right)$ , and the electron spectrum can be used to obtain the distribution of $\eta$ at the collision.

Shot selection We propose to perform the analysis on shots selected to minimize the potential impact of collision parameters excluded from the forward model. For these shots, degeneracies between included and excluded parameters enable the inference procedure to compensate for any contributions to the change in electron beam energy and spectral shape from parameters not included in the forward model.

The laser, electron and collision parameters that were measured, estimated or inferred from previous measurements are summarized alongside their assigned values in the forward models in Tables 2–4, respectively. These parameters are illustrated in Figure 2 for clarity. The temporal displacement of the collision from the focus provided in Table 4 combines the shot-to-shot variation in the timing between the two laser pulses with the additional delay due to the unknown injection point of the electron beam. As the electron beam is accelerated to velocities exceeding the group velocity of the laser in the plasma, ${v}_\mathrm{g}$ , if injection occurs earlier, the colliding beam would need to arrive earlier for the electron beam to interact with the peak laser intensity.

Figures 6 and 7 illustrate the effect of varying the laser, electron beam and collision parameters on the mean energy and width of the post-collision electron spectrum, respectively. This allows the collision parameters with the highest impact on the experimental observables to be identified. We have chosen to neglect spatio-temporal coupling terms such as pulse front tilt and chirp in the laser, as these higher-order effects are expected to have a significantly lower impact on the post-collision observables compared to the parameters considered here. The variation in the number of emitted photons, $A$ , with these parameters has not been shown as this trend is merely the inverse of the trend shown in Figure 6. The variation in the critical energy of the photon spectrum, ${\overline{\epsilon}}_\mathrm{c}$ , is not shown as this is much less sensitive to fluctuations in the collision parameters than the properties of the electron spectrum, which are shown.

Figure 6 indicates that, given the expected uncertainties in the collision parameters, the transverse offset has the most significant effect upon the electron beam energy loss, followed by the laser energy, the longitudinal offset of the collision from focus and the source size of the electron beam. By comparison, the effect of changing the electron beam chirp and duration and the laser duration are negligible. The effect of changing the electron beam divergence also appears insignificant; however, this is because the longitudinal offset of the collision from the laser focus was set to 0, and thus the effect of the changing energy-dependent divergence of the electron beam on its transverse size (and hence the average laser intensity during the collision) is negligible.

In Figure 7, the expected transverse jitter also has the highest expected impact on the electron spectrum; however, the most impactful parameters following this are the electron source size, the longitudinal offset and the electron beam duration. In Figures 6 and 7, the quantum-continuous model has been omitted as it exhibits trends similar to the classical model.

The laser, ${a}_0$ , longitudinal offset of the collision from the laser focus, ${Z}_\mathrm{d}$ , and the electron beam duration, ${\tau}_\mathrm{e}$ , were chosen to be free parameters in the inference procedure. As the laser duration and waist at focus remained constant, a change in laser ${a}_0$ indicates a change in the laser energy. The laser ${a}_0$ was chosen as its expected effect on the electron energy losses was larger than any other parameter, barring the transverse jitter. This also allowed an arbitrary number of shots with differing laser energies to be analysed using the same forward model. The longitudinal offset was selected due to the significance of its impact on both the mean and width of the post-collision electron spectrum. Finally, the electron beam duration was chosen over the source size, which has a greater effect on the width of the post-collision electron spectrum, as degeneracy allowed the effect of changing source size to be compensated for by varying ${Z}_\mathrm{d}$ and ${a}_0$ , as will be discussed in Section 5.2. Finally, only the shots with the highest photon yields were analysed to reduce the probability of analysing a collision transversely offset from the laser focus. For this reason, the transverse offset in the forward models was fixed at 0. This is discussed in greater detail in Section 5.2.

A change in source size alters the relative transverse sizes of the laser and electron beam. However, as illustrated in Figure 8, when the longitudinal position of the collision relative to the laser focus varies, the energy-dependent electron beam divergence, which is included in the forward models as a fixed parameter, changes the relative transverse sizes of the electron beam and laser pulse. Thus, if the electron beam source size is fixed in the forward model but varies in the experiment, the transverse distribution of laser intensity experienced by the electron beam during the collision can be recovered in the inference procedure by varying ${a}_0$ and ${Z}_\mathrm{d}$ , if the electron beam divergence in the forward model is non-zero. Here, the relative size of the electron beam and laser is varied by changing ${Z}_\mathrm{d}$ , and corresponding change in laser intensity due to the change in the laser waist at the collision is compensated for by varying ${a}_0$ . This is further supported by Figures 9 and 10, which indicate that degeneracies exist between ${Z}_\mathrm{d}$ and the initial source size for both $\left\langle {\gamma}_\mathrm{f}\right\rangle$ and ${\Theta}_\mathrm{f}$ .

The transverse jitter is by far the most impactful collision parameter. The combination of the large shot-to-shot variation in the electron beam pointing and the steep radial dependence of the intensity of the colliding laser ensure that the electron beam energy loss becomes negligible if the transverse offset is even $0.5\sigma$ from perfect alignment for ${Z}_\mathrm{d}=0$ . As the transverse offset was expected to affect the electron beam energy loss more severely than any other parameter, both in terms of the magnitude of its effect on the electron beam energy loss and the high probability of a transverse misalignment due to large shot-to-shot variations in the electron beam pointing, the shots that produced the highest gamma yield out of a dataset of successful collisions were more likely to be well-aligned transversely than longitudinally. For this reason, we propose that the Bayesian analysis should be applied to shots that produce the highest photon yields. This increases the probability that the transverse offset is small for the selected shots, and can thus be compensated for by exploiting the degeneracy between the laser ${a}_0$ , ${Z}_\mathrm{d}$ and ${\tau}_\mathrm{e}$ and the transverse offset. Note that while this is true for laser pulses with steep transverse intensity profiles (e.g., Gaussian beams), for laser pulses that vary more gradually transversely (i.e., pulses with Airy rings), it may be necessary to include transverse offset as a free model parameter as shot selection will not guarantee a small transverse offset. This would substantially increase the complexity and computational expense of the forward model.

Figures 11 and 12 demonstrate the existence of degeneracies between ${Z}_\mathrm{d}$ and transverse offset, as was the case for source size and ${Z}_\mathrm{d}$ . However, while degeneracies are evident in $\left\langle {\gamma}_\mathrm{f}\right\rangle$ and ${\Theta}_\mathrm{f}$ when these observables are considered separately, for the inference framework to exploit the degeneracies between free and fixed parameters, it must identify a combination of ${a}_0$ , ${Z}_\mathrm{d}$ and ${\tau}_\mathrm{e}$ that reproduces these observables simultaneously.

We investigated whether the Bayesian inference procedure treats free parameters as effective parameters, that is, uses free parameters to reproduce the effects on post-collision observables of fixed parameters that differ from their set values in the forward models. To this end, the procedure was performed on a series of simulated post-collision electron and gamma spectra for each model of radiation reaction. In these simulations, the longitudinal offset, transverse offset or the electron beam source size was varied. Each inference procedure then fitted the simulated data using the corresponding forward model (i.e., the data simulated using the classical model were fitted using the classical inference procedure). For each inference procedure, the parameters fixed in the forward models had the values given in Tables 2–4. The inference results are summarized in Figure 13.

Degeneracies between electron beam duration and longitudinal offset of the collision from focus allow the quantum-stochastic model to recover the post-collision electron and gamma spectra with a total error (which combines the error in the inferred mean and standard deviation of the electron spectrum with the mean photon energy), which is typically at the few percent level, and is always less than $10\%$ . The maximum total error in the classical model is $25\%$ ; however, for the majority of the parameter values considered in Figure 13, the uncertainty is less than $10\%$ . For both models, the uncertainty in the inferred spectra for varying source size (a fixed parameter) is only a few percent higher than the inferences performed for varying longitudinal offset (a fitted parameter), confirming the high degree of degeneracy between these parameters and motivating the decision to include the latter but not the former. The effect of the transverse offset is harder to replicate: this motivates our decision to down-select our data to reduce the probability that a large transverse offset was present for the collisions we analysed using the Bayesian framework. As the gamma radiation yield scales as $\propto {\tilde{a}}_0Q{\gamma}^2$ ^[5^], where $Q$ is the total electron beam charge, by selecting shots with the high gamma yield normalized by $Q{\gamma}^2$ , the probability that the collision is transversely offset from the laser focus for the selected shots is minimized.

A set of inference procedures was run for mono-energetic electron beams with mean energy 1 GeV, in which the transverse offset in the simulated collision was progressively increased, for the following collision parameters: ${a}_0=21.38$ , ${Z}_\mathrm{d}=0$ , ${\tau}_\mathrm{e}=20\kern0.22em \mathrm{fs}$ . It was found that model differentiation was no longer possible for a transverse offset of $1.5{w}_0=3.3\;\unicode{x3bc} \mathrm{m}$ , at which the reduction in the effective collision ${a}_0$ and the spectral broadening induced by the transverse offset rendered model selection infeasible.

Three further test cases were performed to determine whether the inference procedure is able to extract the correct (input) model of radiation reaction and the correct collision parameters for pre-collision electron spectra and uncertainties representative of experimental data, as these spectra are broadband and have complex, non-normal charge distributions. To achieve this, the post-collision observables were simulated for collisions characterized by different parameters in each test case. Three inferences were then performed on the data produced by each simulation – one per model.

For the first test case, described as ‘ideal’, the stochastic model was used to simulate the collision and the simulation parameters have identical values to their fixed counterparts in the forward models. This means that in both the simulation and the forward model, the transverse offset of the collision from the laser focus is zero, the electron beam source size is 0.68 μm, etc. This test case probes the ability of the inference procedure to retrieve the collision parameters and perform model comparison accurately for an ideal scenario, where the collision conditions are fully described by the forward model. The inferred post-collision electron spectra and photon spectrometer responses for this test case are compared to the simulated data in Figure 14.

In Figure 15, the quantum-stochastic and quantum-continuous models retrieve $\left\langle {\tilde{a}}_0\right\rangle$ , $\left\langle \eta \right\rangle$ and ${\sigma}_{\eta }$ within $1\sigma$ of the input parameters. None of the models infer ${\tau}_\mathrm{e}$ within $1\sigma$ of the input value. For the classical and quantum-continuous models, this indicates that ${\tau}_\mathrm{e}$ has been treated as an effective parameter. Rather than reflecting the true electron beam duration, the inferred value of ${\tau}_\mathrm{e}$ allows the distribution of effective ${a}_0$ and $\eta$ to be most closely reproduced. Of the three models, the quantum-stochastic model infers the value of ${\tau}_\mathrm{e}$ most accurately, but fails to recover the correct value within $1\sigma$ , indicating that the uncertainty in the pre-collision electron spectrum and the lack of sensitivity of the post-collision electron spectrum to small variations in ${\tau}_\mathrm{e}$ inhibit accurate inferences of this parameter. The classical and quantum-continuous inferences both retrieve ${\tau}_\mathrm{e}$ considerably greater than the true value, thereby increasing the range of ${\tilde{a}}_0$ with which the electron beam interacts and hence the range of energy losses it experiences. This allows the classical and quantum-continuous models to reproduce the spectral broadening induced by the quantum-stochasticity inherent in the quantum-stochastic model, which the classical and quantum-continuous models do not predict. As the classical model predicts higher energy losses than the quantum-stochastic model for equivalent collision parameters, it infers a lower $\left\langle {\tilde{a}}_0\right\rangle$ to produce comparable energy losses.

The $\left\langle {\tilde{a}}_0\right\rangle$ and standard deviation ${\tilde{a}}_0$ (the latter results from the spatial and temporal overlap of the electron beam and laser) are reflected in the mean, $\left\langle \eta \right\rangle$ , and standard deviation, ${\Sigma}_{\eta }$ , of the distribution of $\eta$ , which characterizes the collision.

When comparing the quantum-stochastic and classical models, ${r}_{\mathrm{qs},\mathrm{cl}}=1.2$ , while the quantum-continuous and classical models yield ${r}_{\mathrm{qc},\mathrm{cl}}=1.3$ , and the quantum-stochastic and quantum-continuous model comparison gives ${r}_{\mathrm{qc},\mathrm{qs}}=1.1$ , indicating there is insufficient evidence to differentiate among the quantum-stochastic, quantum-continuous and classical models for the simulated collision conditions, given the uncertainties on the predicted pre-collision electron spectrum and the measured post-collision electron and gamma spectra.

For the second test case, the results of which are shown in Figure 16, the classical model was used to simulate a collision that was transversely offset from the laser propagation axis by 1.05 μm. Three inference procedures, one for each radiation reaction model, were performed on the simulated electron and photon spectra. The presence of a finite transverse offset between the electron beam and laser focus at the collision induces electron spectral broadening, which resembles the spectral broadening predicted by the quantum-stochastic model. This test case illustrates the extent to which the inference procedure is able to perform model selection accurately (i.e., select the classical model) if the collision parameters with fixed values in the forward models differ from those values. This is particularly pertinent if the classical model is accurate and the additional collision parameters induce spectral broadening, an effect also predicted by the quantum-stochastic model. This test case also indicates whether and how inference procedures use degeneracy to compensate for collision parameters that differ from the fixed values in the forward model.

The inferred and input collision parameters for the second test case in Figure 17 indicate that only the classical model is able to infer the mean effective ${a}_0$ of the collision and the mean $\eta$ within $1\sigma$ . However, all models, including the classical model, over-estimate the electron beam duration and ${\sigma}_{\eta }$ . This contrasts with the first test case in which the correct model (i.e., the quantum-stochastic model) also inferred ${\Sigma}_{\eta }$ within $1\sigma$ of its true value. In the second test case, all models treat ${\tau}_\mathrm{e}$ as an effective parameter, stretching the electron beam longitudinally to replicate (in so far as possible) the broad distribution of ${a}_0$ with which the electron beam interacted due to the finite transverse offset. This is also reflected in the overly large ${\Sigma}_{\eta }$ inferred. For the second test case (see Figure 16), each of the models appears to infer the post-collision electron and gamma spectra with comparable accuracy. This is substantiated by the Bayes factors; we obtain ${r}_{\mathrm{cl},\mathrm{qs}}=1.0$ , while ${r}_{\mathrm{cl},\mathrm{qc}}=1.4$ and ${r}_{\mathrm{qs},\mathrm{qc}}=1.4$ , and thus there is insufficient evidence to favour any of the models.

In the third test case, the collision was simulated using the quantum-stochastic model and was transversely offset by 3.2 μm from the laser propagation axis. This allowed the model selection capability of the Bayesian analysis to be verified for a collision with parameters that differed from their fixed values in the forward model, in which the quantum-stochastic model was used to produce the test data. The results for the third test case are provided in Figures 18 and 19. For the third test case, the Bayes factor obtained when comparing the quantum-stochastic and classical models, ${r}_{\mathrm{qs},\mathrm{cl}}=1.8$ , the quantum-continuous and classical models yield ${r}_{\mathrm{qc},\mathrm{cl}}=1.7$ while a comparison of the quantum-stochastic and quantum-continuous inferences yields ${r}_{\mathrm{qc},\mathrm{qs}}=1.0$ , and thus there is insufficient evidence to differentiate among the quantum-stochastic, quantum-continuous and classical models.

Each test case yielded $0.33, indicating model selection could not be performed for a single shot given the collision conditions and uncertainties considered for each test case. However, for both test cases in which the quantum-stochastic model was used to generate the test data, ${r}_{\mathrm{qs},\mathrm{cl}}>1.0$ , indicating model selection would be feasible if evidence were combined across multiple shots. For the ideal test, assuming model evidence is consistent for each shot, seven shots would be required to allow model differentiation.

For the test case in which the classical model was used to generate the test data, ${r}_{\mathrm{cl},\mathrm{qs}}=1.0$ and ${r}_{\mathrm{cl},\mathrm{qc}}>1.0$ , indicating that if the collision parameters induce spectral broadening, the stochastic model would not be unduly favoured over the classical model.

We have ascertained that the quantum-stochastic model retrieves $\left\langle {\tilde{a}}_0\right\rangle$ , $\left\langle \eta \right\rangle$ and ${\Sigma}_{\eta }$ within $1\sigma$ of the input parameters when this model is used to generate the simulated test data and the forward model accurately describes the full complexity of the collision. As illustrated in the classical test case, if additional collision parameters that are not present in the forward model are present in the simulated data, the correct model (i.e., the classical model) infers $\left\langle {\tilde{a}}_0\right\rangle$ and $\left\langle \eta \right\rangle$ within $1\sigma$ of the simulation values, but significantly over-estimates ${\tau}_\mathrm{e}$ and hence ${\sigma}_{\eta }$ . This indicates that while the distributions of $\eta$ and $\left\langle {\tilde{a}}_0\right\rangle$ are inferred accurately to the first order, there are insufficient degrees of freedom in the models to accurately infer higher-order moments (i.e., ${\sigma}_{{\tilde{a}}_0}$ , ${\sigma}_{\eta }$ ) if collision parameters that are fixed in the forward models differ significantly from these values.

We have developed a novel Bayesian framework that infers values of unknown collision parameters and predicts corresponding experimental observables for the classical, quantum-continuous and quantum-stochastic models of radiation reaction. We identify challenges associated with the application of a Bayesian approach to this problem, such as over-fitting and insufficiently constraining priors. We address these issues by down-selecting the number of free parameters and the data to be analysed and by exploiting degeneracies between free and fixed parameters. This has motivated the choice of ${a}_0$ , ${Z}_\mathrm{d}$ and ${\tau}_\mathrm{e}$ as fitting parameters and the decision to fix the remaining parameters.

We demonstrate that the Bayesian framework consistently infers $\left\langle {\tilde{a}}_0\right\rangle$ , $\left\langle \eta \right\rangle$ , ${\sigma}_{\eta }$ accurately (within $1\sigma$ ) for the highest-performing model. The $1\sigma$ criterion has been used to provide a rule-of-thumb indication of the accuracy of the inferred parameters. Quantitative comparisons of the relative performance of these models for a single shot yield insubstantial evidence in favour of the correct model over other models, indicating that while model discrimination may not be feasible at the single-shot level, this may be accomplished by combining evidence across multiple shots. For the test cases we have considered, model evidence combined over seven shots or fewer would yield substantial evidence favouring the correct model, assuming each shot produces the same Bayes factor. We find that when the transverse offset differs from its fixed value in the forward models, the standard deviations of the distributions of $\left\langle {\tilde{a}}_0\right\rangle$ and $\eta$ are inferred less accurately; however, accurate model selection with this framework is still robust for transverse offsets up to 3.3 μm, roughly 2.4 times the laser waist, for the electron beam and laser parameters considered.

The inclusion of free parameters such as the transverse offset would be facilitated if the computational expense (and runtime) of the Bayesian analysis were reduced. This could be accomplished if the electron spectra were mono-energetic and if strong priors could be applied to restrict the available parameter space and avoid over-fitting. Strong priors would also facilitate the inference of parameters that accurately represent laser and electron beam parameters, rather than effective parameters that reproduce the collision distribution of $\eta$ . Another route to increase the number of free parameters in the forward model or decrease the computational cost of the procedure could involve substituting the interpolation with a machine learning technique (e.g., a neural network) with a lower associated computational cost. However, it should be noted that a large number of simulations would be required to train a neural network to accurately predict electron and photon spectra for many possible combinations of different collision parameters, which could be computationally expensive to generate. In addition, care would need to be taken to ensure the training dataset fully represented the neural network inputs during the inference procedure.

From an experimental perspective, improved laser stability would facilitate a greater fraction of collisions with good spatial-temporal overlap between the electron beam and colliding laser, increasing the number of collisions for which quantum effects are expected to be substantial.

We anticipate that Bayesian inference will prove to be a powerful analysis tool for the interpretation of future strong-field QED experiments involving colliding lasers and particle beams, and have demonstrated the feasibility and utility of such an analysis for an all-optical radiation reaction experiment.