1. Introduction
Strong-field quantum electrodynamics (SFQED) studies the interaction between matter and intense electromagnetic fields. In recent years, there has been a growing interest in this area due to the availability of high-intensity laser sources, which enable the exploration of novel physical phenomena, with experiments being planned for the near future at ELI Extreme Light Infrastructure, Apollon (Papadopoulos, Zou & Le Blanc Reference Papadopoulos, Zou and Le Blanc2016), CoReLS (Yoon et al. Reference Yoon, Jeon, Shin, Lee, Lee, Choi, Kim, Sung and Nam2019), FACET-II (Meuren et al. Reference Meuren, Bucksbaum, Fisch, Fiúza, Glenzer, Hogan, Qu, Reis, White and Yakimenko2020), LUXE (Abramowicz et al. Reference Abramowicz2019, Reference Abramowicz2021), XCELS (Mukhin et al. Reference Mukhin, Soloviev, Perevezentsev, Shaykin, Ginzburg, Kuzmin, Mart'yanov, Shaikin, Kuzmin, Mironov, Yakovlev and Khazanov2021), ZEUS, NSF-OPAL, HIBEF, among others. Many of the proposed experimental set-ups feature scattering of intense, focused laser pulses with either relativistic electron beams or high-energy photons.
Some experimental evidence for radiation reaction has been observed (Cole et al. Reference Cole2018; Poder et al. Reference Poder2018; Los et al. Reference Los2024). However, additional theoretical and numerical studies are needed to understand the effects of ‘non-ideal’ experimental conditions such as laser jitter, probe beam emittance, synchronisation of the collision, etc. In the future, better control over these parameters will allow precision studies of the radiation reaction, that is, the recoil on the charged particles that emit high-energy photons and electron–positron production in the lab, among other processes. These can then be applied in physics models of plasmas in extreme astrophysical environments, such as black holes and pulsar magnetospheres (Timokhin Reference Timokhin2010; Medin & Lai Reference Medin and Lai2010; Uzdensky & Rightley Reference Uzdensky and Rightley2014; Cruz et al. Reference Cruz, Grismayer, Chen, Spitkovsky and Silva2021; Schoeffler et al. Reference Schoeffler, Grismayer, Uzdensky and Silva2023).
While current laser technology allows us to test strong-field plasma physics in a ‘semi-classical’ regime, next-gen laser facilities will require exploring the fully non-perturbative quantum dynamics of fermions, high-energy photons and the laser field. This regime is expected to require first-principles simulation techniques, thus putting constraints on the more standard particle-in-cell simulations, even if including Monte Carlo routines for the quantum processes. It is thus important to begin translating known simulation methods into the quantum computing framework.
Quantum computing has the potential to handle the complexity of the many-body-physics dynamics in these extreme plasmas. Recently, the plasma physics community has started to adapt standard plasma set-ups and theory to the quantum algorithmic framework (Dodin & Startsev Reference Dodin and Startsev2021; Amaro & Cruz Reference Amaro and Cruz2023; Joseph et al. Reference Joseph, Shi, Porter, Castelli, Geyko, Graziani, Libby and DuBois2023). More specifically, the intersection between plasma physics, quantum field theory and quantum computing has only recently started to be explored (Shi, Qin & Fisch Reference Shi, Qin and Fisch2021).
In this work, we aim to address some of these questions, first by deriving approximate particle distributions and benchmarking with the particle-in-cell code OSIRIS (Fonseca et al. Reference Fonseca, Sloot, Hoekstra, Tan and Dongarra2002) that describes QED effects with a Monte Carlo routine to simulate quantum events and has already been tested in this regime (Vranic et al. Reference Vranic, Grismayer, Fonseca and Silva2016) and later by applying a hybrid quantum algorithm to the simulation of a Fokker–Planck equation relevant for SFQED. We choose the simplest possible field configuration, a strong constant magnetic field background, that could be generalised in the future. We start with a relativistic electron beam, which propagates and loses energy to radiation. We follow the evolution of the electron distribution function over time, accounting for the stochastic nature of the quantum emission.
While prior work (e.g. Kubo et al. Reference Kubo, Nakagawa, Endo and Nagayama2021) suggests that hybrid quantum algorithms could provide sampling efficiency benefits for more complex stochastic differential equations (SDEs) – such as those with high dimensionality or intricate observables relevant to finance and physics – even within the NISQ era, this implementation does not outperform classical methods at this point. As such, our primary aim is not to claim a quantum advantage, but to present a proof of concept demonstrating that quantum hardware can, in principle, be used to simulate SDEs with applications to plasma physics.
This manuscript is structured as follows. In §2, we describe the set-up of an electron beam propagating perpendicularly to a strong magnetic field and radiating energy in the form of photons in the so-called semi-classical regime. We introduce the Fokker–Planck equation and derive solutions for the evolution of the distribution function and the first two energy moments, with details in Appendix B. In §3, we describe the quantum variational simulation of the Fokker–Planck equation, the choice of ansatz and the numerical evolution of the parameters using variational quantum imaginary time evolution (VarQITE). Conclusions are given in §4.
2. Fokker–Planck equation for quantum radiation reaction
In this section, we introduce the main parameters for the physical set-up, the Fokker–Planck equation describing the interaction of electrons with an intense magnetic field, and approximate formulae for the first two moments of the electron distribution functions. The importance of quantum effects of relativistic particles in strong fields is controlled by the quantum nonlinearity parameters, which for leptons and photons are given by
\begin{equation} \chi \equiv \dfrac {\sqrt { (p_\mu \,F^{\mu \nu })^2 }}{E_S \,m c}, \quad \chi _\gamma \equiv \dfrac {\sqrt { (\hbar k_\mu \,F^{\mu \nu })^2 }}{E_S \,m c}, \end{equation}
where
$m$
is the electron mass,
$c$
is the speed of light in vacuum,
$E_S = m^2 c^3/(e\hbar )= 1.32 \times 10^{18}$
V m
$^{-1}$
and
$B_S = E_S/c= 4.41 \times 10^{9}$
T represent the critical Schwinger electric and magnetic fields,
$e$
is the elementary electric charge,
$\hbar$
is the reduced Planck constant,
$F_{\mu \nu }$
is the electromagnetic field tensor, and
$p_\mu , k_\mu$
are the leptonic and photonic 4-momenta. In the case of a relativistic electron moving perpendicularly to a strong, uniform and constant magnetic field, we have
\begin{equation} \chi = \gamma B/B_S. \end{equation}
The electron distribution function
$f=\mathrm{d}N/\mathrm{d}\gamma$
, that is, the number of electrons
$\mathrm{d}N$
per interval of energy
$\mathrm{d}\gamma$
, evolves through a Fokker–Planck equation describing stochastic energy losses
\begin{equation} \frac {\partial f(t, \gamma )}{\partial t}=\frac {\partial }{\partial \gamma }\left [-(\mathcal{A} f)+\frac {1}{2} \frac {\partial }{\partial \gamma }\left (\mathcal{B} f\right )\right ]. \end{equation}
This is the partial differential equation that will be simulated through a quantum algorithm in a later section. The scalar drift and diffusion space-dependent terms can be approximated in the
$\chi \ll 1$
regime as
\begin{equation} \begin{split} \mathcal{A} \sim \dfrac {2}{3} \dfrac {\alpha m c^2}{\hbar } \chi ^2 = a\, \gamma ^2, \quad\mathcal{B} \sim \dfrac {55}{24 \sqrt 3} \dfrac {\alpha m c^2}{\hbar } \gamma \, \chi ^3 = \dfrac {b^2}{2} \, \gamma ^4, \end{split} \end{equation}
where
$a \equiv 2\alpha k^2/(3 \tau _c)$
,
$\alpha \equiv e^2/(\hbar c)$
is the fine structure constant and
$b \equiv \sqrt {55\alpha k^3 /(12\sqrt {3}\tau _c)}$
, with Compton time
$\tau _c \equiv \hbar /(mc^2)$
and normalised magnetic field
$k \equiv B/B_S$
(Neitz & Di Piazza Reference Neitz and Di Piazza2013; Vranic et al. Reference Vranic, Grismayer, Fonseca and Silva2016). In general, there are no closed-form solutions for the electron distribution function valid in all regimes of quantum nonlinearity
$\chi$
. Further details on the theory of quantum radiation reaction are presented in Appendix B.
2.1. Evolution of the electron distribution function
In this section, we use the classical approach of simulating the particle motion coupled with Monte Carlo routines to account for the hard-photon emission and quantum radiation reaction. As such, it can describe the evolution of the electron distribution, including the energy loss and quantum stochasticity-induced widening of the spectrum.
In OSIRIS, the Monte Carlo (MC) subroutine is added to the particle-in-cell algorithmic loop between the steps of integration of equations of motion (e.g. the Lorentz force) and the current deposition on the grid. First, the total probability of the cross-section for the MC process to occur within the simulation time step is computed for all particles. If this probability is above a random uniformly generated value in the interval
$[0,1]$
, then the process is assumed to have occurred. For the particles in which this happens, a second random number is generated to sample the energy of the generated particle from the corresponding differential cross-section. The new particles with their respective energies will then be responsible for updating the currents on the grid. The PIC loop then proceeds as usual.
In the
$\chi \ll 1$
regime, an initially narrow distribution with energy
$\gamma _0$
retains an approximately Gaussian shape throughout the interaction. Following the approach of Torgrimsson (Reference Torgrimsson2024a
,Reference Torgrimsson
b
) and Blackburn (Reference Blackburn2024), which is based on a perturbative expansion in
$\chi$
, the mean energy
$\mu = \langle \gamma \rangle = \int \gamma f\,\mathrm{d}\gamma$
and energy spread
$\sigma ^2 = \int (\gamma -\mu )^2 f\,\mathrm{d}\gamma$
in a constant magnetic field can be obtained as
\begin{equation} \dfrac {\mu }{\gamma _0} \sim \dfrac {1}{1+2 R_c t/3} + \dfrac {165\chi _0}{8 \sqrt {3}(1+2/3R_ct)^2} \log \left (1+\dfrac {2R_ct}{3}\right ), \end{equation}
\begin{equation} \dfrac {\sigma ^2}{\gamma _0^2} \sim \dfrac {\sigma _0^2 + 55 R_c \chi _0 t/(24 \sqrt {3})}{(1+2R_c t/3)^4} ,\\[10pt]\end{equation}
where
$R_c \equiv \alpha \,b_0 \,\chi _0$
is the classical radiation reaction parameter,
$b_0 \equiv e B/(m\,\omega _c)$
is an adimensional normalised magnetic field and
$\omega _c \equiv e B/(m \gamma _0)$
is the synchrotron frequency. In our simulations, the time is normalised as
$t\rightarrow t\, \omega _c$
, and the average electron energy and normalised magnetic field are
$\gamma _0=b_0=1800$
. From (2.6), the maximum energy spread occurs approximately at
$t\sim 1/(2 R_c)$
for
$\sigma _{max} \sim 3^{5/4} \sqrt {55}/64 \,\gamma _0 \sqrt {\chi _0}$
, which recovers the scaling derived by Vranic et al. (Reference Vranic, Grismayer, Fonseca and Silva2016). To test the validity of these analytical expressions, we run 1-D3V simulations with an electron beam moving perpendicular to the constant, uniform magnetic field vector. The time step is
$\mathrm{d}t = 0.001 \,\omega _c^{-1}$
, the cell size
$\mathrm{d}x=0.049\,c/\omega _c$
, and the normalising frequencies and external magnetic field values are chosen to enforce
$\chi _0=\{10^{-3},10^{-2}, 10^{-1}\}$
.
Figure 1. Evolution of the moments of the distribution function of electrons in a constant magnetic field. (a) Average energy and (b) spread for zero initial spread
$\sigma _0=0$
(monoenergetic beam of electrons). Dashed lines , theory; coloured lines, OSIRIS simulation results.
In figure 1, we show results from OSIRIS MC simulations for different values of the external magnetic field (different values of initial, average
$\chi _0$
), and the moments from (2.5) and (2.6). The latter expressions capture the simulation results well, with some loss of accuracy for
$\chi _0=10^{-1}$
.
In figure 2, we show snapshots of Gaussian distribution functions using the parameters
$(\mu ,\sigma )$
from (2.5) and energy histograms from the OSIRIS MC simulations. Although for
$\chi _0=10^{-2}$
there is good agreement with simulation results, it is clear that in the higher
$\chi _0=10^{-1}$
regime, the validity of the Gaussian functional approximation is lost, even though the moments
$(\mu , \sigma )$
predicted analytically are close to what was obtained in the simulation. This deviation is due to an increased skewness of the distributions. These results will serve as a ‘ground truth’ for testing the quantum algorithm presented in the following section.
Figure 2. Snapshots of the distribution functions, taking the analytical formulas for the mean energy and energy spread ((2.5) and (2.6)) for
$\chi _0=\{10^{-2}, 10^{-1}\}$
. Dashed lines, Gaussian approximation; coloured lines, OSIRIS simulation results. For
$\chi _0=10^{-1}$
, there is a visible deviation from the Gaussian approximation.
3. Variational quantum simulation
In this section, we introduce the variational quantum simulation (VQS) method employed in our study. Quantum computing holds the potential to revolutionise simulations of complex systems, encompassing both quantum many-body problems and certain classical physics phenomena. The fundamental unit of quantum computation is the qubit, which, unlike a classical bit, can exist in a superposition of states described by
$|\psi \rangle = a_0|0\rangle + a_1|1\rangle$
, where
$a_0$
and
$a_1$
are complex amplitudes satisfying the normalisation condition
$|a_0|^2 + |a_1|^2 = 1$
(Nielsen & Chuang Reference Nielsen and Chuang2010). When multiple qubits are combined into a register, their joint state can exhibit entanglement – a uniquely quantum mechanical correlation that leads to an exponentially large state space, challenging to simulate efficiently on classical computers. In this work, we use the big-endian convention – the left-most qubit is the most significant.
Currently, quantum computers are in the noisy intermediate-scale quantum (NISQ) era (Preskill Reference Preskill2018), characterised by devices that contain a moderate number of qubits, but are susceptible to errors and decoherence. These limitations restrict the depth and complexity of quantum circuits that can be reliably executed, posing significant challenges for implementing algorithms requiring long coherence times and high gate fidelities.
To address these challenges, variational quantum circuits (VQCs) have emerged as a promising approach suitable for NISQ devices. VQCs are hybrid quantum-classical algorithms that employ parametrised quantum circuits optimised using classical optimisation routines to minimise a cost function, typically the expectation value of an observable. This method reduces the required circuit depth by offloading part of the computational workload to classical processors, making it more practical for current quantum hardware. Variational quantum algorithms have found applications across various fields, including quantum chemistry (Peruzzo et al. Reference Peruzzo, McClean, Shadbolt, Yung, Zhou, Love, Aspuru-Guzik and O’Brien2014), material science (Liu, Li & Yang Reference Liu, Li and Yang2024), biology (Robert et al. Reference Robert, Barkoutsos, Woerner and Tavernelli2021) and others (Farhi, Goldstone & Gutmann Reference Farhi, Goldstone and Gutmann2014). By carefully designing the variational circuits and selecting appropriate ansätze, VQCs leverage the expressive power of quantum systems while operating within the practical constraints of NISQ devices. Extracting meaningful information from the exponentially large quantum state space necessitates meticulous selection and measurement of observables, as quantum measurements yield probabilistic outcomes that collapse the superposed state. The variational approach thus provides a flexible and adaptive framework for quantum simulation, contributing significantly to the advancement of quantum computing applications in science and engineering.
In plasma physics, several phenomena of interest are intrinsically dissipative. Since dissipation is an irreversible process, this is challenging to model on quantum computers. Engel, Smith & Parker (Reference Engel, Smith and Parker2019) studied Landau damping on a quantum framework using a linearised version of the Vlasov equation. Despite the name ‘damping’, the energy of the system is not lost, but rather transferred from the electric field to the particle distribution function in a reversible/unitary manner. Vissers & Bouten (Reference Vissers and Bouten2019) studied a ‘quantum stochastic’ process of a laser driven two-level atom interacting with an electromagnetic field. However, a procedure to generalise beyond two-level systems is not provided. Kubo et al. (Reference Kubo, Nakagawa, Endo and Nagayama2021) proposed a partial/stochastic differential equation (PDE/SDE) solver of the Fokker–Planck equation in the form of an Itô-process in a quantum variational framework. This approach was then generalised by Alghassi et al. (Reference Alghassi, Deshmukh, Ibrahim, Robles, Woerner and Zoufal2022) through the Feynman−Kac formula, which unifies the heat, Schrödinger, Black−Scholes, Hamilton−Jacobi and Fokker–Planck equations under this formalism.
The expressibility of variational quantum circuits is one of the most important figures of merit in the field of VQS. It is defined as the ability of the variational quantum circuit to produce a variety of quantum wavefunctions – the higher this value is means the larger the fraction of the Hilbert space is that is accessible through the circuit. This metric has been studied in detail, and heuristics on which architectures one can reach higher expressibility have been discussed (Sim, Johnson & Aspuru-Guzik Reference Sim, Johnson and Aspuru-Guzik2019). However, few proofs or universal rules have been derived so far. In the case of Kubo et al. (Reference Kubo, Nakagawa, Endo and Nagayama2021), the authors suggest using an ansatz with only
$R_Y$
and CNOT quantum gates, which permits only real-valued amplitudes of the wavefunction, but it also allows for negative values. Dasgupta & Paine (Reference Dasgupta and Paine2022) suggest an architecture to enforce even-symmetry of the real-valued wavefunction around the middle of the computational basis representation. Endo et al. (Reference Endo, Sun, Li, Benjamin and Yuan2020) based the motivation for using compact variational circuits on the intuition that the dynamics of the physical system only span low energy states and therefore is limited to a small portion of the entire Hilbert space.
In this work, we consider wave/distribution functions that are well localised in space, but can have a mean position/energy that can be off-centred, that is, not located at the middle of the computational basis
$|01\ldots 1\rangle , |10\ldots 0\rangle$
. For this, we use a quantum circuit which simulates the advection equation based on the finite differences operator as presented by Sato et al. (Reference Sato, Kondo, Hamamura, Onodera and Yamamoto2024). This operation can be used to shift the entire wavefunction, which allows the exploration of a larger space of states.
In figure 3, we show the hybrid quantum optimisation loop. First, at a time step
$t$
, a quantum wavefunction is produced by the quantum circuit, where some of the quantum unitary operations are parametrised by values
$\theta (t)$
. Then, a cost function
$\langle C \rangle$
is measured from the quantum circuit and minimised/updated through an algorithm run on a classical device. The new parameters are then fed-back to the quantum circuit to produce the wavefunction for the next time step
$\theta (t+\mathrm{d}t)$
. The finite-difference approximation to the Fokker–Planck operator is of the trinomial type, where any grid cell at any time step is only influenced by its two close neighbouring cells. Higher-order finite-difference schemes are possible to apply, but require more complex algorithms.
Figure 3. Variational quantum algorithm. From the quantum processing unit (QPU), we measure a cost function, from which the variational parameters are optimised in the CPU.
Figure 4. (a) Variational ansatz used in this work. Only four qubits are shown for simplicity. The blue box represents a variational block that can be repeated
$k$
times. The red box represents the enforcing of even-symmetry on the wavefunction. In the last step, the advection operator is applied to displace the wavefunction from the middle of the grid. (b) Sketch of the typical evolution of the wavefunction.
3.1. Variational ansatz
In figure 4, we show the ansatz used in this work. In the first variational block (blue), we produce a real-valued wavefunction using parameters
$\boldsymbol{\theta }$
and gates
$\mathrm{CNOT}-R_Y$
in a ring structure (the qubit
$n-1$
connects with the first qubit). This block can be repeated to increase the expressibility of the ansatz. The second block (red) consists of enforcing even-symmetry on the wavefunction; that is, the amplitudes have a mirror symmetry around half of the computational basis
$\psi _{0111} = \psi _{1000}, \psi _{0000}=\psi _{1111}$
. If, instead, we apply an angle
$-\pi /2$
on the
$R_Y$
gate acting on the last qubit, then the wavefunction becomes odd instead of even-symmetric. At
$t=0$
, we fit these parameters so that the wavefunction has an approximate Gaussian shape with a prescribed standard deviation
$\sigma$
. In this way, the parameters
$\boldsymbol{\theta }$
are implicit functions of
$\sigma$
.
In the final step of the quantum circuit, we apply the advection operator
$e^{-i \theta _\mu \partial _x }$
, which is a unitary operation, and which translates the entire wavefunction by some quantity
$\theta _\mu$
. Sato et al. (Reference Sato, Kondo, Hamamura, Onodera and Yamamoto2024) provide a quantum circuit for the finite-difference version of the first order of
$\partial _x \sim (s_+-s_-)/2$
, where
$s_+ = s_-^\dagger$
is the operator that shifts the wavefunction one cell to the right (left) with a periodic boundary condition. Here, we use the second order of this finite-difference operator to enforce higher fidelity of the translation.
3.2. Variational quantum simulation of the Fokker–Planck equation
In figure 5, we show the complete evolution of the distribution functions for
$\mu _0=1800$
,
$\sigma _0=90$
from OSIRIS Monte Carlo simulations, following the same set-up parameters as Niel et al. (Reference Niel, Riconda, Amiranoff, Duclous and Grech2018), which we take as being the ‘ground-truth’ benchmark. In these simulations, the electrons can be seen to cool down quicker for higher
$\chi _0$
, and the distribution functions occupy different fractions of the energy domain throughout the evolution, thus having different requirements on the resolution for the energy grid.
Figure 5. Evolution of electron distribution functions from OSIRIS simulations for different initial average
$\chi _0$
values and spread
$\sigma _0=90$
. In all cases, the maximum simulation time is
$20\,\omega _c^{-1}$
.
The VQS algorithm maps the dynamics of the quantum state, (2.3), to those of the variational parameters
$\theta (t)$
of the ansatz. The mapping is performed using McLachlan’s variational principle (MVP) (McLachlan Reference McLachlan1964; Endo et al. Reference Endo, Sun, Li, Benjamin and Yuan2020). Using the same notation as Kubo et al. (Reference Kubo, Nakagawa, Endo and Nagayama2021), the unitary quantum circuit from figure 4 produces a wavefunction
$|v\{\boldsymbol{\theta }(t)\}\rangle$
. However, diffusive processes generally do not preserve the
$L_2$
norm, that is,
$\partial _t \int |v|^2\,\mathrm{d}x \neq 0$
. Therefore, it is useful to define an overall variational wavefunction
$|\tilde {v}\{\boldsymbol{\theta }(t)\}\rangle = \alpha (t)\,|v\{\boldsymbol{\theta }(t)\} \rangle$
, where
$\alpha (t)$
is a (classical) normalisation parameter. If the temporal resolution is sufficiently high and the ansatz sufficiently expressive, the MVP method evolves the variational parameters such that the correct dynamics is enforced through
\begin{equation} \min _{\theta (t)} \bigg \| \frac {\mathrm{d}}{\mathrm{d} t}|\tilde {v}\{\boldsymbol{\theta }(t)\}\rangle -\hat {L}\,|\tilde {v}\{\boldsymbol{\theta }(t)\}\rangle \bigg \|, \end{equation}
where
$\hat {L}$
is a linear operator that generates the evolution of the system. At each time step, the following matrix equation needs to be solved on a classical computer
\begin{equation} M_{k,j} \,\dot {\theta }_j = V_k, \end{equation}
where
\begin{equation} M_{k,j} \equiv \mathrm{Re} \left ( \dfrac {\partial \langle \tilde {v} \{ \theta (t) \} | }{\partial \theta _k} \dfrac {\partial | \tilde {v} \{ \theta (t) \} \rangle }{\partial \theta _j} \right ), \quad V_k \equiv \mathrm{Re} \left ( \dfrac {\partial \langle \tilde {v} \{ \theta (t) \} | }{\partial \theta _k} \hat {L} | \tilde {v} \{ \theta (t) \} \rangle \right ) \end{equation}
are quantities that can be obtained efficiently through dedicated quantum circuits. One advantage of this method, if run on a quantum device, is that one does not need to have access to the full amplitudes of the wavefunction. Instead, only the variational parameters need to be tracked and their evolution is obtained through efficient measurement of a few observables. Whereas the
$M_{k,j}$
only depend on the chosen ansatz, the
$V_k$
elements depend on both the ansatz and the particular Fokker–Planck equation operators (in this case, (2.3)). In this work, we classically simulate the retrieval of these quantities through numerical differentiation and integration of the wavefunction using Numpy (Harris et al. Reference Harris2020). A complete example of application of analytical solution to the heat equation using this variational method is presented in Appendix A.
As previously explained, we first optimise the parameters (except for
$\theta _\mu$
, which we keep fixed) using a cost function
$C = ||\psi (\theta )-\psi _{target}||^2$
(mean-squared error). We choose a number of iteration steps such that the cost function has decreased significantly and converged such that the quantum variational ansatz produces the centred target distribution (see figure 6). After this, we apply the advection step by changing the parameter
$\theta _\mu$
to match the initial mean value of the distribution on the energy grid.
The number of variational parameters depends indirectly on the value of
$\chi$
. A higher initial
$\chi _0$
corresponds to entering deeper into the quantum regime of radiation reaction, where the resulting electron distribution becomes increasingly asymmetric and skewed. This symmetry breaking typically requires a larger number of variational parameters (or circuit layers) to be accurately captured within our ansatz. Furthermore, for a given fixed simulation time, a higher
$\chi$
will lead to higher energy loss, and consequently requires a wider grid domain. From empirical parameter scans for different values of
$n$
and
$k$
, we find that, in practice, achieving good fidelity in reproducing the initial symmetric (Gaussian) distribution typically requires the number of layers
$k$
to scale approximately linearly with the number of qubits n, specifically in the range
$k \sim n$
to
$k \sim 2n$
.
Figure 6. Fitting of the initial wavefunction for
$\chi _0=10^{-3}$
,
$n=6$
qubits and
$k=5$
layers of variational parameters. (a) A centred wavefunction is obtained through fitting the variational parameters. (b) Evolution of the cost function as the optimiser converges on a good approximation of the target wavefunction.
In figure 7, we show the VarQITE simulation of the distribution function for
$\chi _0=\{ 10^{-3}, 10^{-2}\}$
, and results from solving the PDE with a classical algorithm. This PDE solver employs the same energy-grid discretisation as the VarQITE approach to provide a reference to which to compare. We use a forward Euler pusher in time and a finite-difference scheme on the energy-grid, only accounting for nearest neighbour cells. The discretised versions of the
$\partial _\gamma (\mathcal{A}f)$
advection and
$0.5\,\partial _{\gamma \gamma } (\mathcal{B}f)$
diffusion operators are then the same between the approaches.
Figure 7. Evolution of electron distribution functions for
$\chi _0=\{10^{-3},10^{-2}\}$
, for an initial spread in energy
$\sigma _0=90$
. Full line, OSIRIS simulation; dashed line, PDE solver; circles, VarQITE. In both cases,
$n=6$
qubits, while the number of layers is (a)
$k=5$
for
$\chi _0=10^{-3}$
and (b)
$k=6$
for
$\chi _0=10^{-2}$
.
In figure 8, we show the evolution of the first two moments of the distribution functions. There is a general agreement of the quantum circuit results with the OSIRIS simulations. There is a small deviation due to finite grid resolution, which can be resolved by increasing the number of qubits.
Figure 8. Evolution of the moments of the distribution functions obtained through VarQITE (circles) and compared against OSIRIS (lines). (a) First moment (mean), (b) second moment (spread). The number of layers for each case is the same as in figure 7.
In figure 9(a), we show the evolution of all the
$27$
variational parameters for the
$\chi _0=10^{-3}$
case. The trajectories are smooth and some parameters have minimal deviations, which suggests that the ansatz is over-parametrised and can be made more efficient. Figure 9(b) shows the values of the mean energy computed from the wavefunction produced by the quantum variational ansatz against the variational parameter responsible for the translation
$e^{-i \theta _\mu \partial _x }$
, where a linear scaling between the two is visible.
Figure 9. (a) Evolution of all variational parameters for
$\chi _0=10^{-3}$
, highlighting the ‘average energy’ parameter in red. (b) Linear correlation between average energy from the wavefunction and the parameter of the wavefunction translation.
3.3. A note on the extraction of moments of the distribution function
The global translation of the wavefunction is correlated with the average energy through a simple rescaling. Therefore, we would not need to measure the mean energy of the wavefunction because we could simply read this value from one of the variational parameters. This is only allowed because of our choice of a translated, even-symmetric wavefunction and would not be viable if the ansatz allowed for skewed distributions.
The average and variance of the energy can be measured efficiently (
$n\,Z_i$
and
$n^2/2$
$Z_i Z_j$
gates, respectively), with a Pauli-string decomposition of the observables:
\begin{equation} \hat {x} = \frac {1}{2}(2^n - 1) \, I - \frac {1}{2} \sum _{i=0}^{n-1} 2^{n - i - 1} Z_i, \end{equation}
\begin{equation} \hat {x}^2 = \left ( a^2 + \sum _{i=0}^{n-1} b_i^2 \right ) \, I + 2a \sum _{i=0}^{n-1} b_i Z_i + 2 \sum _{0 \leqslant i \lt j \leqslant n-1} b_i b_j \, Z_i Z_j,\\[18pt] \end{equation}
with
$a = ({1}/{2}) (2^n - 1), b_i = -({1}/{2}) \times 2^{n - i - 1}$
and
$Z_i$
is the Pauli-
$Z$
gate applied to qubit
$i$
. Based on this,
$\hat {x}=(0,1,2,3,\ldots)$
and
$\hat {x}^2=(0,1,4,9,\ldots )$
are diagonal matrices. The commutation relation between the individual terms can be used to make the measurement more efficient and the results of
$\langle x\rangle = \langle \psi |\hat {x}|\psi \rangle , \langle x^2\rangle = \langle \psi |\hat {x}^2|\psi \rangle$
can then be linearly rescaled to be compatible with the wavefunction domain in energy
$\gamma$
.
4. Conclusions
In this work, we present a quantum-hybrid algorithm approach to the quantum radiation reaction, providing similar results when benchmarked against the particle-in-cell Monte Carlo (PIC-MC) code OSIRIS. We have simulated the Fokker–Planck equation describing the stochastic cooling of an electron population in a strong magnetic field using the variational quantum imaginary time evolution (VarQITE) algorithm, showing good agreement with the results from PIC-MC. We have developed a new variational ansatz that mimics a Gaussian trial-function, capturing the first and second moments of the distribution functions while requiring fewer parameters than more general variational ansätze, by enforcing a parity symmetry on the wavefunction in addition to a direct control on its average position through the advection operator. We have also derived closed-form analytical results for moments of the energy distribution functions, their entropy and autocorrelation functions.
Future steps will include extending this approach to other Fokker–Planck equations of interest (for example, the Kompaneets equation or laser cooling of trapped atoms), to the full Boltzmann equation (non-local operator) and to represent the quantum system not as a classical distribution function but as a Fock state (Hidalgo & Draper Reference Hidalgo and Draper2024), where the dynamics would be naturally unitary. Additionally, using specialised circuits for loading specific distribution functions as part of the variational ansatz (Gaussian, uniform, Chi-squared) could lead to a reduced number and more interpretable parameters, which would speed up the VarQITE algorithm and possibly mitigate the barren plateau problem.
Our work can contribute to further development of both classical diffusive plasma physics and to simulate the fully quantum nature of plasma interactions, namely the transition from semi-classical to non-perturbative SFQED regimes.
Acknowledgments
The authors thank the anonymous referees for their insightful comments and suggestions, which have significantly contributed to enhancing the clarity and rigor of this work. The authors also thank the following people for fruitful discussions / proofreading parts of the manuscript: Mr José Mariano on Monte Carlo sampling, Mr Lucas Ansia on classical simulation of the Fokker–Planck equation, Mr Bernardo Barbosa on the quantum radiation reaction in SFQED, Mr Pablo Bilbao on plasma dynamics in strong magnetic fields, Mr Gabriel Almeida on expressibility and complexity of variational quantum circuits, Mr Anthony Gandon on efficient Pauli-String decomposition of Hamiltonians for variational quantum circuits, Mr Efstratios Koukoutsis on Lindbladian/dissipative quantum simulations and Mr Diogo Cruz on quantum algorithms for PDE solvers. Simulations were performed at the IST cluster (Lisbon, Portugal). We have used the quantum frameworks Qiskit (Javadi-Abhari et al. Reference Javadi-Abhari2024) and PennyLane (Bergholm et al. Reference Bergholm2022).
Editor Nuno Loureiro thanks the referees for their advice in evaluating this article.
Declaration of interest
The authors report no conflict of interest.
Funding
This work was supported by the Portuguese Science Foundation (FCT) Grants Nos. PTDC/FIS-PLA/3800/2021, UI/BD/153735/2022. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 – EUROfusion). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
Data availability
The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/OsAmaro/QuantumFokkerPlanck. The repository includes notebooks explaining the classical PDE solver, the VarQITE approach, the analytical models and input decks for OSIRIS simulations.
Appendix A. Applying the variational principle to the heat equation
In this appendix, we provide a step-by-step application of McLachlan’s variational principle (MVP) to the simplest diffusive PDE, the heat equation. Similar to the classical variational principle used in quantum mechanics to obtain approximate ground states of a system, the MVP for imaginary or real-time evolution is useful both numerically and analytically.
The one-dimensional (1-D) homogeneous heat equation
\begin{equation} \partial _t u = k \,\partial _{xx} u \end{equation}
with thermal diffusivity coefficient
$k$
, admits a kernel solution (when
$u(t=0,x)=\delta (x)$
)
\begin{equation} u(t,x) = \dfrac {1}{\sqrt {4 \pi k t}} \exp \left ( -\dfrac {x^2}{4 k t} \right ), \end{equation}
from which analytical solutions to more general initial conditions can be constructed.
This suggests using an ansatz
\begin{equation} v(A,\sigma ) = A \exp \left ( -\dfrac {x^2}{2\sigma ^2} \right ), \quad \partial _t v = \hat {L}_H\,v, \end{equation}
with
$\hat {L}_H$
the linear operator generating the evolution of the system (in this case, it is the Laplacian
$\partial _{xx}$
) and variational parameters
$\theta _i(t)=(A(t),\sigma (t))$
for the analytical solution of the MVP. We take the thermal diffusivity to be
$k=1$
for simplicity. Since this function is
$L_1$
-normalised (that is,
$\int v \,\mathrm{d}x = 1$
), the normalisation
$A$
will be correlated and a function of the parameter
$\sigma$
.
The following derivatives and integrals are used in the calculations:
\begin{equation*} \partial _A v = A^{-1}\,v, \quad \partial _\sigma v = x^2\,\sigma ^{-3} \,v, \quad \partial _{xx} v = v\,(x^2-\sigma ^2)\,\sigma ^{-4}, \end{equation*}
\begin{equation} \int v \,\mathrm{d}x = A \sqrt {2\pi } \sigma , \int v^2 \,\mathrm{d}x = A^2 \sqrt {\pi } \sigma, \end{equation}
\begin{equation*} \int x^2 v^2 \,\mathrm{d}x = A^2 \sqrt {\pi } \sigma ^3 /2, \int x^4 v^2 \,\mathrm{d}x = 3 A^2 \sqrt {\pi } \sigma ^5 /4. \end{equation*}
From the MVP (3.3), we need to compute the matrix elements
$\int (\partial _k v )(\partial _j v )\,\mathrm{d}x$
, which require either the
$\int v^2\,\mathrm{d}x$
or
$\int x^2\,v^2\,\mathrm{d}x$
integrals. We obtain
\begin{equation*} M_{A,A} = \sqrt {\pi } \sigma , \,M_{A,\sigma } = M_{\sigma ,A} = \sqrt {\pi } A/2, \,M_{\sigma ,\sigma } = 3\sqrt {\pi } A^2 /(4\sigma ). \end{equation*}
For the column vector, we need to compute
$\int (\partial _k v) (\hat {L}_H v)\,\mathrm{d}x$
. We obtain
\begin{equation*} V_A = -A \sqrt {\pi } / (2\sigma ), \quad V_\sigma = A^2 \sqrt {\pi } / (4\sigma ^2). \end{equation*}
We now need to solve the matrix equation
$M_{k,j} \dot {\theta }_j = V_k$
. Factors of
$\sqrt {\pi }$
will cancel out:
\begin{equation} \begin{bmatrix} \sigma &\quad A/2 \\ A/2 &\quad 3A^2/(4\sigma ) \end{bmatrix} \begin{bmatrix} \dot {A} \\ \dot {\sigma } \end{bmatrix} = V_k \leftrightarrow \begin{bmatrix} \dot {A} \\ \dot {\sigma } \end{bmatrix} = \begin{bmatrix} 3/(2\sigma ) &\quad -1/A \\ -1/A & \quad 2\sigma /A^2 \end{bmatrix} \begin{bmatrix} -A/(2\sigma ) \\ A^2 /(4\sigma ^2) \end{bmatrix} = \begin{bmatrix} -A/\sigma ^2 \\ 1/\sigma \end{bmatrix},{}_{\vphantom{\Big(}} \end{equation}
where we explicitly computed the matrix inverse of
$M_{k,j}$
. The second equation (for
$\sigma$
) is separable. We can then replace
$\sigma (t)$
in the first equation to obtain the amplitude
$A(t)$
. The parameters thus become
\begin{equation} \sigma (t) = \sqrt {2t + \sigma _0^2}, \,A(t) = A_0 \sigma _0/\sqrt {2 t + \sigma _0^2} ,\end{equation}
where
$\sigma _0, A_0$
are the initial spread and amplitude, respectively. This is precisely the exact analytical solution (A.2). In figure 10, we show a comparison between a standard PDE solver solution of the heat equation, a numerical solution of the MVP equations using the ansatz (A.3) (without a quantum circuit), and the analytical solution (A.6). The numerical grid had
$2^6$
cells, the number of time steps was
$200$
for a maximum simulation time of
$40$
and a small initial spread of the distribution function
$\sigma _0=0.01$
.
Figure 10. Numerical solution of the 1-D heat equation. (a) Snapshots of the distribution
$u(t,x)$
from a PDE solver and numerical evolution of a Gaussian ansatz in the MVP. (b) Evolution of the variational parameters: colour, numerical MVP; dashed, analytical solution to the MVP (A.6).
Appendix B. The Fokker–Planck equation
B.1. Review of theoretical study of radiation reaction
In this section, we review past works on the radiation reaction, introduce the main parameters that describe the interaction of electrons with an intense magnetic field and derive analytical solutions for some observables of interest. In the past decades, the theoretical modelling of electron radiation losses in intense fields has been extensively studied. Shen & White (Reference Shen and White1972) derive the transport partial differential equation (PDE) for electrons in a constant uniform magnetic field. Analytical results of multi-photon scattering of electrons in periodic structures (e.g. oriented-crystals) have been derived by Khokonov (Reference Khokonov2004) and Bondarenco (Reference Bondarenco2014); however, without simple explicit closed-form expressions for the observables, nor immediate application in electron-laser scattering. Lobet et al. (Reference Lobet, d’Humières, Grech, Ruyer, Davoine and Gremillet2016) numerically studied the coupled dynamics of photons and leptons, including radiation reaction and pair production in a constant magnetic field. Vranic et al. (Reference Vranic, Grismayer, Fonseca and Silva2016), Ridgers et al. (Reference Ridgers2017) and Niel et al. (Reference Niel, Riconda, Amiranoff, Duclous and Grech2018) derive the ordinary differential equations (ODEs) for the first two momenta of the energy distributions of electrons interacting with a laser pulse, but do not present their corresponding closed-form explicit solutions. Niel et al. (Reference Niel, Riconda, Amiranoff, Duclous and Grech2018) show connections between the different regimes of radiation reaction and the corresponding particle trajectories and stochastic PDEs. Artemenko et al. (Reference Artemenko, Krygin, Serebryakov, Nerush and Kostyukov2019) numerically solved the electron transport PDE for different magnetic field values using the ‘global constant field approximation’ to map between a constant magnetic field set-up and a laser pulse with a finite duration envelope. Jirka et al. (Reference Jirka, Sasorov, Bulanov, Korn, Rus and Bulanov2021) derive the average energy decay in a pulsed plane-wave laser in the highly ‘nonlinear quantum’ high-
$\chi$
regime (the quantum nonlinear parameter,
$\chi$
, is defined in §2).
More recently, some papers have contributed with new analytical results for the study of classical (CRR) and quantum (QRR) radiation reaction regimes. Bilbao & Silva (Reference Bilbao and Silva2023) derive analytically the evolution of an electron momentum ‘ring distribution’ in a maser set-up, relevant for astrophysical plasma scenarios. The formulae are derived in the CRR regime using the ‘method of lines’ technique but appear to hold even for
$\chi \sim 0.1$
, and suggest that this class of ring distributions is a ‘global attractor’. Zhang, Zhang & Zhou (Reference Zhang, Zhang and Zhou2023) derive approximate analytical formulae for the evolution of a lepton distribution function in a constant magnetic field during multiphoton scattering in the QRR
$\chi \lt 1$
regime. The authors identify and analytically characterise the ‘quantum peak splitting’, which occurs when an initially peaked distribution function evolves into a doubly peaked distribution. Bulanov et al. (Reference Bulanov, Grittani, Shaisultanov, Esirkepov, Ridgers, Bulanov, Russell and Thomas2024) derive CRR closed-form expressions for distribution functions for a set of initial conditions of interest, again using the ‘method of lines’, and find an equilibrium (time-independent) solution for the QRR Fokker–Planck equation of an electron beam in an LWFA set-up. Kostyukov et al. (Reference Kostyukov, Nerush, Mironov and Fedotov2023) derive short-time evolution expressions for the electron wave packet in a strong EM field, including the expectation value of the electron spin. The approach relies on Volkov functions and the Dyson–Schwinger equation, and naturally leads to the damping (radiation reaction) of the wave packet. Torgrimsson (Reference Torgrimsson2024a
,Reference Torgrimsson
b
) derives explicit analytical time-dependent electron distribution functions in momentum and spin through resummation techniques. Following a similar approach, Blackburn (Reference Blackburn2024) derives explicit analytical QRR formulae for the final electron mean energy and spread as functions of initial electron energy, laser
$a_0$
and pulse duration. These expressions are derived in the low-
$\chi$
regime, but remain approximately valid for a larger range of parameters.
B.2. Regimes of radiation reaction
Here, we follow the notation of Vranic et al. (Reference Vranic, Grismayer, Fonseca and Silva2016). When
$\chi \ll 1$
, the equation for energy loss can be derived from the Landau–Lifshitz equation, leading to
\begin{equation} \dfrac {\mathrm{d}\gamma }{\mathrm{d}t} = -c_{rr} \,\gamma ^2, \end{equation}
with
$c_{rr} \equiv 2 \,\omega _c e^2 b_0^2/(3 m c^3)$
a constant defining the effect of classical radiation reaction. In the set-up considered, three physical parameters completely describe the dynamics: the initial electron energy
$\gamma _0$
, the external magnetic field
$B$
and the evolution time
$t$
. The solution for the energy is
\begin{equation} \gamma _c(t)\sim \gamma _0/(1+c_{rr}\,t\,\gamma _0). \end{equation}
This scaling for the average energy can be corrected in the
$\chi \lesssim 1$
regime.
From a distribution point of view, the relativistic Vlasov equation for the electron distribution function
$f$
is modified. The radiation reaction no longer preserves phase-space volume, corresponding to the operator on the right-hand side of the equation
\begin{equation} \frac {\mathrm{d} f(t, \vec {x}, \vec {p})}{\mathrm{d} t} = \frac {\partial f}{\partial t} + v \cdot \nabla _x \,f + F \cdot \nabla _p f = \mathcal{C}[f_e] ,\end{equation}
where
$\mathcal{C}$
is a collision operator whose form depends on the model/regime of radiation reaction considering Niel et al. (Reference Niel, Riconda, Amiranoff, Duclous and Grech2018). In our simplified set-up, the coordinate gradient and the force terms vanish, so the only term of interest is the collision operator.
For higher values of
$\chi$
, the emission of photons becomes stochastic, in general, leading to a spread in energy. Conceptually, the electrons can be thought of as being slightly scattered through a cloud of photons of the field, where each photon has an energy much lower than that of the electron.
There are several approximation methods for the theoretical analysis and numerical sampling of the cross-sections contained in the collision operator
$\mathcal{C}$
. One of the most commonly used is the locally constant field approximation (LCFA), which is exact for a plane wave laser field, and where the electric, magnetic and Poynting vectors are mutually orthogonal.
Within this approximation, the probability of emission of a photon with
$\chi _\gamma$
by an electron with
$\chi$
per unit time, also known as the nonlinear Compton scattering (nCS) differential cross-section/rate, is
\begin{align} \frac {\mathrm{d}^{2} N_\gamma }{\mathrm{d} t \ \mathrm{d} \chi _\gamma } (\chi ,\chi _\gamma ) =\frac {\alpha m c^{2}}{\sqrt {3} \pi \hbar \gamma \chi }\left [ \left (1-\xi +\dfrac {1}{1-\xi } \right ) K_{2 / 3}(\tilde {\chi }) - \int _{\tilde {\chi }}^{\infty } \ K_{1 / 3}(x) \,\mathrm{d} x \right ] \end{align}
with
$\xi \equiv \chi /\chi _\gamma$
,
$\tilde {\chi } \equiv 2\xi /(3\chi (1-\xi ))$
and
$K_n$
is the modified Bessel function of second kind (Abramowitz & Stegun Reference Abramowitz and Stegun1964). This rate in
$\chi _\gamma$
can be mapped to a rate in particle energy
$\mathrm{d}^2 N_\gamma /\mathrm{d}t\, \mathrm{d}\gamma _\gamma \,(\gamma ,\gamma _\gamma )$
.
In the LCFA approximation, the Lorentz invariant is
$E^2-B^2=0$
. In the case of a constant magnetic field, the boosted electromagnetic field in the electron rest frame is almost a perfect plane wave, with a very slight deviation from this approximation. The result is that the nCS and the quantum synchrotron emission spectra for ultra-relativistic leptons are in practice very similar, as can be confirmed by Niel et al. (Reference Niel, Riconda, Amiranoff, Duclous and Grech2018) and Zhang et al. (Reference Zhang, Zhang and Zhou2023).
The Fokker–Planck (FP) equation can be seen as an extension of the standard Vlasov equation to kinetically and stochastically model collisions between plasma species and laser–plasma interaction. An important application is simulating the energy loss of an electron beam as it interacts with electromagnetic fields (Neitz & Di Piazza Reference Neitz and Di Piazza2013; Vranic et al. Reference Vranic, Grismayer, Fonseca and Silva2016). In this case, the particle distribution evolves through
\begin{equation} \frac {\partial f(t, \vec {p})}{\partial t}=\frac {\partial }{\partial p_{l}}\left [-\mathcal{A}_{l} f+\frac {1}{2} \frac {\partial }{\partial p_{k}}\left (\mathcal{B}_{l k} f\right )\right ] \end{equation}
with
$\mathcal{A}_{l} \equiv \int q_{l} w(\vec {p}, \vec {q}) \mathrm{d}^{3} \vec {q}$
,
$\mathcal{B}_{l k} \equiv \int q_{l} q_{k} w(\vec {p}, \vec {q}) \mathrm{d}^{3} \vec {q}$
, the drift and diffusion coefficients, respectively, and
$ w(\vec {p}, \vec {q})\,\mathrm{d}^3\vec {p}$
is the probability per unit time of momentum change of the electron
$\vec {p}\rightarrow \vec {p}-\vec {q}$
, with
$\vec {q}$
the momentum of the photon, and where Einstein index contraction was assumed.
In the co-linear approximation of radiation emission,
$\chi _\gamma /\chi \sim \hbar k/(\gamma m c)$
, this problem becomes essentially 1-D, with drift and diffusion coefficients
\begin{equation} \mathcal{A} \equiv \dfrac {\gamma m c}{\chi } \int _0^{\chi } \chi _\gamma \,\dfrac {\mathrm{d}^2 N_\gamma }{\mathrm{d}t\,\mathrm{d}\chi _\gamma } \, \mathrm{d}\chi _\gamma , \quad \mathcal{B} \equiv \dfrac {(\gamma m c)^2}{\chi ^2} \int _0^{\chi } \chi _\gamma ^2 \,\dfrac {\mathrm{d}^2 N_\gamma }{\mathrm{d}t\,\mathrm{d}\chi _\gamma } \, \mathrm{d}\chi _\gamma. \end{equation}
In the regime of
$\chi \ll 1$
, the
$\mathcal{A}$
and
$\mathcal{B}$
coefficients have polynomial approximations (Vranic et al. Reference Vranic, Grismayer, Fonseca and Silva2016)
\begin{equation} \begin{split} \mathcal{A} \sim \dfrac {2}{3} \dfrac {\alpha m^2 c^3}{\hbar } \chi ^2 \propto \gamma ^2, \quad \mathcal{B} \sim \dfrac {55}{24 \sqrt 3} \dfrac {\alpha m^3 c^4}{\hbar } \gamma \, \chi ^3 \propto \gamma ^4. \end{split} \end{equation}
The Fokker–Planck equation is sparse (only two local, differential operators), while the Boltzmann equation is non-local and its corresponding Hamiltonian evolution matrix is dense. Consequently, both numerical (either a PDE solver or Monte Carlo sampling) and analytical solutions of the latter equation are often more challenging than for the former.
B.3. Entropy and autocorrelation
Having derived the evolution of the spread (2.6) and assuming a Gaussian functional form of the distribution function, the Shannon entropy of the electron beam evolves as
$S(t) = -\int f \log (f) \,\mathrm{d}\gamma \sim \log (\sigma (t)) + c^{te}$
up to an additive constant (which we choose such that the simulation and theory curves match at late times). Since this is a monotone function of the spread, it has the same qualitative behaviour and peaks at
$\log (\sigma _{\max })$
. From a physics point-of-view, the entropy changes due to a competition between the drift and diffusion processes. At late times, the entropy of the electrons decreases, which can be interpreted as a transfer of entropy to the radiation that is being emitted.
One can also compute an approximate auto-correlation function
$g(\tau )$
(not to be confused with the gaunt factor) from the classical solution
$\gamma _c(t) \sim \gamma _0/(1+2R_c t/3)$
as
\begin{equation} \begin{split} g(\tau ) \equiv \dfrac {\langle \gamma (t) \gamma (t+\tau ) \rangle }{\langle \gamma ^2(t) \rangle } = \int _0^\infty \dfrac {1}{(1+2 R_c t/3)(1+2 R_c (t+\tau )/3)}\, \mathrm{d}t \\ \propto \dfrac {\log (1+2 R_c \tau /3)}{2 R_c \tau /3}, \end{split} \end{equation}
where normalisation
$g(0) = 1$
is enforced. This quantifies how correlated two electron trajectories are on average if measured with a time delay of
$\tau$
. Using instead the second-order expansion for the average energy (that is, the full (2.5)) leads to a more accurate closed-form solution, albeit with more terms (not shown here). Considering the stochastic component of the trajectories would require computing expectation values of nonlinear functions of Wiener processes
$W(t), W(t+\tau )$
. The solution to the Fokker–Planck equation
$f(t,\gamma )$
does not give us information on the auto-correlation directly, and a marginal distribution
$f(t_1,\gamma _1; t_2, \gamma _2)$
would have to be obtained.
In figure 11, we show the comparison between these two quantities
$S(t)$
,
$g(\tau )$
from OSIRIS simulations computed from distribution histograms and individual particle trajectories, and the approximate analytical results. Again, for
$\chi _0=10^{-1}$
(blue curve), the
$t$
-axis is again multiplied by a factor of
$2$
to have the same range as the other curves.
B.4. McLachlan’s variational principle approach
Here, we apply the technique described in Appendix A to the Fokker–Planck equation. Since the number of particles (
$L_1$
norm) is conserved and we are mostly interested in the first moments of the distribution function
$\mu =\int \gamma f\,\mathrm{d}\gamma$
,
$\sigma ^2 = \int (\gamma -\mu )^2 f\,\mathrm{d}\gamma$
, we enforce this into the ansatz
\begin{equation} v(\mu ,\sigma ) = \dfrac {1}{\sqrt {2\pi \sigma ^2}} \exp \left ( -\dfrac {(\gamma -\mu )^2}{2\sigma ^2} \right ),\,\,\partial _t v = \hat {L}_{FP}\,v \!=\! -\partial _\gamma (a \,\gamma ^2 \,v) + 0.5 \,\partial _{\gamma \gamma } (b\, \gamma ^4\, v) ,\end{equation}
with
$\hat {L}_{FP}$
the linear operator generating the dynamics, and variational parameters
$\theta _i(t)=(\mu (t),\sigma (t))$
, for the analytical solution of the MVP. We could have used an ansatz with the amplitude
$A(t)$
, which would have led to a
$3\times3$
matrix system of equations.
The following derivatives will be used in the following calculations:
\begin{equation} \dfrac {\partial _\mu v}{v} = \dfrac {\gamma -\mu }{\sigma ^2}, \quad \dfrac {\partial _\sigma v}{v} = \dfrac {(\gamma -\mu )^2-\sigma ^2}{\sigma ^3}. \end{equation}
From the MVP (3.3), we need to compute the matrix elements
$\int (\partial _k v )(\partial _j v ) \mathrm{d}\gamma$
, which require either the
$\int (\gamma -\mu )^2 \,v^2 \,\mathrm{d}\gamma$
or the
$\int (\gamma -\mu )^4\,v^2 \,\mathrm{d}\gamma$
Gaussian integrals. In general, the change of variables
$(\gamma -\mu )/\sigma \rightarrow \gamma$
simplifies the calculations. We obtain
\begin{equation*} M_{\mu \mu } = 1/(4 \sqrt {\pi } \sigma ^3), \quad M_{\sigma \sigma } = 3 / (8 \sqrt {\pi } \sigma ^3 ), \quad M_{\mu \sigma } = M_{\sigma \mu } = 0, \end{equation*}
where the cross-term is null due to the odd-symmetry of the integrand function. The diagonal structure of this matrix simplifies the calculations considerably. The inverse matrix becomes
\begin{equation} M_{kj}^{-1} = 4 \sqrt {\pi } \sigma ^3 \begin{bmatrix} 1 & \quad 0 \\ 0 & \quad 2/3 \end{bmatrix}. \end{equation}
For the column vector, we obtain
\begin{equation*} V_\mu = \dfrac {a(-2\mu ^2+\sigma ^2)}{8 \sqrt {\pi } \sigma ^3} -\dfrac {3 b \mu (2 \mu ^2 +\sigma ^2)}{8 \sqrt {\pi } \sigma ^3}, \end{equation*}
\begin{equation*} V_\sigma = -\dfrac {3 a \mu }{4 \sqrt {\pi } \sigma ^2} - \dfrac {3b(-4 \mu ^4 +36 \mu ^2 \sigma ^2 +13 \sigma ^4 )}{64 \sqrt {\pi } \sigma ^3}. \end{equation*}
Applying the inverse matrix
$M_{kj}^{-1}$
to the column vector
$V_k$
, we obtain
\begin{equation} \begin{bmatrix} \dot {\mu } \\ \dot {\sigma } \end{bmatrix} = \begin{bmatrix} \underline {-a \mu ^2} + a \sigma ^2/2 - 3 b \mu (2\mu ^2 + \sigma ^2)/2 \\ \,\underline {-2 a \mu \sigma + b \mu ^4/(2\sigma ) } - b ( 36\mu ^2 \sigma ^2 +13 \sigma ^4)/(8 \sigma ) \end{bmatrix}. \end{equation}
Changing notation
$2 \alpha _{rr} = a$
and
$b \mu ^4 = B/(m^2 c^2)$
, the underlined terms, which are leading order assuming
$\mu \gg \sigma$
, recover the results of (8) and (14) from Vranic et al. (Reference Vranic, Grismayer, Fonseca and Silva2016).
In figure 12, we show a comparison between a PDE solver solution of the Fokker–Planck equation, a numerical solution of the MVP-VarQITE equations using the ansatz (B.9) (without a quantum circuit), and the analytical solution (2.5). The numerical grid had
$2^7$
cells, the number of time steps was
$900$
for a maximum simulation time of
$20$
,
$\chi _0=10^{-2}$
and initial momenta
$\mu _0=1800, \sigma _0=20$
. We have thus shown the potential of MVP for retrieving the distribution moments in the quantum radiation reaction.
Figure 12. Numerical simulation of the Fokker–Planck equation with
$\chi _0=10^{-2}$
. (a) Evolution of the distribution function
$f(t,\gamma )$
. The dashed line represents the result from a standard PDE solver, while circles represent the numerical integration of the MVP equations. (b) Evolution of the moments and variational parameters
$(\mu , \sigma )$
. The coloured lines represent the moments of the distribution function obtained from the PDE solver, while dashed lines represent the analytical results of (2.5) and (2.6).
Open access
