LICENTIATE THESIS
Modeling concepts for Laser Wakefield Acceleration
FRIDA BROGREN
Department of Physics
University of Gothenburg
Gothenburg, Sweden 2025
Modeling concepts for Laser Wakefield Acceleration
Frida Brogren
This thesis is electronically published, available at
https://hdl.handle.net/2077/90143
Department of Physics
University of Gothenburg
SE-412 96 Gothenburg
Sweden
Telephone: +46 (0)31-786 00 00
Printed by Stema Specialtryck AB
Borås, Sweden 2025
ABSTRACT
Laser wakefield accelerators have strong potential as future compact sources
of high-charge relativistic electrons. However, their practical application is
currently limited by significant shot-to-shot variations, large energy spreads,
and high emittance. These fluctuations primarily arise from the high sensi-
tivity of the acceleration process on the properties of the driving laser. To
increase the knowledge of this interaction, I consider two complementary
methods to model laser wakefield acceleration: data-driven modeling and
Particle-In-Cell (PIC) simulations.
The data-driven method employs a neural network to model relationships
between laser properties, electron properties and secondary X-ray generation
in a laser wakefield accelerator. In experimental setups, measurement errors
are prevalent and in particular the developed model mitigates the influence
of these by employing a novel optimization criteria. Analysis of the model
reveals strong consistency with physics-based simulations and theoretical
expectations and, additionally, reveals unknown correlations between laser
driver and system output. Thus, the model shows strong potential as a
foundation for the stabilization and analysis of laser wakefield accelerators
in experimental settings.
On the other hand, for the design and study of laser-plasma interactions,
physics-based models, such as particle-in-cell (PIC) simulations are essential.
These simulations come at a significant computational cost, which can be
levitated by performing them in a Lorentz-boosted frame. However, this
approach is often limited by the Numerical Cherenkov Instability (NCI). In
this thesis, a lightweight PIC framework, π-PIC, is introduced to facilitate
easy testing and integration of new PIC algorithmic developments. Using
this framework, the influence of energy and momentum conservation on
NCI growth is examined.
Keywords: Laser-plasma accelerator, Neural network, Particle-In-Cell, Laser
Wakefield acceleration, Numerical Cherenkov Instability

LIST OF PAPERS
This thesis serves as an introductory text to the following appended papers:
Paper [I]
F. BROGREN, S. JALAS, L. HÜBNER, P. MESSNER, M. SCHNEPP, M. TRUNK, C.
WERLE, P. WINKLER, M. MARKLUND, A. GONOSKOV, W. P. LEEMANS, A. R. MAIER,
M. KIRCHEN Data-driven modeling of a laser-plasma accelerator-based X-ray
source Published in Physical Review Accelerators and Beams.
Contribution: I developed the models, performed the analysis and PIC sim-
ulations as well as provided the initial draft of the article, which was finalized
together with M. Kirchen and S. Jalas. The collection of data as well as the
design of the laser-plasma accelerators was carried out by the LUX team.
Paper [II]
F. BROGREN, J. MAGNUSSON, C. OLOFSSON, A. GONOSKOV π-PIC: a framework
for modular particle-in-cell developments and simulations. Preprint.
Contribution: I developed the code for the Moving Window and Absorbing
Boundaries extensions, as well as for the Electrostatic solver. I conducted the
code benchmarking and generated all figures in the article, except Figure 3,
which was prepared by C. Olofsson. I wrote the initial draft of the manuscript,
which was subsequently revised by the other authors.
To my Grandmother, who has stayed curious, happy, and kind throughout
the years. It is truly inspiring. I hope you like it!
ACKNOWLEDGEMENTS
I would like to thank my parents, who through my upbringing installed the
confidence in me that I can do it if I want to; my supervisor and colleagues,
who show me ways to make it possible; my partner, who supports me when
it feels too much; and my sisters, who remind me that there is more to life. It
would not be possible without you. Thank you!

CONTENTS
1 Introduction 1
2 Physics of laser wakefield acceleration 4
2.1 Plasma dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.1 Conservational properties of plasma . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.2 Plasma waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Laser wakefield acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Ponderomotive force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 Wake excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.3 Physical limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.4 Electron trapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 Data-driven modeling of laser wakefield acceleration 18
3.1 Prediction and estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Orthogonal distance loss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 Experimental setup and data selection . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4 Electron and X-ray dependence on laser properties . . . . . . . . . . . . . . . . 23
3.5 Prediction of X-ray spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.6 Conclusion and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4 Particle-in-cell modeling of laser wakefield acceleration 28
4.1 The PIC algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.1.1 An electrostatic PIC algorithm: a detailed example . . . . . . . . . . . . . . . 33
4.2 Numerical stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.2.1 The Courant-Friedrichs-Lewy condition . . . . . . . . . . . . . . . . . . . . . . 37
4.2.2 Energy and momentum conservation . . . . . . . . . . . . . . . . . . . . . . . . 38
4.2.3 Charge conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2.4 Numerical dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2.5 Numerical Cherenkov instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.3 π-PIC: a framework for modular PIC simulations . . . . . . . . . . . . . . . . . 42
4.3.1 Architecture of π-PIC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.3.2 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.4 Study of numerical Cherenkov instability . . . . . . . . . . . . . . . . . . . . . . . 45
4.4.1 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.4.2 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.5 Conclusion and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5 Research papers 50
References 91
A Derivation of numerical dispersion relation 99
1
1 Introduction
99.9% of the observable universe is believed to consist of plasma, the state of
matter where electrons have enough energy to escape the attraction of the
nuclei, and are moving freely. Apart from being abundant in the universe,
plasma features in many interesting phenomena, such as fusion, star and
galaxy formation, solar physics and lightning. Plasmas can create extremely
strong electromagnetic fields[1], violent instabilities and shocks[2], while it
also admits stable states consisting of periodic waves. Astonishingly, much of
this physics is governed by the same set of self-consistent equations, which,
despite their simplicity, can be extremely difficult, or even impossible, to
solve analytically.
One plasma phenomenon of interest to the particle accelerator community
is the interaction between a laser pulse and plasma. In this interaction the
pulse excites a plasma oscillation which can be use to accelerate electrons to
high energies. This is the set up of the laser wakefield accelerator[3]. These
sources can deliver high-energy, high-charge electron beams at centimeter
acceleration lengths[1]. As such, the technique enables generation of rela-
tivistic electron beams suitable for applications in medicine, bio-chemistry
and fundamental physics at a fraction of the cost of conventional accelerator
facilities.
While promising, laser wakefield-based electron sources are, for many prac-
tical applications, limited by high shot-to-shot variability as well as large
energy spread and emittance [4]. The main reason for this is that the excita-
tion of a plasma wake, its amplitude, and the electron injection are heavily
affected by the wavefront and temporal distribution of the driving laser. This
interaction involves a complex non-linear relationship, which generally lacks
a detailed analytical description. The subject of this thesis is alternative
ways of modeling laser wakefield acceleration (LWFA) with the motivation to
increase detailed knowledge and improving stability.
To model the laser-plasma interaction, this work presents two complemen-
tary approaches: (i) data-driven modeling and (ii) physics-based particle-in-
cell (PIC) simulations.
In the first approach, presented in chapter 3, a data-driven method is devel-
oped to model the relationship between properties of laser, electrons and
X-ray radiation in an experimental study of a laser wakefield accelerator
2 INTRODUCTION
with subsequent X-ray generation using an undulator. The study focuses on
extracting the underlying dependence between input and output properties
of the system. To achieve this a neural network model is designed and op-
timized based on minimizing the orthogonal distance between model and
data points in order to avoid bias, known as attenuation bias, originating
from errors in input variables. Using the model the dependence between
input and output variables of the system is extracted and analyzed in an
experimental setting.
The second approach, presented in chapter 4, concerns modeling LWFA
using PIC codes. PIC codes are an invaluable tool for the studies of plasma in-
teractions. The usefulness of this method is partly due to the code’s capacity
to maintain consistency while incorporating physical processes not covered
by the PIC plasma description, such as the ionization of atoms, radiation
reactions, and collisional effects. Moreover, several modifications to both
particle and field updates can be implemented to maintain numerical stabil-
ity across various simulation conditions. In this thesis, a framework for PIC
simulations, called π-PIC, is presented which is designed as a lightweight
PIC code for extending simulations to multi-physics processes as well as
implementing and testing new algorithms.
Although there has been large progress over the past decades in computing
power and algorithmic efficiency, PIC codes (still) comes at a large computa-
tional cost, often rendering simulations of long timescales or large spaces
infeasible. In a Lorentz-boosted frame the laser wavelength is elongated
and plasma contracted – in theory admitting longer time and space steps
compared to the lab frame [5, 6]. However, practically the implementation is
hindered by a disruptive numerical instability, called Numerical Cherenkov
Instability (NCI) [7], which originates from the interaction of relativistically
moving particles with the discretized fields. In this study, π-PIC, equipped
with PIC algorithms that preserve energy exact and momentum approxi-
mately, is employed to investigate whether preserving these conservation
laws mitigates the numerical Cherenkov instability.
This thesis is structured as follows. In chapter 2 the general theory of LWFA
is presented. The chapters 3 and 4 present the two modeling approaches,
each beginning with an introduction and topic-specific theory not covered
in chapter 2. The chapter on data-driven modeling introduces the distinc-
tion between prediction and estimation (section 3.1), while the PIC chapter
describes the algorithm (section 4.1) and numerical stability (section 4.2).
3
Subsequently, both chapters detail the methods and tools used: modified
optimization criteria (section 3.2) and experimental setup (section 3.3) for
the data-driven approach, and the π-PIC framework (section 4.3) and sim-
ulation setup for studying numerical Cherenkov instability (section 4.4.1)
for the PIC simulations. Results are presented in sections 3.4 and 3.5 for
the data-driven modeling approach and section 4.4.2 for the PIC modeling
approach. Each chapter ends with a conclusion and outlook.
4 PHYSICS OF LASER WAKEFIELD ACCELERATION
2 Physics of laser wakefield acceleration
As stated earlier, plasma is the state of matter in which some, or all, electrons
have enough energy to escape the pull of the nuclei and are moving freely.
To simplify the derivation of the evolution of this state, some assumptions
are commonly made: (i) there is no recombination of electrons and atoms;
(ii) the plasma is quasi-neutral, meaning that the total charge of ions equals
the charge of electrons; (iii) the ions are immobile, serving as a stationary
neutralizing background. The last assumption is justified by the difference
in mass between ions and electrons, effectively making ions immobile on
relevant timescales.
The governing equations of plasma dynamics can be derived using the distri-
bution function of particles in six dimensional phase space: f (t , x⃗ , p⃗ ), which
in the kinetic picture can be expressed as a sum of delta functions
Np
∑
f (t , x⃗ , p⃗ ) = δ3(x⃗i (t )− x⃗ )δ3(p⃗i (t )− p⃗ ), (2.1)
i=0
where Np is the number of particles in the system.
Particle density (n) and density velocity (u⃗) are defined as moments of the
distribution function
∫
n (t , x⃗ ) = f d 3p⃗ , (2.2)
∫
p⃗
u⃗ (t , x⃗ ) = f d 3p⃗ , (2.3)
m
where m is the mass of the electron.
Incompressibility of phase space (Louisville’s theorem) states that the num-
ber of particles in a volume element following particle trajectories is con-
served in time. Mathematically, this is expressed as [8, Ch. 5.2, pp. 140-144]
 ‹  ‹
∂ f ∂ x⃗ ∂ p⃗
=−∇
∂ t x
· f −∇p · f . (2.4)∂ t ∂ t
Setting ∂ x⃗∂ t = p⃗/m and
∂ p⃗
∂ t = F⃗ where F⃗ is the force excreted on particles, we
can evaluate the right hand side
∂ f p⃗
+ ·∇
∂ t m x
f + F⃗ ·∇p f =− f ∇p · F⃗ . (2.5)
5
Note that∇x ·p⃗ = 0 since p⃗ is a coordinate in this space. The term on the right-
hand side accounts for contributions from momentum-dependent forces,
such as collisional effects and those arising from quantum electrodynamics
(QED).
By assuming a collisionless plasma and neglecting QED effects, that is setting
f ∇p · F⃗ = 0, we obtain the Vlasov equation [9],
∂ f p⃗
+ ·∇x f + F⃗ ·∇p f = 0. (2.6)∂ t m
With the insertion F⃗ = q (E⃗ + u⃗ × B⃗/c ), the equation describes a collisionless
electromagnetic plasma, where additionally, the evolution of the electromag-
netic field is governed by Maxwells equations
∂ E⃗
= c∇× B⃗ −4π j⃗ , (2.7)
∂ t
∂ B⃗
=−c∇× E⃗ , (2.8)
∂ t
∇· E⃗ = 4πρ, (2.9)
∇· B⃗ = 0, (2.10)
where j⃗ is the charge current and ρ charge density.
The assumption of a collisionless plasma is valid when particles are mainly
interacting through fields generated by collective motion, such that indi-
vidual electrons do not significantly experience the fields of other single
electrons. This condition can be quantified by considering a single immo-
bile electron inside a charge-neutral thermalized plasma. The influence of
electron causes the plasma to re-arrange according to the Boltzmann dis-
tribution ne = n̄e exp(qΦ(r )/kb T ), where n̄e is the mean density, q is the
(negative) charge of the electron, Φ is the electric potential generated, kb is
the Boltzmann constant and T is the temperature. At relevant time scales
we assume that the ion number density, ni , is static. Insertion into eq. (2.9)
and using the electric potential E⃗ =−∇·Φ, yields a Poisson equation
∇2Φ=−4πq (n −n ) =−4πq (n − n̄ e qΦ(r )/kb Ti e i e )≈
≈−4πq (ni − n̄e (1+qΦ(r )/kb T )). (2.11)
6 PHYSICS OF LASER WAKEFIELD ACCELERATION
Using the condition of charge neutrality n̄e = ni the equation simplifies to
√
2 4πq
2n̄ √e 4πq 2n̄∇ Φ= Φ→Φ∝ e −r /λD where eλD = (2.12)kb T kb T
where λD is the, so called, Debye length, and the characteristic length scale
for perturbation of single particles. For particle-particle interaction to be
negligible many particles must lie within the Debye sphere, meaning nλ3D ≫
1. As such, the displacement of one single particle will not significantly effect
the generated electromagnetic field.
This condition applies to a wide range of plasmas from astroplasmas, like
the solar wind and the magnetosphere, to laboratory plasmas such as those
generated in fusion reactors and laser wakefield accelerators. Exceptions
are high-temperature plasmas such as the interior of stars, the ionosphere
and gas discharge (such as lightning). Additionally, non-classical phenom-
ena such as quantum effects induced by intense fields (strong-field QED)
are fundamentally not interactions where particle-particle effects can be
neglected. However, by simply adding an additional step in simulations, cap-
turing these effects, they can often be incorporated using similar methods as
for uncorrelated systems. In this thesis, however, I will focus on collisionless,
charge-neutral plasmas satisfying nλ3D ≫ 1.
2.1 Plasma dynamics
In this section, I will outline the dynamics of plasma, focusing on the conser-
vational properties of plasma and waves in plasmas. As stated in the previous
section, the plasma is assumed to be charge-neutral and collisionless. In
addition, when deriving wave modes in plasmas, I will adopt the so-called
cold fluid approximation, in which thermal effects are neglected.
2.1.1 Conservational properties of plasma
The conservative properties of the Vlasov system, the so called fluid equa-
tions, can be derived by taking moments of the Vlasov equation. It should
∫
not come as a surprise that the zero-order moment of eq. (2.6): d f d 3d t p ,
yields the continuity equation, as this relation is inherently linked to the
conservation of density in phase space used to derive the Vlasov equation.
PLASMA DYNAMICS 7
Note, however, that in this case the integral is carried out over momentum
space only as opposed to the full phase space
∫
 ‹
d 3
∂ f p⃗ ∂ p⃗
p + ·∇x f + ·∇ f = 0. (2.13)∂ t m ∂ t p
Each term in the integral is evaluated as
∫
 ‹
3 ∂ f ∂d p = n (t , x⃗ ), (2.14)
∂ t ∂ t
∫ ∫
 ‹  ‹
d 3
p⃗ p⃗
p ·∇x f =∇ 3x · d p f =∇x · (n (t , x⃗ )u⃗ (t , x⃗ )), (2.15)m m
∫
 

d 3p q (E⃗ + u⃗ × B⃗/c ) ·∇p f =
∫
 

= q d 3p ∇p · f (E⃗ + u⃗ × B⃗/c ) = 0. (2.16)
To evaluate the last equation, Gauss’s theorem is applied, and it is assumed
that any physically realistic distribution function vanishes as pi →∞. The
final expression becomes
∂ ∂
n =−∇x · (n u⃗ )→ ρ =−∇x · j⃗ , (2.17)∂ t ∂ t
which is the law of charge conservation. Similarly the fluid momentum
equation can be derived, by taking the first order moment of eq. (2.6) [10].
The conservation law becomes
 
 ‹
∂ u⃗ × B⃗
+ u⃗ ·∇x m u⃗ =−q E⃗ + . (2.18)∂ t c
∫
The second order moment d 3p p
2
2m f gives the energy conservation property
[8]
∂
We k +∇·Q = j⃗ · E⃗ (2.19)∂ t
where We k is the total kinetic energy of the particles and Q is the kinetic
energy flux. The right hand side can be rewritten using eqs. (2.7) and (2.8)
resulting in
 
∂ c E⃗ × B⃗
(W
∂ t e k
+We m ) =−∇·Q −∇· (2.20)4π
8 PHYSICS OF LASER WAKEFIELD ACCELERATION
where W = 1e m 8π (E
2 + B 2) is the electromagnetic energy. Thus the change
in total field energy and particle kinetic energy in a given volume equals
the electromagnetic and kinetic energy flux. Note that this is valid for a
plasma with negligible thermal effects (cold approximation), thus, the en-
ergy conservation law does not include any contribution from pressure or
temperature.
2.1.2 Plasma waves
The waves which can propagate in a plasma is given by the dispersion relation,
it is derived as follows. Starting by applying a Fourier transform to eqs. (2.7)
and (2.8) as
− iωE⃗ = c i k⃗ × B⃗ −4π j⃗ , (2.21)
− iωB⃗ =−c i k⃗ × E⃗ , (2.22)
where, for brevity, any notation for the Fourier transform:
E⃗ (x⃗ , t )→F{E⃗ }(k⃗ ,ω) has been omitted. Inserting eq. (2.22) into eq. (2.21)
and evaluating the triple cross product yields
(ω2− c 2k 2)E⃗ + c 2k⃗ (k⃗ · E⃗ ) =−iω4π j⃗ . (2.23)
From the distribution function f , in eq. (2.6), the current is given as
∫
q
j⃗ (ω, k⃗ ) = d 3p p⃗ f (ω, k⃗ , p⃗ ). (2.24)
m
The resulting equation does not, in general, have an analytical solution, but
in the weak perturbation limit, we can get a bit further by assuming a small
perturbation around a ground state f ≈ f0(p⃗ ) + ϵ f1(t , x , p⃗ ), where F⃗ = ϵF⃗1
[11]. To first order eq. (2.6) becomes
∂ f1 p⃗+ ·∇ f1 + F⃗1 ·∇p f0 = 0. (2.25)∂ t m
Applying a Fourier transform yields a expression for f1 as
 
i E⃗ + p⃗ × (k⃗ × E⃗ )/ωm ·∇
− p
f0
f1 = . (2.26)
(ω− p⃗ · k⃗/m )
PLASMA DYNAMICS 9
Altogether, this results in
(ω2− c 2k 2)E⃗ + c 2k⃗ (k⃗ · E⃗ ) =
 
∫
4πn 20q 3 ωE⃗ + p⃗ × (k⃗ × E⃗ )/m ·∇p f0= d p p⃗ . (2.27)
m (ω− p⃗ · k⃗/m )
For a thermalized plasma f0 is the Maxwell-Boltzmann distribution. However,
for cold plasmas
f0 =δ(px )δ(py )δ(pz ) =δ3(p⃗ ). Evaluating the integral in the cold limit yields
4πn q 2
(ω2− c 2k 2)E⃗ + c 2k⃗ (k⃗ · E⃗ ) = 0 E⃗ . (2.28)
m
This is an eigenvalue equation which can be written in matrix form DE⃗ =
2
ω2p E⃗ , where ω =
4πn0q
p m is the, so called, plasma frequency and, clearly, an
eigenvalue to the matrix equation. Only considering two dimensions the
resulting matrix equation becomes
 
ω2− c 2k 2y −c
2kx ky
2 E⃗ =ω
2 E⃗ , (2.29)
−c kx ky ω2− c 2k 2 px
This equation has a linearly independent solution if det[D−ω2p 1⃗] = 0. This
condition can, after some algebra, be rewritten as
(ω2− c 2k 2−ω2p )(ω
2−ω2p ) = 0 (2.30)
The two modes described by ω2− c 2k 2−ω2p = 0 and ω
2−ω2p = 0 are called
the light-like wave and the plasma oscillation (or Langmuir wave). Insert-
ing these two equations into eq. (2.28) gives the direction of propagation
compared to the electric field
ω2− c 2k 2−ω2p = 0→ k⃗ · E⃗ = 0→ k⃗⊥E⃗ , (2.31)
ω2−ω2p = 0→ k⃗ (k⃗ · E⃗ ) = k
2E⃗ → k⃗ ||E⃗ . (2.32)
Thus we have one electromagnetic wave with transverse direction of propaga-
tion and one with parallel direction of propagation. In magnetized, thermal,
collisional and multi-fluid plasmas there are several other possible waves,
however, to describe LWFA these two types of waves will suffice.
10 PHYSICS OF LASER WAKEFIELD ACCELERATION
2.2 Laser wakefield acceleration
Laser wakefield acceleration utilizes the plasma oscillation described in the
previous section, to create strong accelerating field for particle acceleration.
Laser wakefield acceleration were proposed in the 1970s [3] as a promising
method to accelerate particles up to relativistic speeds. Since then, signifi-
cant progress has been made, taking the technique from idea to realization.
Today, electrons have been accelerated up to 10 GeV [12, 13]. This section out-
lines the theoretical foundations of this technique and discusses its inherent
limitations.
2.2.1 Ponderomotive force
The driving force behind laser wakefield interaction is the ponderomotive
force, which is a force emerging from fast oscillations together with a, slowly
varying, averaged gradient. For example, suppose A⃗ = A0(x||− vp t )cos(k x||−
ωt )ê⊥ is the vector potential of a polarized laser pulse with an envelop A0
traveling at the group velocity vp in direction ê||. By using E⃗ = − 1 ∂ A⃗c ∂ t , B⃗ =
∇× A⃗, eq. (2.18) becomes
 
d u⃗ 1 ∂ A⃗
m = q − + u⃗ × (∇× A⃗) (2.33)
d t c ∂ t
u⃗ can be divided into a transverse (to the direction of propagation of the
laser pulse) and a longitudinal component: u⃗ = u⃗|| + u⃗⊥ and similar with
∇=∇||+∇⊥. Expanding eq. (2.33) in these terms and using
∂ u⃗⊥ d u⃗+ ||
 

+ (u⃗⊥ ·∇⊥) (u⃗⊥+ u⃗||) + u⃗|| ·∇|| u⃗⊥ = (2.34)∂ t d t
 
q ∂ A⃗ q 3∑ ∂ A
− j= +  u j êi − (u⃗⊥ ·∇⊥+ u⃗|| ·∇||)A⃗ . (2.35)m c ∂ t m c ∂ x
i , j=1 i
Assuming that the transverse motion is dominated by the electric quiver
motion, that is m ∂ u⃗⊥∂ t = q E⃗ → u⃗ =−
q
⊥ m c A⃗ and that the pulse width is much
longer than the wavelength such that∇⊥u⃗|| ≈ 0 the equation reduces to
d u⃗ q 2||
m =− ∇ A2, (2.36)
d t 2m c 2 ||
LASER WAKEFIELD ACCELERATION 11
∑
where I used that ∂ (A )A ê = 1∇A2i , j i j j i 2 . Since 〈A2〉 ≈ A20 on longer timescales
and A20 is proportional to the intensity, the ponderomotive force can be in-
terpreted as a radiation pressure.
2.2.2 Wake excitation
Through the ponderomotive force, the laser pulse will push the electrons
ahead and transversely, leaving a positively charge area behind the pulse, and
excite a plasma oscillation, trailing the laser pulse. There are two relevant
regimes of this interaction. In relation to the described ponderomotive force,
it is convenient to introduce the normalized vector potential a = q A00 m c 2 . For
a0 ≪ 1 the ponderomotive force is linear and relativistic effects are negligible.
In contrast, when a0 ≳ 1 the interaction is non-linear and relativistic.
The linear interaction can be derived in one dimension by using eq. (2.36),
with the addition of a longitudinal electric field as
d u⃗ 2|| m c
m = q E⃗ 2
d t ||
− ∇
2 ||
a , (2.37)
where a⃗ = a0A⃗/A0. Using ∂t E⃗|| = −4π j⃗|| = −4πq ne u⃗||, since c∇⊥ × B⃗ = 0 in
one dimension, and assuming ne = n0 +δne ≈ n0 for small linear perturba-
tions, yields
m ∂ 2E⃗ 2|| m c+q E⃗ = ∇ a 2. (2.38)
q n04π ∂ t 2
|| 2 ||
Identifying ω2p and rearranging
 
d 2 ω
2
p m c
2
+ω2p E⃗|| = ∇ a
2 → (2.39)
d t 2 q 2 ||
 
d 2 ω
2
2 p+ω ϕ = a 2, (2.40)
d t 2 p 2
where ϕ = −qφ/m c 2 is the normalized electric potential. This is clearly
corresponding to a plasma oscillation with a source term coming from the
ponderomotive push. The solution can be expressed using Greens functions
as
∫
ω m c 2 tp ∂
E||(x||, t ) = d t
′ sin(ωp (t − t ′)) a 2(x||, t ′). (2.41)q 2 0 ∂ x||
12 PHYSICS OF LASER WAKEFIELD ACCELERATION
Figure 2.1: (a) and (b) shows the wakefield excitation for a one dimensional
cold plasma in the linear and non-linear limit, respectively. Displayed are the
normalized vector potential a = q |A|/m c 2 (blue), the density perturbation
δn/n0 (orange) and the normalized longitudinal electric field E /E0 (green)
where E0 = cωp m/q .
This solution applies in the non-relativistic limit, however, in typical acceler-
ation scenarios a0 > 1. The plasma response in the non-linear (relativistic)
limit is given by [1, 14]
c 2 ∂ 2ϕ (1+a 2) 1
2 = 2 − . (2.42) 
ωp ∂ ξ2 2 1+ϕ 2
The expression is derived using comoving-coordinates ξ= x||− c t and τ= t
such that ∂t →−vp∂η+ ∂τ, where vp is the group velocity of the laser pulse,
and ∂x|| → ∂η. Furthermore, the relativistic motion of the particles must
be taken into account. Therefore, the expression for the ponderomotive
force derived above must be modified. In this frame the interaction can
be considered static, thus the derivative with respect to τ is zero and a is
constant in time τ. This is called the quasi-static approximation in the limit
vp ≈ c , corresponding to a low density plasma, which is typically the medium
LASER WAKEFIELD ACCELERATION 13
for LWFA. In fig. 2.1 the linear and non-linear wake excitation is shown. In
the non-linear relativistic regime the wave steepens and takes on a saw-tooth
form, notably, the wavelength elongates.
2.2.3 Physical limitations
For LWFA there are a couple of important scalings that influence the maxi-
mum energy gain possible.
Maximal field
Primarily, the wave breaking limit influence the maximal field that a plasma
can sustain. Heuristically, this corresponds to the point at which the elec-
trons in the strongest field move faster than the wave itself. In the linear
limit the maximal amplitude can be estimated by assuming E||∝ sin(kp x −
ωp t ) where kp = ωp/vp ≈ ωp/c and using eq. (2.9) with t = 0 as Emax =
4πq n0 cω m
k =
p
q . In the relativistic limit the corresponding limit is Ep max,rel =
Æ
2(γp −1)Ema x , where γp = (1− v 2 2 −1/2p /c ) [1]. The phase velocity of the
plasma wake is the group velocity of the laser pulse that is vp = vg ,l = ∂k (ω(k )).
While the wakefield consists of plasma oscillations, the laser propagates by
transverse light-like plasma waves, i.e the derived dispersion relation in
eq. (2.31). Therefore
∂ q
v 2 2 2
c
p = c k +ωp = q . (2.43)∂ k 1+ω2p/c 2k 2
When γp ≫ 1, one can approximate γp ≈ c k/ωp , with this approximation
Æ Æ cω m
E = 2(c k/ω −1)E ≈ 2c k/ω p ∝ n 1/4max,rel p ma x p q 0 . Thus by increas-
ing the density the accelerating field can be increased. Notably, this is differ-
ent from the linear approximation, where the maximal electric field scales
p
like Ema x ∼ n0.
Dephasing
When the accelerated electrons have gained enough energy, they will, eventu-
ally, outrun the wake. The electrons start to decelerate roughly after covering
14 PHYSICS OF LASER WAKEFIELD ACCELERATION
half a plasma wavelength relative to the wake. Assuming the electrons move
at the speed of light, we have
2 c
(c − vp ) = . (2.44)λp Ld
Thus,
λp/2 2πcω2
Ld =
−3/2
q ≈ 3 ∝ n0 . (2.45)
1− 1−ω2p/ω2 ωp
Since Wmax ≈−q EmaxLd , or in the relativistic limit Wmax ≈−q Emax,relLd , this
length sets the limit for the maximum achievable energy gain.
Pump depletion
In considering the limitations of LWFA the depletion of the driving laser pulse
should be taken into account. The length at which the driver is depleted can
be estimated by equating the energy of the wake at some time td p with the
∫ L ∫
initial energy of the pulse: p d E 2|| d x|| = E
2
l d x||. For a gaussian pulse, in0
the quasi-static approximation, this becomes [15]
¨
cω2 2/a 2, a ≪ 1
L 0 0p d ≈ 8.7 3 (2.46)ωp 1, a0 ≫ 1.
Therefore, as a rule of thumb, L −3/2d p ∝ n0 . p
For this reason, and since the maximal energy only scales as Emax ∝ n at
best, high energy electrons are typically accelerated using long acceleration
lengths and low density plasmas. However, there is a second limitation: the
diffraction of the laser pulse. The diffraction length is typically on the order
of the Rayleigh length, however, relativistic self-focusing in the plasma can
significantly extend the distance over which the laser remains focused [1].
Nonetheless, to accelerate electrons to energies ≳ 1GeV, optical guiding is
generally employed [12, 13].
2.2.4 Electron trapping
So far, the discussion has focused only on wake excitation and propagation.
However, to accelerate electrons, they must also be properly positioned to
LASER WAKEFIELD ACCELERATION 15
become trapped by the wave. To consider this, the one dimensional descrip-
tion lacks detail, and since the three dimensional analytical descriptions in
this regime is limited, this section considers simulations to explain the LWFA
interaction. In 3D, the laser pushes the electrons forward and longitudinally,
which creates a current around a positively charged area behind the pulse.
The electrons meet up at the back of the wake and create the density peaks
show in fig. 2.1.
Figure 2.2: (a) shows the normalized transverse vector potential of the laser
pulse. (b) shows the longitudinal normalized electric field, E /Ema x , where
Ema x = cωp m/q . (c) shows the transverse electric field. To become trapped
and accelerated electrons needs to reside in the blue region of panel (b). The
fields were generated though PIC simulations.
This motion generates not only the longitudinal field described in the pre-
vious sections but also a transverse focusing field, as shown in fig. 2.2(b–c).
Consequently, to become trapped in the wake an electron has to be placed
at the back of the wake (blue region in fig. 2.2.(b)) to experience a simultane-
ously transverse focusing force and longitudinal accelerating force.
Furthermore, to be come trapped, an electron needs to have the same ve-
locity as the local phase velocity of the wake. In a thermal plasma this can
16 PHYSICS OF LASER WAKEFIELD ACCELERATION
happen spontaneously, or be triggered by wave breaking. These processes
are however, random and uncontrolled. A more controlled injection method
is downramp injection [16], where the phase velocity of the wave is decreased
by a density downramp. When the density decreases the plasma wave elon-
p
gates, since λp ≈ωp/c ∝ n0 , thus the velocity of the density peak relative
to the laser pulse is decreasing, momentarily allowing more electrons to
enter the trapped orbits of the wake.
However, with this scheme electrons are injected over a wide range of initial
positions in the wake, causing them to experience different accelerating
fields, which results in large energy spread. Minimizing energy spread re-
quires that the electrons are injected at closely matched positions and times
in phase space. A combination of down-ramp injection and ionization injec-
tion [17] provides a good method to achieve this.
Figure 2.3: Shows particle orbits where green is trapped orbits. The orbits
were generated using the non-linear one dimensional theory for Hamiltonian
particle orbits.
In ionization injection electrons are injected from zero momentum by further
ionization of the medium. To study the injection, particle orbits can be
derived by considering the one dimensional Hamiltonian of a particle, in a
frame traveling with the wave
q
H = c m 2c 2e + (p⃗ −q A⃗)2 +qφ− vp p||. (2.47)
Using the expression for ϕ =−qφ/m c 2 in eq. (2.42), results in the particle
orbits shown in fig. 2.3. There exists a region close to the laser in which
LASER WAKEFIELD ACCELERATION 17
particles with zero initial momentum are able to become trapped, see fig. 2.3.
However, this region does not overlap with the peak of the laser, where ioniza-
tion of the inner-shell electrons occurs. By introducing a density downramp,
the trapped particle orbits shift forward, providing a window for the newly
ionized electrons to be captured.
Additionally, the injected electrons are themselves influencing the wakefield
by generating a secondary wake. This phenomenon is called beamloading
[18, 19]. The injected electrons pushes the electrons streaming around the
pulse outward by the coulomb force, locally weakening the electric field.
A high injected charge locally inverts the electric field, causing the trailing
electrons to accelerate less than the leading electrons, resulting in large
energy spread. In contrast, an insufficient charge also leads to large energy
spread, since the unperturbed accelerating field is positively correlated. An
optimal charge distribution exists in between these two extremes, which
flattens the accelerating field, ensuring that all electrons experience uniform
acceleration and minimizing energy spread[19].
18 DATA-DRIVEN MODELING OF LASER WAKEFIELD ACCELERATION
3 Data-driven modeling of laser wakefield ac-
celeration
As discussed in the preceding chapter, the dynamics of particle acceleration
in laser wakefields are highly non-linear and influenced by several interacting
processes. The laser wavefront, energy, spectrum and beam profiles as well
as variations in plasma density may all affect the wake excitation, injection
and subsequent acceleration. Furthermore, in experimental setups, these
parameters often show significant shot-to-shot fluctuations, exceeding the
typical sensitivity of LWFA. For instance, peta-terawatt laser system can have
focus fluctuations on the order of µm-mm [20, 21], while optimal injection,
happens over the course of a few wavelengths (∼ 10µm), for minimal en-
ergy spread and optimal beamloading. In downramp-ionization injection,
any variation in the laser intensity within this region strongly affects the
properties of the accelerated electrons.
An intriguing application of LWFA is to use it as a seed for undulators for
generating X-rays. Using the LWFA technique would substantially lower the
costs of applications, such as molecular imaging, today carried out using
km-long linear accelerators. However, the instability of todays laser-based
accelerators still lacks compared to conventional accelerators[4].
For optimization and stabilization of these systems, it is essential to under-
stand how variations in system parameters affect the acceleration process.
Large-scale statistical modeling of these systems has long been limited by the
sparsity of data in the field. However, recent advances in long term stability
[21, 22] and advancement towards high-repetition-rate systems, pave the
way for data-driven modeling.
Paper A presents a study on how laser properties in an LWFA experiment
influence the electrons and the resulting X-rays produced when the acceler-
ated electrons are guided through an undulator. Specifically, a Feedforward
Neural Network (FNN) model is developed to reduce biases caused by mea-
surement errors in the laser parameters. The standard regression models
(including FNN), employs an optimization criteria based on the difference
between the model and the dependent variable. This is adequate when there
is only measurement error on the output of a system, however, will introduce
bias when there is an additional error on the input. To mitigate the bias one
has to employ a minimization criteria based on the orthogonal distance be-
PREDICTION AND ESTIMATION 19
tween data and model. This approach is very similar to orthogonal distance
regression, however, so far not widely implemented in the context of neural
networks. In this work, I develop such an optimization criteria and apply it
the problem of modeling the dependence between laser, electron and X-ray
properties.
For many applications of X-rays, such as absorption and emission spec-
troscopy, the knowledge of the precise X-ray spectrum radiating the sample
is critical [23], in particular when the spectrum exhibit large shot-to-shot vari-
ations. In such cases, a so-called virtual diagnostic, where a model predicts
the spectrum through non-invasive measurements, can ease the constraints
on stable operation. This work presents a model achieving high accuracy
prediction of the X-ray spectrum using laser and electron diagnostics.
3.1 Prediction and estimation
When modeling a process f , characterized by some input x and output y ,
there can be two distinct motivations. Either, one is interested in achieving
a good prediction of y given x or one is interested in a good estimation
of some internal variable of the process θ such that f (x ,θ ) = y . In some
cases, these two objectives align, however when there is a measurement
error on x the distinction becomes important. Here, this framework is used
to create two FNN models: an explanatory model and a predictive model.
The explanatory model uses a optimization criteria based on the orthogonal
distance, described in the following section, to mitigate the influence of
errors in input variables and, thus, achieves a more correct estimation of the
dependence between laser, electron and X-ray properties. The predictive
model employs a FNN to estimate the X-ray spectrum from the laser and
electron parameters, with its primary application being virtual diagnostics.
3.2 Orthogonal distance loss
To illustrate the problem of measurement error modeling a synthetic exam-
ple is show in fig. 3.1. Here, two networks are trained: one on minimizing
∑
the mean squared error between model and data ( | f (x )− y |2) and one
on minimizing the orthogonal distance between model and data. By not
taking the measurement error of the input into account the model based on
20 DATA-DRIVEN MODELING OF LASER WAKEFIELD ACCELERATION
minimizing the mean squared error is biased towards lower dependencies
between input and output.
Figure 3.1: Shows a example problem where two models (blue and orange) are
trained to replicate the underlying function (black). One model was trained
using the orthogonal distance loss (OD) and the other using mean square error
(MSE). The data was generated with equal error variance on x and y . The 2D
histogram shows the data distribution and the black dots randomly selected
points of this distribution.
There exists numerous models which corrects for measurement errors [24–
27], where the most well used is Total-Least-Square regression [28]. However,
most implementations set limitations on the base functions used or is only
implemented in lower dimensionality [29]. Implementations for neural net-
works are still largely missing. Here, I present a implementation, in parts
based on the method presented in [29].
The measurement error problem can be mathematically expressed as
xi = x
∗
i + ϵx ,i
yi = f (x
∗
i ) + ϵy ,i (3.1)
where x ∗i are the true measurement values and ϵx ,i ,ϵy ,i are noise parameters
with the same dimensionality as xi , yi . Note that, in general xi , yi has the
dimensionality of the number of output and input respectively. The goal is to
EXPERIMENTAL SETUP AND DATA SELECTION 21
find a FNN fv with weights v that approximates f . To do this a proxy is used
to first model x → x ∗. This proxy is constructed by a second neural network
fw (x , y ). The optimization criteria becomes (for derivation see Paper A)
 2  2
∑ yi j ′ − fv ( fw (xi , yi )) j ′ xi j − fw (xi , yi ) j
L (xi , yi ) = min + , (3.2)v,w σ2 σ2i , j , j ′ y , j ′ x , j
where σx , j ′ ,σy , j are the standard deviations corresponding to ϵx ,i ,ϵy ,i , re-
spectively and j , j ′ runs over all input and output, respectively. To optimize
fv , fw given data pairs xi , yi the gradient with respect to v, w is and used in
Stochastic Gradient Decent (SDG) (specifically Adam). In practice, this is
done in two steps, where fv is first optimized with one SDG step and, sequen-
tially, fw with one step. Note that, since the weights w features in both terms,
fw is prevented from simply taken the form fw (x , y ) = x , which would be
similar to not having this term in the optimization criteria.
For the explanatory model this architecture was used to extract a model that
approximate the underlying function relating laser properties to output of
the system.
3.3 Experimental setup and data selection
The analysis was performed on a data set from the LUX beamline (fig. 3.2)
from a day-long measurement, comprised of up to 80, 000 data points. The
high-power laser system Angus provided a 2.6 J laser pulse at 1 Hz repetition
rate. The laser was focused by an off-axis parabola (OAP) into a plasma target
with three gas inlets. By injecting different mixtures of hydrogen and nitrogen
into the inlets a tailored plasma density distribution, through the ionization
of the laser pulse, was achieved. The electrons were injected during the den-
sity downramp, and accelerated during the subsequent density plateau. After
acceleration, the electrons were focused, using a quadropole doublet, into
the BEAST II undulator where they emit spontaneous undulator radiation in
the X-ray range. After the undulator the electrons are deflected by a dipole
spectrometer, measuring the electron spectrum. The X-rays were focused
using a grazing mirror onto a transmission grating which disperse the radi-
ation onto a X-ray Charge Coupled Device measuring the X-ray spectrum.
Along the laser beamline, diagnostics were collecting measurements of the
22 DATA-DRIVEN MODELING OF LASER WAKEFIELD ACCELERATION
Figure 3.2: The laser was focused into the plasma target, where variable laser
gas settings creates a tailored plasma profile. The accelerated electrons were
subsequently focused by a quadrupole doublet into an undulator where they
emitted spontaneous undulator radiation. The electron spectrum was mea-
sured downstream of the undulator using a dipole spectrometer, while the X-ray
radiation was dispersed by wavelength using a transmission grating and sub-
sequently recorded by a Charge-Coupled Device (CCD). Along the beamline
several diagnostics were monitoring the laser properties, including a spec-
tral measurement, energy sensor, wavefront sensor and beam far-field and
near-field profile. Additionally, the electron beam position and charge were
non-invasively measured using two beam position monitors (BPMs).
laser energy, spectrum, beam profile and wavefront. Additionally, the elec-
tron position and charge were monitored using two Beam Position Monitors
(BPMs).
The laser diagnostic system provided up to 70 different parameters. To build
an accurate model capable of identifying correlations between specific laser
parameters and system output, a subset of these parameters needed to be
selected. To determine which parameters that influence the acceleration
most, permutation feature importance [30], was employed. Furthermore, to
avoid multicollinearity among the input variables, in order to allow indepen-
dent variation of each variable in the analysis, the variance inflation factor
was used. This is a measure of mutual correlation within a set of variables,
ELECTRON AND X-RAY DEPENDENCE ON LASER PROPERTIES 23
calculated as VIF = 1/(1−R2i ) where R
2
i is the R2-score of a linear regression
model that predicts the i :th variable based on all other variables. The se-
lection procedure resulted in a set of 14 laser parameters. The explanatory
model used only this set of parameters as input, while the predictive model,
in addition to the laser parameters, used the electron data: BPM measure-
ments and electron spectrum as input. For training and validation the X-ray
and electron spectrum were decomposed into 20 variables using principal
component analysis.
3.4 Electron and X-ray dependence on laser properties
After training the explanatory model on the experimental data, it was ana-
lyzed by varying one or more input parameters, as shown in fig. 3.3, while
keeping the remaining inputs constant at their sample mean. This approach
allows the effect of individual parameters to be examined, which would
otherwise be hidden due to the high dimensionality of the dataset.
The analysis showed that there is a strong connection between the laser focus
position and the electron spectrum through the process of beamloading. The
focus position influences laser intensity in the injection region (at the plasma
downramp), which results in an increase or decrease of injected electrons. A
low degree of injected charge, leads to a positively correlated accelerating
field, in turn generating a negative correlation in the longitudinal electron
phase space, leading to a large energy spread, and a large mean energy.
On the other hand, a high degree of injected charge leads to a negatively
correlated electron field and a positive correlation of electron energy phase
space. The optimal energy spread is found in between these two extremes,
where the generated wakefield of the bunch cancels the positive correlation
of the acceleration field. This dynamic explains the correlation between
focus position and electron energy E , charge Q and energy spread ∆E in
fig. 3.3.(a-c).
Furthermore, a good overlap was found between the model surface as a func-
tion of focus position and laser energy when compared with PIC simulations,
see fig. 3.3.(d-i). This procedure can be used in tuning the system to find
global optimal points outside the range of the explanatory model.
By analysis of the correlation between laser focus position, laser central wave-
length and electron energy spread, a source of measurement error could be
24 DATA-DRIVEN MODELING OF LASER WAKEFIELD ACCELERATION
Figure 3.3: (a-c) shows the charge (Q ), electron energy (E ) and energy spread
∆E as a function of focus position (z f o c ) for the explanatory model (blue) and
PIC simulations (black). Additionally, (a-c) shows the distribution of the ex-
perimental data. Panel (d-f) shows the function surface extracted from the
explanatory model as a function of laser energy and focus position. This func-
tion is compared by the function surface in the same dimensions, given by PIC
simulations in panel (g-i). This surface were generated by randomly sampling
points in the space, performing a PIC simulation and then modelling the result-
ing surface using a cubic spline. The maximal overlap of the surface given by
the explanatory model and the PIC simulations was found in the region marked
by the red box in panel (g-i).
identified. Laser central wavelength is directly influencing the focus position
by chromatic aberrations of the optical elements along the beamline. These
effects are included in the wavefront measurement, however, additional
chromaticity introduced by the optical elements in the wavefront diagnos-
PREDICTION OF X-RAY SPECTRUM 25
tics, behind the OAP, is not present in the focus position delivered inside the
plasma target. The resulting measurement error was identified by comparing
the optimal focus position, given by the explanatory model, with the depen-
dence on laser central wavelength and focus position. Furthermore, the
chromaticity of the wavefront diagnostics, calculated through ray-tracing,
showed good agreement with the chromaticity predicted by the model.
For higher order wavefront aberrations of the laser, such as astigmatism and
coma, analysis of the model indicated that they affect the electron accel-
eration similarly as changes in focus position, that is, through changes in
effective intensity in the injection region leading to different injection and
beamloading.
Additionally the explanatory model showed that focus position had a strong
correlation with generated X-ray spectrum. This result is expected, as undula-
tor radiation depends directly on the electron energy spectrum. Nonetheless,
the clear model relationship between the focus position, which can be ac-
tively controlled in the beamline, for example by using a beam expander, and
the resulting X-ray spectrum provides a good foundation for tuning system
output in future experiments.
3.5 Prediction of X-ray spectrum
The predictive model achieved high accuracy predictions, with R2-scores as
high as 0.98 for the prediction of the central X-ray wavelength. This shows
promising potential for using the model as a non-invasive diagnostic. In
fig. 3.4 the model performance over 50 consecutive shots from the validation
set is shown. This model used the electron spectrum, BPM measurements
and the 14 selected laser properties as input. Additional, two other models
were trained to test the performance when the electron spectrum was not
provided as an input. The first model used only the laser properties as input
and the second model used the laser parameters as well as the BPM measure-
ment. These models achieved a R2-score of 0.61 and 0.93, respectively. The
competitive R2-score of the second model suggests that it could be employed
in scenarios where electrons are involved in subsequent experiments, as the
model would then remain entirely non-invasive.
26 DATA-DRIVEN MODELING OF LASER WAKEFIELD ACCELERATION
Figure 3.4: (a-c) shows the predicted, measured and residual error of the X-ray
spectrum over 50 shots of the validation set for the predictive model.
3.6 Conclusion and outlook
This study showed that data-driven models, when carefully designed, can
be used to extract the complex dynamics that governs the LWFA process.
By analysis of the multidimensional explanatory model, dependence be-
tween focus properties and electron and X-ray properties were extracted,
uncovering correlations, previously hidden by the high dimensional space
and identifying sources of measurement errors. Moreover, the predictive
model showed high accuracy predictions of the X-rays spectrum, which in
the future can be used to provide prospective users valuable details on the
X-ray spectrum. In addition, both models can be used to provide a base for
active feedback in future high-repetition rate systems.
CONCLUSION AND OUTLOOK 27
Laser wakefield acceleration is a compact way of generating relativistic elec-
trons for secondary generation sources. However, their operational stability
and electron quality lacks compared to conventional linear accelerators. The
presented models can be used as tools to identify sources of variations and
stabilize schemes by feedback loops, to create more stable plasma-based
accelerators.
With the recent push toward high-power high-repetition-rate laser systems,
the growing volume of available data will further enhance the role of data-
driven methods in the field of LWFA. Combined with advances in high-
performance computing, these developments will be crucial in transitioning
plasma-based accelerators from experimental demonstrations to practical
applications.
28 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
4 Particle-in-cell modeling of laser wakefield
acceleration
The particle-in-cell (PIC) method was introduced to the field of plasma
physics in the late 1950s and early 1960s by John Dawson and Oscar Bune-
man among others [31]. Since then it has evolved into an invaluable tool for
studies of a variety of plasma phenomenons ranging from astrophysics [32]
to strong-field quantum electro dynamics [33].
The widespread adoption of the PIC algorithm can be largely attributed to its
versatility. PIC simulations have the ability to model a broad range of plasma
regimes while consistently incorporating physical processes beyond basic
plasma theory, such as ionization, radiation reactions, external fields, and
boundary effects.
However, one should keep in mind that PIC codes, through the discretiza-
tion, represent a reduced model of reality and are therefore subject to various
numerical limitations. PIC simulations are constrained by stability and ac-
curacy conditions that depend strongly on the chosen resolution and the
system’s physical parameters (e.g., plasma frequency). However, there is
also numerical inaccuracies that cannot be reduced simply by enhancing
resolution, but may be addressed with alternative discretization methods.
Such an instability is Numerical Cherenkov instability (NCI) [7]. This instabil-
ity is a persistent challenge to the laser wakefield community since it hinders
a potential speed up of LWFA simulations by performing the simulations in
a Lorentz boosted frame [5, 34]. In this frame the plasma is contracted and
laser wavelength elongated, as illustrated in fig. 4.1. Since the minimal time
and space step are typically restricted by the laser wavelength, this transfor-
mation enables a decrease in resolution and thereby reduces computational
costs.
Mitigating strategies of numerical Cherenkov instability (NCI) entails filter-
ing of fields in coordinate space (digital filtering) [35, 36] and in spectral space
[37–39]. Arguably, a less invasive solution is provided by using co-moving
coordinates in a Galilean transform [40]. These efforts have, to a large degree,
made boosted-frame simulations feasible, though some caveats remain. For
instance, both digital filtering and current scaling can introduce unwanted
damping of physical phenomena. Additionally, Galilean transformations
only mitigate the NCI for particles moving with the Galilean-transformed
29
Figure 4.1: In a Lorenzt-boosted frame (γ = 3) the wavelength of the pulse
denoted as the blue-red colormap (arbitrary scale) is elongated and the plasma
(gray) is contracted enabling a decrease in resolution.
frame. In this thesis, I investigate a different approach based on the conser-
vation of physical properties.
Numerical stability frequently manifests as violations of physical principles.
This occurs because the discretization breaks underlying symmetries, like
spatial or temporal translation symmetry. While the connection between
symmetries and conservation laws (via Noether’s theorem) only holds for
continuous symmetries, it is reasonable to expect that breaking a symmetry
through discretization alters or weakens the associated conservation law.
In this section, I compare the NCI growth of three different PIC algorithm: (i) a
standard PIC algorithm with a spectral solver [41] and a Boris particle pusher
[42], (ii) an energy conserving solver [43] and (iii) the energy conserving
solver together with a momentum correction (described in Paper B). To
compare different PIC algorithms, the modular PIC code π-PIC was used.
While most PIC codes share a common structure, few offer the flexibility
needed to easily incorporate and test different algorithmic approaches. The
π-PIC framework is intended as a step toward providing such flexibility and
supporting decentralized development of PIC algorithms.
The content of this section is ordered as follows. An overview of the PIC
algorithm is presented in section 4.1, followed by a discussion of the related
30 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
stability criteria in section 4.2, where additionally the connection between
numerical instabilities and conservation laws is examined in greater detail.
Theπ-PIC framework is described in section 4.3. Finally, in section 4.4,π-PIC
is used to study how the NCI is affected by different PIC implementations.
4.1 The PIC algorithm
The evolution of a collisionless plasma is governed by eq. (2.6). However,
when electrons reach relativistic velocities, the relativistic form of this equa-
tion must be used
∂ f p⃗
+ ·∇ f + F⃗ ·∇p f = 0, (4.1)∂ t γ
where f (x⃗ , p⃗ , t ) is the distribution function and F⃗ = q (E⃗ − p⃗ × B⃗/(γm c ))
is the Lorentz force. The electromagnetic field (E⃗ , B⃗ ) evolves according to
Maxwells equations eqs. (2.7) to (2.10). This system can, in principle, be
solved numerically on a grid, however, the phase-space dependence of f
is six dimensional, which creates large computational costs. Nevertheless,
such methods do exist and are mainly used in low-resolution or reduced-
dimensionality implementations, commonly referred to as Vlasov codes [44].
A more computationally efficient option is to reduce the number of particles
by considering macro-particles (or quasi-particles) as
Np
∑
f (x⃗ , v⃗ , t ) = wi S (x⃗i (t )− x⃗ )δ(v⃗i (t )− v⃗ ), (4.2)
i
where wi encodes the number of physical particles per macro-particle, Np
is the number of macro-particles and S is a shape-function encoding the
distribution of charge within the macro-particle (with dimension 1/L3). In
general, the PIC algorithm consist of four steps, shown in fig. 4.2 and detailed
in the upcoming paragraphs.
THE PIC ALGORITHM 31
Figure 4.2: The PIC algorithm can be divided into four steps: (i) particle push,
where the particles are advanced using the Lorentz equation, (ii) current and
density deposition, where the particle position and velocity is interpolated to
current and density, (iii) field advancement and (iv) field interpolation, where
the electromagnetic field is calculated at the position of the particles.
The particle push
The particle movement is governed by the equation of motion of particles in
electromagnetic fields, namely
 
∂ p⃗ v⃗ × B⃗
= q E⃗ + , (4.3)
∂ t c
∂ x⃗
= v⃗ , (4.4)
∂ t
p
where p⃗ = γm v⃗ and γ(v⃗ ) = 1/ 1− v 2/c 2 is the Lorentz factor. The most
used discretization scheme is the Boris pusher [42], which advances the
32 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
particles by solving the equations
p⃗ n−1/2− p⃗ n+1/2  ‹p⃗ n
= q E⃗ n + × B⃗ n , (4.5)
∆t m c
x⃗ n+1 = x⃗ n +∆t v⃗ n+1/2, (4.6)
where the superscript denotes time indices. The term p⃗ n is constructed as
p⃗ n γn+1/2v⃗ n+1/2 +γn−1/2v⃗ n−1/2
= , (4.7)
m 2c γn
where γn = (γn−1/2 +γn+1/2)/2. The half timestep can be implemented by at
first iteration advancing the velocity half a timestep using a forward differ-
encing scheme: v 1/2 = v 0 + ∆t q2m (E
0 + v 0×B 0).
Current deposition
In this step the current density is interpolated from the movement of the
macro-particles on the grid. Since the density velocity is given by eq. (2.3)
the discretized form of the current becomes
Np
∑
j⃗ (x⃗ ) = q v⃗i wi S (x⃗i − x⃗ ). (4.8)
i
In some codes, the charge density ρ is deposited in a similar manner. In
general only the current is needed to advance the field using eq. (2.7) and
eq. (2.8), therefore, the charge density is mostly used when codes employ
divergence cleaning to ensure charge conservation. Other codes uses, so
called, charge conserving deposition schemes. This topic will be discussed
further in section 4.2.3
Field advancement
The field advancements can be done using a variety of field solvers. In the
context of PIC simulations, the two most commonly used classes of field
solvers are spectral solvers and finite-difference methods. Spectral solvers
advance the fields in time within spectral space, whereas finite-difference
methods employ discretization schemes, such as the Finite-Difference Time-
Domain (FDTD) method, to evolve the fields. The topic of choosing field
solvers is further discussed section 4.2. Additionally a description of an
electrostatic field solver is given in section 4.1.1.
THE PIC ALGORITHM 33
Field gathering
The last step in the PIC update is the gathering of fields to particle locations.
Since the fields are defined on a grid nodes x⃗ ′j they are interpolated to the
particle locations x⃗i using the shape function as
Nn
∑
E (x⃗i ) = S (x⃗ − x⃗ ′i j )E (x⃗
′
j ), (4.9)
j
Nn
∑
B (x⃗i ) = S (x⃗
′ ′
i − x⃗ j )B (x⃗ j ), (4.10)
j
where Nn is the number of adjacent grid nodes considered.
Another important aspect to consider is the implementation of the algo-
rithms on modern computer architectures. Simulating realistic plasma in-
teractions can require several days, even with supercomputers. Therefore
the evolution of PIC algorithms has greatly benefited of the development
of computing power which together with more efficient implementations
has pushed the boundaries for systems which are feasible to analyze via
simulations.
4.1.1 An electrostatic PIC algorithm: a detailed example
As an illustrative example of the design of a PIC algorithm, this section out-
lines the construction of an electrostatic solver, providing the reader with an
overview of the considerations and the limitations associated with a specific
algorithm.
The electrostatic assumption states that the time derivative of the magnetic
field is negligible. Consequently, Maxwell’s-Faraday’s law, eq. (2.8), reduces
to
∇× E⃗ ≈ 0, (4.11)
while, Maxwell’s-Ampere’s law, eq. (2.7), remains unchanged. Similarly, the
electrostatic assumption states that the push of electrons is solely electric.
Consequently, an electrostatic one dimensional PIC algorithm can be imple-
34 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
mented by advancing the state according to
∂
p
∂ t x
= q Ex , (4.12)
∂ Ex +4π jx = c (∇× B⃗ )x . (4.13)∂ t
Using a simple finite difference method (namely Forward-Euler) the advance-
ment of Ex can be implemented as
€ Š
j+1 j j+1
Ex ,i = Ex ,i −∆t 4π jx ,i +Ci (4.14)
where upper indices denotes time indices and lower denotes space indices,
moreover, since B⃗ is constant in time, the contribution from the magnetic
field has been grouped into Ci . The particle equation of motion can be
implemented using a forward leap-frog differencing scheme
j j
v j+1/2
∆t q E
= v j−1/2 + x
(x )
, (4.15)
m
x j+1 = x j +∆t v j+1/2. (4.16)
With these two advancement strategies, all that remains is the interpola-
tion from particle distribution to the current on the grid. The current is
interpolated to the grid from the particle position and momentum as
Np
∑
j+1
jx ,i = w qS (x
j+1 j+1/2
p − xi )vp , (4.17)
p
where S is the shape function, w is the macro-particle weight and the lower
index denotes particle index or grid point. Here, we use the Cloud-In-Cell
shape function which corresponds to a uniform distribution of charge in
each macro particle with the extent of one grid cell. The electric field is
interpolated back to the particle position similarly
Ng
∑
E jx (x ) = S (x − x
′
i )
j
Ex ,i , (4.18)
i
where Ng are the number of adjacent grid nodes to x
′. With this simple set of
equations is is possible to simulate quasi-electrostatic phenomenons such
THE PIC ALGORITHM 35
as plasma oscillations and the initial evolution of charge-current neutral
two-stream instability as shown in fig. 4.3.
However, this description have its limitations. To consistently simulate elec-
trostatic systems Gauss equation, eq. (2.9), has to be satisfied, which, together
with eq. (4.11) leads to Poisson equation as
«
∇· E⃗ = 4πρ 2
−1→D
∂ φ
4πρ+ = 0. (4.19)
∇× E⃗ = 0↔ E⃗ =−∇φ ∂ x2
If the simulation is initialized in a state that is consistent with these equations
then the state will conserve the relation to numerical accuracy (see Paper
B). The eq. (4.19), implies that in a periodic space the total charge in the
∫
simulation box must sum up to zero, since ∇ · f d 3 x = 0 for an arbitrary
periodic function f . Moreover, integrating eq. (4.13) yields
∫ ∫ ∫ ∫
∂ 1 ∂ 
Ex d x + 4π jx d x = (−∂ φ)d x + j d x =∂ t 4π∂t x x
∫ ∫

= c (∇×
B⃗ )x d x = 0→ jx d x = 0, (4.20)
where, Stokes (integral) theorem, was used to evaluate the integral over∇× B⃗ .
In summary: to satisfy eq. (4.19) (in a periodic space) both the spatial average
of current and charge over the simulation box has to be zero. Consequently,
the solver is also restricted to these plasma states. The effects of these lim-
itations are visible in the evolution of the two-stream instability in fig. 4.3
simulated using the electrostatic and a electromagnetic algorithm. Initially
the phase space of the electrostatic and electromagnetic solver are identical,
however, as the simulation progresses they begin to deviate. Simultaneously,
the total energy of the electrostatic solver start to increase rapidly. One can
assume that these deviations are connected to the different descriptions of
physics: while the electromagnetic solvers allows time-varying fields and
net charge, the electrostatic solver does not. It should be noted that neither
of these algorithms conserves energy, although, the electromagnetic solver
performs considerably better in this regard than the electrostatic one.
36 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
Figure 4.3: Two-stream instability simulated with the electrostatic and a elec-
tromagnetic solver. The electromagnetic solver uses a spectral field solver and
the Boris particle pusher. The solvers initially simulate similar dynamic (a-b)
but when the instability grows they deviate (c). The growth of total energy
signals the breakdown of the electrostatic approximation, as shown in the bot-
tom figure, where the times corresponding to the plasma states in (a–c) are
indicated by green lines.
NUMERICAL STABILITY 37
4.2 Numerical stability
While the previous section described the limitations of a PIC algorithm based
on the used description of physics, or its regime of validity, this section focus
on the limitations which are inherent to the numerical scheme itself. This is
the topic of numerical stability.
More specifically, numerical stability refers to an algorithm’s capacity to
yield consistent results when subjected to small perturbations in the initial
conditions. Stability requires that numerical errors remain bounded over
time. Conversely, if they grow or oscillate without limit, the algorithm is said
to be numerically unstable. Naturally, stability is a necessity for any useful
algorithm.
As previously stated: numerical stability frequently manifests as violations
of physical principles. In the following sections, sections 4.2.2 to 4.2.5, I
elaborate on this statement by outlining several common types of numerical
instabilities, discussing their connection to physical principles and symme-
tries, as well as mitigation strategies. In particular, section 4.2.5 presents a
derivation of the numerical dispersion relation used to characterize the NCI.
4.2.1 The Courant-Friedrichs-Lewy condition
For simulation of electromagnetic fields a fundamental stability constraint
is the Courant-Friedrichs-Lewy (CFL) condition. This condition states that
the numerical domain used for propagating any field point in time must
include the physical domain, or, put differently: during one timestep a signal
should not be able to travel between two spatial grid points which are not
coupled through the numerical scheme. For FDTD, which uses at most the
neighboring points (in space and time) for evaluating the field (in 1D) this
translates to
∆x
∆t ≤ . (4.21)
c
This condition is a minimal requirement for field solvers using finite differ-
encing schemes, meaning that there might be other less relaxed conditions
for achieving a stable algorithm.
For spectral solvers there are no explicit CFL-condition since the fields are
solved on a global basis in Fourier space [40]. However, the simulations must
still resolve the relevant wavelengths and frequencies, which, for simulations
38 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
of electromagnetic waves in vacuum or low density plasma, ultimately leads
back to the constraint given by eq. (4.21).
4.2.2 Energy and momentum conservation
A part from the CFL condition, the spatial resolution of most PIC algorithms
are constrained by the Debye length. If the Debye length is not resolved
aliased modes will appear as thermal fluctuations, increasing the energy
and temperature until the condition ∆x < λD is fulfilled[45–47]. The phe-
nomenon, referred to as numerical heating, arises from the Finite Grid Insta-
bility (FGI).
This instability can be traced back to the interaction between particles and
the fields on the grid, namely the force excreted depends on the position
of particles relative to the grid [48, Ch. 8.2, pp. 159]. Consequently, the is-
sue results from spatial non-invariance, and it is reasonable to expect that
momentum conservation would mitigate it. In fact, one way to combat FGI
is by using smoother shape functions, since this acts as a low-pass filter
on the spurious small wavelength oscillations [31, Ch. 8, pp. 167][49]. Ex-
act momentum conservation, however, may come at the expense of energy
conservation [48, Ch. 5.5, pp. 149][50]. This appears contradictory, as the
non-conservation of energy (numerical heating) was precisely the issue that
momentum-conserving codes were designed to address. However, there are
multiple ways in which symmetries can be broken, and a discretization that
preserves one conservation law may easily violate another.
To mitigate numerical heating also energy-conserving PIC schemes can be
used, however, at the expense of momentum conservation [51]. The first
energy conserving code was developed by Lewis in the early 70s [52], followed
by many others [53–56].
The problem of finding a scheme that conserves energy and momentum
exactly is, to this day, subject to active research. However, there exist electro-
magnetic PIC algorithms that conserve energy and momentum well, though
not exactly. For many practical purposes, this level of conservation is suffi-
cient. Examples of such schemes are found in [57, 58]. Another example is
the algorithm, later used in this thesis, which conserves energy exactly and
momentum approximately.
NUMERICAL STABILITY 39
4.2.3 Charge conservation
Another quantity which one would like to preserve is charge, or equivalently,
ensure that the continuity equation is satisfied
∂ ρ
+∇· J⃗ = 0. (4.22)
∂ t
As mentioned in section 4.1 charge is generally not preserved by the PIC
algorithm since Gauss law∇·ρ is typically not used in the advancement of
the fields. Historically there has been two main approaches to ensure charge
conservation: (i) using a, so called, charge conserving deposition scheme
for the interpolation of the current [59–63] (ii) correcting the fields retroac-
tively according to eq. (2.9), this strategy is commonly known as divergence
cleaning [31, 64, 65, Chap. 15.6, pp. 359-360].
4.2.4 Numerical dispersion
Numerical dispersion is the spurious dispersion "felt" by the electromag-
netic waves when interacting with grid nodes. It causes light in simulated
vacuum to move in velocities different from the speed of light (generally
slower [66]). Moreover, since it is caused by the interaction with grid nodes,
the dispersion is generally direction dependent. All maxwell solvers based
on Finite-Difference Time-Domain (FDTD) [67] shows numerical dispersion.
Several modifications of the PIC algorithm has been suggested to reduce the
numerical dispersion [68–72]. A practical strategy to eliminate numerical
dispersion involves utilizing a spectral solver, as it fully bypasses the spatial
grid.
4.2.5 Numerical Cherenkov instability
An issue connected to numerical dispersion is numerical Cherenkov insta-
bility (NCI). This is an instability generated when particles are moving faster
than their electromagnetic radiation [40]. Consequently using a spectral
method effectively suppresses NCI, however, in simulations with relativisti-
cally streaming plasmas new modes re-emerge in the intersection of electro-
magnetic modes and numerical aliased modes [39, 40, 73].
To better understand the origin of NCI, we will derive the dispersion relation
for a cold relativistic streaming plasma. The following derivation is presented
40 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
in more detail in [7, 73]. We begin with the expression for the current given in
eq. (2.24). In practice, the discretization of the current and electromagnetic
fields causes aliasing, which can be represented by the substitutions ω→
ω′ =ω+n2π/∆t and kx → k ′x = kx + l 2π/∆x , where n , l ∈Z. The current
and fields are discretized using different schemes, in different steps of the PIC
algorithm. To account for this, we make this substitution in the expression for
the current, while keeping ω and k unchanged for the fields. Consequently,
in the relativistic case, the current takes the form:
∫
i q 3 F⃗ ·∇p f− 0j⃗ (k⃗ ,ω) = d p p⃗ . (4.23)
mγ (ω′− p⃗ · k⃗ ′/mγ)
Assuming the ground state is a cold uniform streaming plasma: f0 = n δ30 (p⃗−
p⃗0) = n0δ(px−p0)δ(py )δ(pz ), the resulting dispersion relation, using eq. (2.23),
becomes
(ω2− c 2k 2)E⃗ +c 2k⃗ (k⃗ · E⃗ ) =
2 ∫ωω F⃗ ·∇ δ3p p (p⃗ − p⃗0)
= d 3p p⃗ . (4.24)
γq ω′− p⃗ · k⃗ ′/mγ
∫
Inserting F⃗ = q (E⃗ +p⃗×B⃗/c γm )and evaluating the integral using d x f (x )δ′(x−
x0) = f ′(x0) results in
(ω2− c 2k 2)Ei + c 2ki k j E j =
 
2 ω ∂ Ek +εk l mεmnp pl kn Ep/γmω =ω p p , (4.25)γ ∂ p i ′k ω −p k ′ /γm n n p=(p0,0,0)
where εi j k is the Levi-Civita symbol, and the magnetic field Bi has been
replaced by Bi =
c
ωεi j k k j Ek given by eq. (2.22). Reducing the system to two
dimensions and assuming kz = 0, leads to a matrix equation εE⃗ = 0, where ε
NUMERICAL STABILITY 41
has the elements
– 2 ™ω
2− 2 2 − pε11 = ω c ky ′ (ωω
′+ v 20 ky k
′
y ) , (4.26)γ(ω′− v 20kx )
– ™
ω2p v0k
′
2 − yε12 = c kx ky (ω− v k ) , (4.27)γ(ω′− v0k ′ 2
0 x
x )
– ™
ω2p v0ky
ε21 = c
2ky kx − ′ , (4.28)γ(ω′− v0kx )
– ™
ω2
ε = ω2− c 2k 2 p22 x − ′ (ω− v k ) . (4.29)γ(ω′− v 0 x0kx )
As in the general case (eq. (2.28)) the system has a linearly independent solu-
tion if det(ε) = 0. By evaluating this equation, one can examine how plasma
waves are coupled. The dispersion relation is given by ε11ε22 − ε12ε21 = 0
and, after some algebra (see appendix A), this equation can be written on
the form
ω2 2€ p ω′ Š€ ω Š
(ω′− v k ′ )2− ω2− c 2k 2 p
ω− v0kx ′ ′
0 x − ′ = g (ω,ω , k⃗ , k⃗ ) (4.30)γ3 ω γ ω′− v0kx
In the non-discretized case ω = ω′, k⃗ = k⃗ ′, v0 = 0 and the first term corre-
sponds to a longitudinal Lorentz transformed plasma oscillation ω− v0k =
ωp/γ
3/2. This can be seen by Lorentz transforming the frequencies as ω→
(ω− pv0kx )γ and ωp →ωp/ γ. While the second term corresponds to the
light-like plasma wave ω2 =ω2 /γ+ c 2k 2p . The coupling term is
2  ′ 2 ωp €ω Š v0c ky
g (ω,ω′, k , k ′) = v 2k k −k ′0 y y y + (k
′ k −k ′y x x ky ) , (4.31)γ ω ω
which reduces to zero when k = k ′ and ω = ω′. However, in discretized
numerical schemes, the discretization introduces aliasing in both frequency
and wavenumber, which causes overlaps between the physical transverse
and plasma oscillation mode, that is when
q
ω2 =ω2p/γ+k
2c 2 and ω′−k ′ 2 3x v =± ωp/γ . (4.32)
Conclusively, the coupling between modes provides a gateway for the system
to exchange longitudinal momentum for transverse when these two modes
overlap. We shall use this information to characterize the NCI in section 4.4.
42 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
4.3 π-PIC: a framework for modular PIC simulations
To investigate the NCI growth, the PIC codeπ-PIC was used. With the motiva-
tion to extend this code for simulating LWFA, two extensions were developed
which is described in this section. Furthermore, the architecture of π-PIC is
described, which enables the testing of different algorithms within the same
framework.
Although, there exists numerous implementation of the PIC algorithm, they
all share a similar underlying structure, outlined in section 4.1. The similarity
of methods provides the foundation for developing a unified framework that
facilitates decentralized development. Within such a framework, benchmark-
ing and the integration of new multiphysics extensions could be performed
in a consistent and collaborative manner. The goal ofπ-PIC is to provide this
foundational framework. This section will briefly describe the architecture
of π-PIC, as well as two extensions useful for simulating LWFA: absorbing
boundaries and moving window.
4.3.1 Architecture of π-PIC
To allow maximum flexibility in modifying and extending π-PIC three in-
teraction levels are available. The first is the Python interface which allows
easy parallelized initialization, modification and retrieval of particles and
fields. The second interaction level are the extensions, which consists of
precompiled fields and particle handlers. Two examples of extensions are
described in section 4.3.2.
The third interaction level is the solver. A solver consists of a PIC solver
which contains a field solver. Although the interface is flexible enough to
accommodate other structures, the PIC solver is intended to handle particle
processing, while the field solver is responsible for advancing the fields. Both
has a standardized interface, meaning that there are a number of mandatory
methods with predefined arguments. These three interaction levels are called
and managed by the ensemble, where the execution order of these methods
is defined.
π-PIC: A FRAMEWORK FOR MODULAR PIC SIMULATIONS 43
Figure 4.4: π-PIC features a modular architecture organized into three inter-
action levels: (i) the Python interface, where simulations are initialized and
advanced, and where particles and fields can be accessed or modified; (ii) ex-
tensions which allows additional functionality, such as moving window; (iii)
the solver where different PIC algorithm can be implemented. A fourth compo-
nent is the ensemble, which executes the advance of the simulation state and
manages particle multi-threading.
4.3.2 Extensions
Two extensions relevant to LWFA are presented in Paper B: the moving win-
dow and absorbing boundaries. The moving window is a widely used tech-
nique in beam–plasma interaction studies that reduces computational cost
by simulating only the region surrounding the driver beam. The remainder
of this section focuses on the implementation of absorbing boundaries, as
this topic is more technically involved.
In general, various boundary conditions are applicable to plasma simula-
44 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
tions, including periodic, reflecting and absorbing. For most electromagnetic
solvers implementing reflecting and periodic boundaries are relatively simple
since they prescribe a predefined value at the boundary. Open-boundaries,
however, are more challenging since they impose a condition on the kind of
electromagnetic waves allowed at the boundary, namely only waves propa-
gating into the boundary, or equivalently, no reflected waves.
For this reason, open boundary conditions are commonly emulated by intro-
ducing an absorbing medium at the boundary. A popular implementation is
Perfectly Matched Layers (PML) [74], where, so called, complex coordinate
stretching is used to emulate an absorbing space within a few grid points.
Here, I take on a minimalistic approach known as masking, where incoming
electromagnetic fields are gradually depleted by iterative attenuation [75].
Figure 4.5: Three simulations of plasma oscillations initialized at the state
shown in (a). The black lines represent the longitudinal electric field. Panels
(b–d) show the plasma state after one plasma period. Simulation (b) employs
absorbing boundaries, simulation (c) uses periodic boundaries, and simulation
(d) features a simulation box four times larger, where the initial state is extended
on both sides by uniform plasma.
A limitation of this method, in contrast to other techniques like perfectly
matched layers (PML), is the necessity for a relatively large damping region.
On the other hand, an advantage of this method is its versatility; it can be
applied to almost any type of Maxwell solver.
To fully emulate open boundaries in plasma simulations particles needs to
be added and removed to maintain a constant density and velocity distri-
bution at the boundary. Similar to fields, this must be conducted gradually
and smoothly to avoid introducing artifacts and high-frequency noise at the
STUDY OF NUMERICAL CHERENKOV INSTABILITY 45
boundary. In π-PIC the removal of particles is done at the same rate as fields,
while simultaneously adding particles with velocity distributions determined
by temperature. To evaluate this method, three simulations with different
boundary conditions were performed, as shown in fig. 4.5: absorbing (b),
periodic (c), and an emulation of an infinite system represented by a simula-
tion box four times larger while at the same resolution (d). The phase-space
structures are similar in the simulations with absorbing boundary and large
simulation box, although slight differences appear in the corresponding elec-
tric fields. The absorbing boundary covered one quarter of the box from each
side, consequently, no significant correspondence is expected within this
region.
4.4 Study of numerical Cherenkov instability
In this section, π-PIC is employed to compare the NCI growth of three PIC
algorithms with different PIC solvers: (i) a Boris pusher combined with a
spectral solver (hereafter referred to as Boris), (ii) an energy-conserving solver
(EC), and (iii) an energy-conserving solver with an additional momentum
correction (EMC). A detailed description of the implementation of the en-
ergy conserving solver can be found at [43] and the implementation of the
momentum correction in Paper B.
4.4.1 Simulation setup
To evaluate simulations in different Lorentz frames it is important to find a
system that has a simple Lorentz transformation. A cold plasma is evidently
such a system. However, introducing a perturbation can accelerate instability
growth, thereby reducing the required simulation time. An electromagnetic
plane-wave traveling through plasma transforms into a new plane wave
under a Lorentz transform and is therefore a good candidate with a simple
evolution. To avoid non-linear displacement in the propagation direction
circularly polarized light is considered, which creates a standing wave. The
46 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
corresponding electric and magnetic field is
A
E⃗ = 0
ω
(−sin(k x −ωt ) ŷ + cos(k x −ωt )ẑ ), (4.33)
c
A0ωB⃗ = (−sin(k x −ωt )ẑ − cos(k x −ωt ) ŷ ), (4.34)
c
where A0 is the amplitude of the vector potential. To initialize the system the
momentum of the particles is needed. By assuming that the electromagnetic
field is much stronger than the field generated by the collective motion of the
particles, as well as E⃗ ≫ v⃗ × B⃗/c in the force term the momentum simplifies
to
d p⃗ q A
= q E⃗ → p⃗ =− 0 (cos(k x −ωt ) ŷ + sin(k x −ωt )ẑ ). (4.35)
d t c
In the simulation density was set to under-dense scales n0 = 1016 cm−3, tem-
perature to small but non-zero: T = 10−12m c 2, A0 = 100 StatV and λ= 1µm.
The choice of these parameters were guided by λ≫λD to mitigate plasma
2
self-interaction and A c m0 < q in order to keep the transverse velocity below
the relativistic limit. The frequency of the oscillations is governed by the
dispersion relationω2 =ω2 +k 2c 2p . The simulation box was set to the length
of one wavelength and the simulation time to one plasma period.
4.4.2 Results and discussion
To examine the dependence of the NCI growth on resolution, simulations
were performed for various resolutions shown in fig. 4.6.(a,b). As seen in
fig. 4.6.(c), the NCI manifests as high-frequency oscillations in the electro-
magnetic field. Notably, the NCI pattern changes its characteristics with
varying resolution. There is a moderate difference in instability growth be-
tween EC and Boris. The growth of NCI differed marginally between the EC
and EMC pushers, the EC solver had on average lower instability amplitude
by less than 1% (this is not shown in fig. 4.6). The instability is more severe
for higher resolutions, which is uncommon for numerical instabilities, but
in this case, likely a consequence of the coupling term eq. (4.31) varying with
different resolutions.
To further characterize the instability, the Ey field was Fourier transformed
in space, shown in see fig. 4.7. The spectrum clearly shows that the insta-
bility predominantly grows at the intersection of the two modes derived in
STUDY OF NUMERICAL CHERENKOV INSTABILITY 47
Figure 4.6: (a-b) shows the instability amplitude of the EC and Boris solver
as Root-Mean-Square (RMS) of the difference between the simulated Ey field
and the analytical solution (Ẽy ). The simulation was performed with γ = 10
for different resolutions and evaluated at the plasma period Tp . The different
resolutions of the simulations corresponding to the black dots and the function
surface is interpolated using a cubic spline. The right figure shows the Ey for
the nine simulations enclosed by the black rectangle in (b).
q
section 4.2.5, namely ω2 =ω2p/γ+ k
2c 2 and ω′ − k ′x v = ± ω
2 /γ3p , where
k ′x = kx +l 2π/∆x andω
′ =ω+m2π/∆t . While there exists instability growth
between intersection of higher order modes, l = [±2,±3,±4], m = 0, the am-
plitude are several orders of magnitude lower than those at l =±1, m = 0.
The instability pattern is similar across all pushers: EC, EMC and Boris solver,
which implies that the instability does not fundamentally violate momentum
or energy conservation. However, the maximum amplitude of the instability
for EC and EMC pusher is ∼ 50% lower than the same value for the Boris
pusher – indicating that energy conservation is beneficial. The analysis also
shows that the Boris pusher exhibits slightly higher growth levels across the
entire k -space, as seen from the elevated background levels in fig. 4.7.
The figure appears to capture to separate effects: (i) the NCI growth and (ii)
numerical heating. When the NCI modes are subtracted, the background
48 PARTICLE-IN-CELL MODELING OF LASER WAKEFIELD ACCELERATION
Figure 4.7: Shows the frequency pattern of the instabilities in k-space for
the solvers energy conserving, momentum and energy-conserving and Boris
pusher. The solution ofω2 =ω2p +k
2c 2 andω′−k ′x v = 0 whereω
′ =ω+n2π/∆t
and k ′x = kx + l 2π/∆x with n = 0, l =±1,±2 is shown in green (l =±1) and gray
(l =±2).
pattern for the EC and EMC solvers shows a gradual decline of numerical
noise with increasing ky , whereas the Boris solver exhibits a more evenly
distributed disturbance. This is likely due to numerical heating which is
prevented by the energy conservation.
4.5 Conclusion and outlook
PIC simulations are essential for the design and analysis of laser wakefield
accelerators. However, these simulations are computationally expensive,
making research into faster algorithms and methods that reduce numerical
inaccuracies essential. An important step in that process is providing frame-
works, which simplifies the comparison and development of algorithms. In
this work, the PIC code π-PIC was used to study the NCI which hinders a
potential speedup of LWFA simulations by evaluating the simulation in a
Lorentz boosted frame.
The origin of the NCI growth was located to the overlap between the light
mode and the Lorentz transformed plasma oscillation. The unphysical cou-
pling between these modes is a consequence of the discretizations scheme
which allows for momentum transfer between longitudinal and transverse
modes.
It was found that, energy conservation does not eliminate the instability, nei-
CONCLUSION AND OUTLOOK 49
ther does energy conservation in conjunction with approximate momentum
conservation. However, conserving energy does reduce the amplitude of the
instabilities, albeit, not to a level that enables boosted frame simulations.
Promising directions for future work include exploring alternative current
deposition schemes. Recent studies [76–78] suggest that the numerical
Cherenkov instability (NCI) can be mitigated by modifying density and cur-
rent deposition schemes. A common assumption in PIC schemes is that a
particle’s density contribution remains constant (or varies linearly) over a
timestep. However, this assumption becomes increasingly inaccurate when
particles move at relativistic velocities relative to the grid, leading to inaccu-
racies in the current deposition.
Additionally, The presence of the coupling term suggests that the violated
conservation property is associated with Lorentz boost symmetry. Due to
discretization, it becomes possible for ω ≠ ω′, a condition that would be
forbidden in a continuous system preserving Lorentz (boost) invariance.
Following this reasoning, the conserved quantity that is broken is not energy
or momentum, but rather the boost current. However, as discussed in the
introduction of this chapter, the correspondence between conserved quanti-
ties and symmetries is not as strong for discretized systems as in continuous
symmetries. Nevertheless, further investigation along these lines would be a
promising continuation of this work, which would deepen the understanding
of the generation and growth of NCI, thereby enabling more computationally
efficient schemes for PIC simulation of LWFA,
5 Research papers
In this final chapter the reader can find the two appended papers:
Paper A
F. BROGREN, S. JALAS, L. HÜBNER, P. MESSNER, M. SCHNEPP, M. TRUNK, C.
WERLE, P. WINKLER, M. MARKLUND, A. GONOSKOV, W. P. LEEMANS, A. R.
MAIER, M. KIRCHEN. "Data-driven modeling of a laser-plasma
accelerator-based X-ray source". Phys. Rev. Accel. Beams 28, 101302 .
DOI: https://doi.org/10.1103/dnnw-lfn5
Paper B
F. BROGREN, J. MAGNUSSON, C. OLOFSSON, A. GONOSKOV. "π-PIC: a
framework for modular particle-in-cell developments and simulations".
arXiv preprint: 10, 48550.
DOI: https://doi.org/10.48550/arXiv.2511.09950
Paper A
F. Brogren, S. Jalas, L. Hübner, P. Messner, M. Schnepp, M. Trunk, C. Werle, P.
Winkler, M. Marklund, A. Gonoskov, W. P. Leemans, A. R. Maier, M. Kirchen.
"Data-driven modeling of a laser-plasma accelerator-based X-ray source".
Phys. Rev. Accel. Beams 28, 101302

PHYSICAL REVIEW ACCELERATORS AND BEAMS 28, 101302 (2025)
Data-driven modeling of a laser-plasma accelerator-based x-ray source
F. Brogren ,1 S. Jalas ,2 L. Hübner ,2 P. Messner,2 M. Schnepp,3 M. Trunk ,2 C. Werle ,2
P. Winkler ,2 M. Marklund ,1 A. Gonoskov ,1 W. P. Leemans,2,3
A. R. Maier ,2 and M. Kirchen 2
1Department of Physics, University of Gothenburg, Origovägen 6B, 412 96 Gothenburg, Sweden
2Deutsches Elektronen-Synchrotron DESY, Notkestraße 85, 22607 Hamburg, Germany
3Department of Physics, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany
(Received 28 March 2025; accepted 26 September 2025; published 21 October 2025)
Laser-plasma accelerators (LPAs) enable compact, bright x-ray sources, but their practical application
demands a significant reduction of beam instabilities that originate from drive laser fluctuations. Moreover,
the complexity of the laser-plasma interaction makes it difficult to disentangle and quantify the impact of
individual parameters on machine performance. To address this challenge, we develop a data-driven
modeling strategy and apply it to an extensive dataset collected during a daylong operation of the LPA-
based x-ray source LUX. By making use of an orthogonal-distance-based training objective, our approach
reduces bias originating from measurement errors, which allows us to estimate the functional dependencies
between laser, accelerated electrons, and generated undulator radiation. In addition, we demonstrate
accurate prediction of the x-ray spectrum based on noninvasive measurements, showcasing the potential of
our approach for building virtual diagnostics.
DOI: 10.1103/dnnw-lfn5
I. INTRODUCTION Carrying out this task with traditional data analysis has
Laser-plasma accelerators (LPA) [1,2] can deliver high- proven difficult due to inherent noise and measurement
energy [3] and high-current [4] electron beams within only uncertainties in these experiments. Yet, recent progress in
a few centimeters. These compact sources are therefore long-term stability and the push toward higher repetition
attractive for applications such as radiation therapy [5] and rates have significantly improved the quality and avail-
medical imaging [6], and, as shown in recent experiments, ability of experimental data [12,13].
capable of driving free-electron lasers down to the vacuum Given these developments, data-driven approaches have
ultraviolet regime [7–9]. Despite these milestones, wide- emerged as a powerful new tool for understanding and
spread adoption of laser-plasma accelerators is still limited controlling the physics of plasma accelerators [14–16]. By
because they lag behind conventional technology in beam leveraging machine learning algorithms, noisy experimen-
quality and shot-to-shot stability. This is primarily a tal data can be converted into high-dimensional models that
consequence of the high sensitivity of the acceleration reveal the underlying system dynamics. Such models not
process to small variations of the drive laser pulse [10,11]. only help in identifying the sources of electron variability
It is widely accepted that resolving this issue primarily but also provide the foundation for autonomous optimiza-
depends on advancements in high-power laser technology. tion or can serve as a virtual diagnostic [17]. Particularly for
However, the complexity of the nonlinear laser-plasma virtual diagnostics, interest is growing as they could enable
interaction makes it difficult to disentangle the influence of applications where maintaining stability is less critical, as
different laser parameters and identify those that have the long as accurate information about the generated electron
beam or secondary radiation is available.
largest impact on the electron beam properties. To enable In recent years, data-driven methods have become
targeted laser development and enhance the LPA perfor- increasingly popular in the field of plasma acceleration.
mance, it is essential to develop a deeper understanding of For instance, several experiments have demonstrated the
the causal relationship between laser, electron, and secon- use of Gaussian process modeling to improve the electron
dary radiation parameters. beam quality via Bayesian optimization [18,19]. Other
studies have shown that it is possible to accurately predict
the electron beam properties from measured laser data,
Published by the American Physical Society under the terms of using different model architectures, such as fully connected
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to neural networks [11], variational neural networks [15], and
the author(s) and the published article’s title, journal citation, random forests [20]. In addition, there are several reports
and DOI. of successful reconstruction of x-ray spectra using
2469-9888=25=28(10)=101302(14) 101302-1 Published by the American Physical Society
F. BROGREN et al. PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
noninvasive measurements, although these studies were parameters and output quantities. Uncertainty in these
conducted at conventional accelerators [21–23]. variables arises from diagnostic inaccuracies and from
So far, these data-driven methods have focused primarily unmeasured fluctuations in the system state, which effec-
on predictive accuracy rather than explainability, in regard tively appear as noise on the observed output.
to estimating the correlations of the underlying physical Using models that assume error-free inputs, e.g., ordinary
processes. Although prediction and estimation are used least squares regression, on data with input errors generally
interchangeably in many fields, there is cause to differ- leads to biased estimates, a phenomenon known as attenu-
entiate between the two. Prediction generally refers to the ation bias [27,28]. Several methods have been developed to
“guessing” of the (future) response of some process, while address this issue, including Deming regression and its
estimation refers to the guessing of model parameters, both extensions [29,30], total least squares [31], instrumental
with the support of data [24,25]. variables [32], and repeated measurements [32,33].
To provide a concrete example, consider a process However, most existing measurement error models rely on
modeled by y ¼ fðx; θÞ, where θ is a model parameter. In specific assumptions about the functional form, such as
prediction tasks, the objective is to minimize the residuals linearity in regression parameters. In contrast, models with
between fðx; θÞ and y by adjusting θ. Conversely, for more relaxed assumptions, e.g., spline-based methods,
parameter estimation, the objective shifts to minimizing require the selection of multiple hyperparameters, which
the error of θ itself [26]. This is often motivated by some is increasingly problematic with higher dimensionality [34].
physicalmeaning of θ that describes the true dependence of y Feedforward neural networks (FNNs) can approximate
on x. Inmany cases, these two objectives align. However, the complex, nonlinear functions in high-dimensional spaces,
differentiation becomes important when there is a measure- making them a natural choice for the type of data-driven
ment error in the input variable x, which is a common inference considered here [35]. Although they make no
scenario in experiments. explicit assumptions about the functional form, they still rely
In this paper, we follow this line of thought and develop a on implicit assumptions about the data-generating process,
versatile approach to data-driven modeling of experiments which can be critical for accuratemodeling. This includes, in
based on feedforward neural networks. Depending on the particular, assumptions about the error distribution.
task, our framework allows targeting either one of the To illustrate this point, we present a synthetic example in
above objectives to derive an explanatory model or pre- Fig. 1, where two neural networks are trained to recover a
dictive model. known underlying function (black) from noisy data. The
Both models are employed to gain specific insights into model trained with mean squared error (MSE) (orange) fails
the physics of a LPA-driven x-ray source. The analysis is to capture the correct relationship, since it is minimizing the
based on a large dataset generated during a daylong
operation of the LUX laser-plasma accelerator [12], where
high-quality electron beams were accelerated to 260 MeV
and subsequently used to generate spontaneous undulator
radiation centered at 12 nm wavelengths.
The paper is structured as follows. In Sec. II, we motivate
our approach and present a detailed derivation of the
modeling process. In Sec. III, we describe the experimental
setup and present our feature selection process, which is
quantified by permutation importance and low correlation
between input variables. In Sec. IV, we demonstrate the
usefulness of the explanatory model as a tool in finding the
causal relation between input and output variables in the
presence of measurement errors. Throughout several exam-
ples, we uncover the influence of individual laser properties
on the LPA process and the generated undulator radiation.
Finally, we shift our attention to the predictive model and
demonstrate its effectiveness in building virtual diagnos- FIG. 1. Two models were trained to estimate the underlying
tics. By using only noninvasive measurements, we recreate dependence ftrue (black) using different optimization criteria:
the x-ray spectrum, achieving an R2 score up to 0.98 for the orthogonal distance (OD) marked in blue and mean squared error
central wavelength. (MSE) marked in orange. The randomly generated data points x
and y had a normally distributed error with μ ¼ 0 and σ ¼ 0.7.
II. DATA-DRIVEN MODELING The black dots represent a few random samples of the validation
set and the 2D-histogram shows the entire dataset. The shaded
As many experimental systems, the LPA process can be area represents the 95% quantile interval across the predictions
characterized by a mapping between measurable input from an ensemble of 50 estimators.
101302-2
DATA-DRIVEN MODELING OF A LASER-PLASMA … PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
FIG. 2. Data is split into a training and a validation set. The training set is randomly sampled with replacement to create subsamples.
Each subsample is given to one network, which is then trained by repeatedly propagating the data, evaluating the orthogonal distance
and updating the weights through back propagation. The final model is calculated from the median of all networks in the ensemble and
the model uncertainty from the quantiles of the distribution of models. Each network (upper right) consists of two FNNs. The first takes
in measurement data x and y and generates an estimation of the x without measurement noise (x). The second FNN takes in x and
generates an estimation of y without measurement noise (y). Once trained fv will be an estimation of the underlying dependence
between x and y.
vertical distance between model and data points. In cases where xi are the true measurement values and εx;i; εy;i are
where there are errors on both input (x) and output (y), the assumed to be normally distributed (multivariate) noise
optimization criterion should instead minimize the orthogo- parameters with the same dimensionality as xi; yi. The goal
nal (or total) distance (OD) between model and data, as is to estimate f. We model the system using two FNNs
demonstrated by the improved reconstruction obtained with fw; fv with weights w; v, where fwðxi; yiÞ is an estimation
OD minimization (blue) in Fig. 1. of xi and fv(fwðxi; yiÞ) is an estimation of fðxi Þ. When
In the following Sec. II A, we present the implementa- optimized, the function fv will be an estimation of f. A
tion of such an OD optimization criteria in the context of schematic figure of the model is shown in Fig. 2 (upper
neural networks. A similar approach as the one developed right). To find the values of v; w such that fv; fw align with
in this paper can be found in [34]. Note that the goal of this the data, the likelihood is maximized, that is, maximizing
model is to estimate the dependencies between input and the probability of obtaining the data x ¼ fxig; y ¼ fyig
output of the system, not solely to predict the system given v; w: Pðx; yjv; wÞ. The presence of measurement
output. This is an important distinction since it affects the errors on output and input implies that for optimally trained
optimization condition and consequently the modeling. In networks the data Di ¼ ½xi; yi should be distributed
Sec. II B, we turn to the problem of prediction and elaborate according to a multivariate Gaussian distribution centered
on the distinction between prediction and estimation. We around μwv;i¼½fwðxi;yiÞ;fvðfwðxi;yiÞÞT . This assumption
call the models resulting from these two approaches the results in the likelihood
explanatory model and predictive model.
A. Explanatory model Pð 1x; yjv; wÞ ¼ Π pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e−1½Di−μ T −12 wv;i Σ ½Di−μwv;ii ;k
For formulating a mathematical description of the ð2πÞ detðΣÞ
system, we define xi and yi as the vectors of measured (2)
input and output variables, where i∈ ½1; N and N are the
number of measurements. Note that x ¼ ½x ;…; x Ti i1 in andx where Σ is the covariance matrix of the errors εx and εy with
yi ¼ ½yi1;…; y Tin are vectors where nx; ny are the numbery dimension k × k, where k ¼ nx þ ny. Assuming that Σ is
of input variables and output variables, respectively. diagonal, that is, assuming no covariance between errors,
The problem of measurement errors can mathematically the expression is simplified. Moreover, by taking the
be expressed as negative log of the expression and discarding all constant
¼  þ ¼ ð Þ þ terms a more computer-friendly optimization condition forxi xi εx;i; yi f xi εy;i; (1) v; w is obtained
101302-3
F. BROGREN et al. PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
X ð − Þ X ½ − ð ð ÞÞ 2 example that the OD loss outperforms the MSE loss exceptD μ 2 yij0 fv fw xi; yi 0 2 2
min ij wv;ij ¼ min j in the limiting case where σx;j ≪ σ
2 2 y;j
0 .
v;w v;w
i;j σj σi;j;j0 y;j0 Under optimal training, fw is expected to approximate an
− inverse of fv, in the sense that f
0 0 0 0
½ ð Þ  w
ðx ; fvðx ÞÞ ¼ x where x
2
þ xij fw xi; yi j is an arbitrary vector. However, the underlying function (f); (3)
σ2x;j may not be bijective, leading to ambiguity when the
gradient is close to zero, as for x ¼ 0 in Fig. 1.
where the index j and j0 are running over all input and We found that this flexibility can result in sharp changes
2 2 of fw and large variation of function shapes betweenoutput variables, respectively, and σx;j; σy;j0 are the diagonal different training sessions. A common reason for high
elements of Σ, that is the variance of the errors εx;j; εy;j0 . variance are large neural weight values; a widely utilized
Consequently, maximizing the likelihood is equivalent to method to address this issue is weight regularization. In L2-
minimizing the orthogonal distance between the model and regularization the weights are pushPed toward lower valuesthe data. We refer to the term in Eq. (3), prior to by introducing a penalizing term η 2i jwij in the loss. In a
minimization, as the OD loss. Note two things about the Bayesian framework, it can be shown that L2-regulariza-
above expression: (i) The variances σ2 ; σ2x;j y;j0 , which in the tion is equivalent to imposing a Gaussian prior on the
experiment corresponds to the variance of the measurement weights [37]. Note that η is a regularization parameter and
error, are free variables and need to be estimated. (ii) If xi is must be chosen based on convergence and model perfor-
measured with a small error compared to y 2i, i.e., σx;j ≪ mance. The training loop can be summarized as (i) One
σ2 0 for all j; j0, it follows that f ðx ; y Þ ≈ x since the last batch of data is propagated through the networks, whichy;j w i i j ij
term in Eq. (3) would drastically increase otherwise. gives fwðxi; yiÞ; fvðfwðxi; yiÞÞ. (ii) The OD loss L ¼
Consequently, the optimization criteria effectively trans- P ðD 2 P Pij−μwv;ijÞ − η jw j2 − η 2form into the MSE. i;j 2 w i i v i jvij is calculated.σj
As a result, the OD loss includes the MSE optimization (iii) The partial derivative with respect to w; v of L is
criteria. For MSE, however, assumptions regarding error calculated. (iv) The network is updated by wþ ¼ s ∂Lw ∂w þ w,
distributions remain implicit, whereas they are explicitly vþ ¼ s ∂Lv ∂v þ v, where sv; sw are the step lengths.stated in the OD loss. From a probabilistic view the training data is sampled
Ideally, the measurement error for each input and output from a larger distribution of data, therefore a model
variable should be estimated individually and incorporated optimized on a subsample might differ from a model
into the OD loss. However, obtaining accurate error created on another subsample. This uncertainty is called
estimates is often challenging in practice. In the context model uncertainty or epistemic uncertainty and will
of experimental data, assuming equal measurement uncer- decrease with the amount of data [38,39]. Model uncer-
tainty for both input and output variables (i.e., setting
¼ 0 tainty can be approximated by creating subsamples of theσx;j σy;j0 for all j; j ) is arguably a more reasonable dataset using random sampling (with replacement) and
approximation than the one implicit in MSE minimization, training models separately on these subsets, see Fig. 2. This
which presumes σ2x;j ≪ σ2y;j0 . technique is known as bootstrapping [40]. The variation of
It is important to note that only the ratio between σx and these models approximates the model uncertainty [41],
σy affects the model outcome, not their absolute values. furthermore this uncertainty will grow when the data is
This follows directly from the structure of the OD loss sparse, which makes the modeled uncertainty a good
Eq. (3), where uniform scaling of all terms leaves the indication of the region where the model is trustworthy.
position of the minimum unchanged, whereas changes in
their relative weights do not. B. Predictive model
Based on this reasoning, we assume that the variance of
the measurement error is on the same scale as the root- The previous section outlined a method for estimating
mean-square deviation (RMSD) of the data. Since neural the function describing a process, done by maximizing the
networks benefit from inputs and outputs being similarly likelihood Pðx; yjv; wÞ. Shifting focus to prediction of y
scaled, the data are normalized to zero mean and unit given x, the likelihood to maximize is Pðyjx; uÞ, where u
RMSD prior to training. Under this normalization, the are the weights of the model fu. Furthermore, y is no longer
assumption of equal measurement errors is implemented by allowed as an input to the model, since in this case it would
setting σ2 ¼ σ2 ¼ 1 for all j; j0. not be predicQtive. Therefore, the likelihood is expressed asx;j y;j0
Further details on the choice of measurement error Pðyjx; uÞ ∝ exp½
jyi−fuðxiÞj
i 2  and the optimization criterion2σy
parameters are provided in Supplemental Material [36], reduces to the MSE. Conclusively, for prediction, meas-
where we also discuss the consequences of mismatched urement errors of input variables are ignored and the
error assumptions. In particular, we show in an illustrative predictive model consists of only one FNN.
101302-4
DATA-DRIVEN MODELING OF A LASER-PLASMA … PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
III. EXPERIMENT AND DATA resolution) that analyzed the scattering signal from the
A. Experimental setup pyroelectric sensor. The beam profile and wave front
properties of the laser were analyzed using the leakage
The experimental setup, shown in Fig. 3, consisted of the through the off-axis parabolic mirror. The measured wave
LUX laser-plasma accelerator with a subsequent electron front data were decomposed into Zernike coefficients for
beamline to generate undulator x-rays. The laser-plasma detailed characterization. The focus position of the laser
interaction was driven by the Ti:sapphire laser system inside the plasma source was inferred from the defocus
Angus, operating at 1 Hz and delivering around 2.6 J pulses value of this measurement. The wave front sensor was
with a duration of 34 fs. The laser pulses were focused by a calibrated relative to the focus position inside the plasma
f=25 off-axis parabolic mirror (OAP) to a spot size of source by shifting a lens to introduce defocus and observ-
25 μm FWHM in the plasma source. ing the plasma recombination light at the focus. The
The laser energy was monitored by a pyroelectric sensor repeatability of the wave front sensor (< 4 nm) corre-
(< 0.5% repeatability) that measured the rejected beam sponds to a focus position uncertainty < 100 μm.
from a polarizer inside a tunable attenuator at the end of the The plasma source consisted of a sapphire plate with a
laser chain. The spectral properties of the laser were capillarylike structure and three gas inlets, providing a
measured using a fiber-coupled spectrometer (1 nm region with a mixture of nitrogen and hydrogen in the front
0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000
FIG. 3. (a) Shows the experimental setup and (b) the shot-by-shot x-ray and electron spectra over the 24 h long experiment. The laser
energy and spectrum were measured before temporarily compressing the laser pulse. Beam profile and wave front were measured from
the leakage of the OAP. The OAP focused the laser into the plasma source, where electrons were accelerated using downramp-assisted
ionization injection, the electrons were thereafter focused by two quadrupoles into the undulator. The electron spectrum was measured
through the deflection of a dipole magnet and the x-rays were focused onto a transmission grating, where the resulting spectrum was
captured by a x-ray CCD. Charge and beam position were measured before and after the undulator with beam position monitor (BPM) 1
and 2. The input variables of the explanatory model were measurements of energy, laser spectrum, beam profile, as well as wave front,
and output variables were statistics from the electron and x-ray spectrum, as well as the charge measured by BPM1.
101302-5
xray E (MeV)
F. BROGREN et al. PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
for ionization injection, followed by pure hydrogen region B. Feature and model parameter selection
for acceleration. The design enabled electron generation Strong correlations between input variables, a phenome-
from downramp-assisted ionization injection [11]. This non known as multicollinearity, can obscure the individual
process involves ionizing the outer shell electrons of contribution of features and reduce model interpretability.
nitrogen, which creates a density peak within the plasma To address this, we used the variance inflation factor (VIF)
at the front of the plasma source. After this density peak, and permutation feature importance [42] to select a subset
the resulting density transition temporarily reduces the of informative and independent features.
wakefield phase velocity, allowing electrons to be locally The VIF is calculated as VIF ¼ 1=ð1 − R2Þ, where R2i i is
injected into the plasma wave. These electrons are then the R2 score of a linear regression model that predicts the ith
accelerated in the density plateau region formed by the pure input variable based on all other input variables. Permutation
hydrogen present in the latter part of the plasma source. feature importance quantifies a variable’s importance by the
This setup delivered on average electrons with median decrease in model accuracy when that variable is shuffled.
energy E ¼ 260 MeV (11 MeV RMSD), median absolute For this, any model with enough expressibility can be used;
deviation ΔE ¼ 5.5 MeV (2.6 MeV RMSD), and charge here we used feedforward neural networks (FNNs). Each
Q ¼ 60 pC (16 pC RMSD). FNN consisted of two hidden layers of 32 neurons each. The
Electron beams exiting the plasma were focused into the accuracy of the model was quantified by the root-mean-
center of the BEAST II undulator by a quadrupole doublet. squared error (RMSE) for the output variables: charge (Q),
The undulator, consisting of 100 periods with 5 mm period electron median energy (E), and electron energy median
length, a fixed 2 mm gap, and undulator parameter deviation (ΔE).
K ¼ 0.29, generated radiation from the accelerated electron By demanding VIF below 5 and high importance the set
beams. Before and after the undulator, the charge and of potential laser diagnostics and statistics used as input
centroid of the electron beam were measured using beam variables could be reduced from 70 to 11. For physical
position monitors (BPMs). Following the undulator, the interpretability, wave front aberrations with order 9 or
electron beams were directed into a dipole spectrometer, higher were excluded from this set. For completeness, all
which achieved an energy resolution of 0.1% within a pairs of variables (e.g., beam profile x- and y-positions)
20 MeV range around the central energy. were included, even in cases where only one exhibited high
The spontaneous undulator radiation was focused using importance. This yielded a final input variable set of 14,
a grazing incidence toroidal mirror, which imaged the shown in Fig. 4. In contrast, incorporating all available
center of the undulator into a radiation spectrometer. This diagnostics only increased the R2 score by 3.4%, sug-
spectrometer consisted of a free-standing gold grating that gesting that the additional variables offered little new
dispersed the radiation by wavelength, with the resulting information to the model. The laser focus position was
spectrum measured on an x-ray CCD. The x-ray spectrum the most critical input variable, causing the RMSE to
measured on average had a central wavelength of 12 nm increase by up to 200% when shuffled, demonstrating its
(0.8 nm RMSD) with bandwidth 1.9 nm (0.18 nm RMSD). substantial influence on model accuracy. Unexpectedly, the
These diagnostics allowed the simultaneous measurement method identified the peak signal of the laser beam profile
of the electron and x-ray properties for all individual shots. as an important feature with a low VIF, suggesting that it
The accelerator was operated continuously for more than captures unique information although no direct physical
24 h, producing more than 80,000 shots. For the analysis, interpretation could be identified.
the first 40,000 shots were used, which allowed us to The output variables of the explanatory model were
designate a single training set (shots 0–20,000) and a chosen as Q, E, ΔE as well as x-ray mean wavelength
subsequent validation set (shots 20,000–40,000). The split (λxray) and standard deviation (Δλxray).
was performed chronologically to ensure that the training For the virtual diagnostic, models were trained using
period preceded the validation period, thereby avoiding different combinations of three input sets: the laser param-
information leakage between the datasets that could com- eters shown in Fig. 4, the electron bunch position and
promise the integrity of the validation performance. Using charge measured by the beam position monitors, and the
the full dataset would have required multiple training and electron energy spectrum, which was reduced to 20
validation intervals to address long-term laser drifts that principal components using PCA. The models were trained
necessitate periodic model retraining—a strategy that, to predict the x-ray spectra, which were likewise decom-
while feasible, was set aside for simplicity. posed via PCA and reconstructed during inference using
Measurements of distributions and profiles were sum- the inverse transform.
marized in statistics, which, together with scalar diagnos- In addition to selecting an appropriate set of input
tics, resulted in 70 defining laser parameters per shot. This variables, there are also a number of model parameters
extensive dataset allowed for robust analysis of the sys- affecting performance, such as number of neurons, number
tem’s performance and stability, providing a comprehen- of hidden layers, weighting, and learning rate. A rigorous
sive foundation for the results presented in the paper. optimization of these parameters was outside the scope of
101302-6
DATA-DRIVEN MODELING OF A LASER-PLASMA … PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
FIG. 4. Feature importance for predicting the median energy (E), median energy absolute deviation (ΔE), and charge (Q) of the
electron bunch, shown as a percentage increase of RMSE in log scale. Multicollinearity is quantified by the variance inflation factor
(VIF). The bar uncertainty marks the 95% quantile of 50 models.
this work; however, the relevant parameters were varied one inevitably introduce small systematic offsets, a residual
by one and the lossmonitored to find parametersmaximizing mismatch between experimental inputs and simulated ones
performance. The explanatory model consists of networks is expected. To enable a fair comparison with the explana-
featuring two FNNs with two hidden layers containing 32 tory model, the discrete simulation results were interpolated
neurons each and Tanh activation. For training, an Adam using cubic splines to create a smooth response surface,
optimizerwith aweight decay ofη ¼ 10−5 anda learning rate which was then translated in input space by applying
of 10−3 was used. The model was trained for 600 epochs. constant offsets in laser focus position and energy. These
Details on the optimization of theweighting parameter can be offsets were optimized by minimizing the difference
found in Supplemental Material [36]. The predictive model between function surfaces resulting from simulations and
consists of one neural network with two hidden layers explanatory model. No rescaling was applied, preserving
containing 64 neurons each with Tanh activation. No weight the relative spacing of simulation points and the physical
decay was used since the convergence was already satisfac- trends predicted by the simulations.
tory. This model was trained for 100 epochs.
The explanatory and predictive model are ensemble IV. RESULTS
models, that is, they consist of 50 networks with identical A. Influence of focus position and pulse energy on
architecture, each trained on a different dataset created by electron injection and acceleration
random sampling of the training data. The final output is
given by the median of all networks, see Fig. 2. The models To gain insight into the mechanisms governing the
were built using the Python package Ensemble-PyTorch [43], electron and subsequent x-ray generation, it is essential
modified to accommodate the OD optimization criteria. to extract the underlying relationships between laser and
resulting beam properties. However, the laser parameters
were not scanned systematically but varied simultaneously
C. Particle-in-cell simulations due to the natural shot-to-shot fluctuations in the experi-
To compare our approach with conventional LPA mod- ment. Consequently, their individual influence is obscured
eling, we performed numerical simulations using the when simply looking at marginalized representations of the
quasicylindrical particle-in-cell (PIC) code FBPIC [44]. data and any observable correlations are further reduced by
In these simulations, the laser focus position and energy noise. To address this, we use the explanatory model that
were randomly sampled over ranges of 1.4 J and 2.7 mm, disentangles these coupled effects and helps reveal the
respectively. The plasma profile was derived from compu- underlying correlations otherwise hidden by the multidi-
tational fluid dynamics simulations. The laser was modeled mensional parameter space.
with a flattened Gaussian beam profile [45] and a Gaussian The effectiveness of this approach is illustrated in
temporal profile. Figs. 5(a)–5(c), where the training data is shown as a 2D-
The simulations were conducted in a boosted Lorentz histogram and compared with the dependencies of charge,
frame with γboost ¼ 5 to accelerate computation. A longi- energy, and energy spread on focus position extracted from
tudinal resolution of 50 cells per micrometer and a trans- the explanatory model.
verse resolution of 3 cells per micrometer were used, with The model shows that shifting the focus position down-
three azimuthal modes to capture the essential physics. stream, i.e., to larger values, reduces the beam charge. This
Because uncertainties during the laser characterization is a result of the drop in laser intensity at the upstream
and the idealized conditions imposed in the PIC setup region of the plasma source, where the mixed gas is located.
101302-7
F. BROGREN et al. PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
foc foc foc foc
FIG. 5. Left column (a)–(c) shows the explanatory model and PIC simulations as a function of focus position for charge Q, electron
energy E, and energy spread ΔE, as well as the distribution of the training data. The middle two columns (d)–(i) show a comparison
between modeled surface (left) and PIC simulated (right) in a matched region. The PIC simulated surface (g)–(i) has been extracted from
a larger range of simulations, shown in (j)–(l). The simulations, marked by black dots, generate a function surface through a cubic spline.
The red rectangle is the region shown in (g)–(i) and is extracted by subtracting the model surface of (d),(f), corresponding to Q and ΔE,
while moving over the PIC simulated surface in (j), (l) to find the point of maximal correspondence.
The lower intensity reduces the ionization rate and thereby observation is in good agreement with PIC simulations
decreases the injected charge. The opposite trend is found within the same experimental parameter range, see
for the electron beam energy, which increases when the Figs. 5(g)–5(i).
focus is shifted downstream. This effect can be attributed to To directly quantify the impact of each of these laser
beam loading [46]. In general, the electron beam absorbs parameters on beam stability, we used the explanatory
energy from the wakefield and thereby suppresses the model to relate shot-to-shot fluctuations of focus position
accelerating field. As fewer electrons are injected (due to and laser energy to electron energy. For instance, a 1 MeV
lower intensity at the downstream focus), the wakefield RMSD fluctuation in electron energy corresponds to a
remains strong and the beam gains more energy. In 7 μm RMSD shot-to-shot variation in focus position. To
contrast, a higher injected charge would extract more cause the same 1 MeV RMSD variation would require a
energy from the wake, weaken its accelerating field and 0.03 J fluctuation in laser energy. In our experiment, 7 μm
reducing the final beam energy. Finally, we observe a V- represents only 4% of the measured RMSD focus-position
shaped dependence of the energy spread on the focus shot-to-shot variation, whereas 0.03 J is about 160% of the
position, also caused by beam loading. As the electron RMSD laser-energy shot-to-shot variation. Consequently,
bunch depletes the wakefield along its length, the accel- within the typical operational stability of today’s laser
erating field weakens toward the tail. Near zfoc ¼ 0, this systems, the electron beam properties remain far more
depletion balances the natural wakefield gradient so that all sensitive to changes in focus position than to variations in
longitudinal slices experience similar accelerating fields, laser energy. This finding stands in contrast to the usual
yielding a minimal energy spread. Moving away from this emphasis on laser-energy stability in discussions of next-
optimum alters the field experienced by the bunch tail, generation systems, suggesting that wave front fluctuations,
increasing the energy spread and giving rise to the V-shaped such as focus-position shot-to-shot variation, may be the
dependence. more pressing issue for improving beamquality and stability.
In addition to focus position, the model lets us Besides validating the explanatory model, the PIC
examine how laser energy affects the beam. As shown simulations allow exploration of parameter regimes not
in Figs. 5(d)–5(f), electron beam parameters remain largely covered by our measurements. The broader range of laser
unaffected by variations in laser energy, in stark contrast to energy and focus position available through PIC simula-
their pronounced dependence on focus position. This tions [Figs. 5(j)–5(l)] reveals that laser energy becomes
101302-8
DATA-DRIVEN MODELING OF A LASER-PLASMA … PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
more influential once varied on a larger scale. As expected,
the beam charge increases with higher laser energy and
when the focus is closer to the mixed-gas region.
Meanwhile, electron energy and energy spread remain
primarily dictated by focus position. A plausible explan-
ation is that while focus position predominantly sets the
injected charge, laser energy simultaneously affects both
the injected charge and the wakefield strength. As a result, a
stronger wakefield can sustain a higher charge without
overloading, thus leaving the beam energy and energy
spread largely unchanged.
Note that the explanatory model and simulations have a
maximal overlap at a point where the difference between
modeled and simulated laser energy is as low as 5%. This
indicates that the model, despite the low significance of
laser energy indicated by feature importance, is able to
extract a likely dependence between electron properties and
laser energy. The overall good agreement between simu- FIG. 6. (a) Shows the model surface for median energy spread
lation and model, indicates the future possibility of multi- (ΔE) as a function of the central laser wavelength deviation λ − λ̄,
fidelity modeling, where simulated data is integrated with where λ̄ is the mean central wavelength over all data samples.
(b)–(f) show ΔE as a function of the wave front aberrations:
experimental models to identify optimal points beyond the vertical astigmatism, horizontal coma, primary spherical wave
range of experimental data. front as well as focus position (zfoc). The orange dashed line
shows the chromaticity of the wave front measurement that is the
B. Chromatic effects and higher-order total chromaticity of the optics in the beam path from the OAP to
the wave front diagnostics.
wave front aberrations
Additional laser variables of importance, see Fig. 4,
included the spectral properties of the laser pulse as well as Fig. 3). The induced measurement error Δzfoc;err ¼
higher-order wave front aberrations. Variations in the cdiagðλ − λ̄Þ effectively leads to equipotential lines on the
spectral properties have a direct impact on the LPA process. model surface with slope equal to cdiag, as can be seen in
For example, both the laser vector potential and the pulse Fig. 6(a). For a derivation, we refer to Supplemental
duration—which are critical for plasma wave excitation— Material [36].
depend on the central wavelength and spectral bandwidth. We estimate that approximately 30% of the total variation
In our measurements, the RMSD of these parameters was in the measured focus position can be attributed to chromatic
0.3 and 0.1 nm, respectively. While the variations in effects. In addition, the observation cdiag ≈ 2csys indicates
bandwidth could result in a percent-level change in laser that measurement errors constitute a significant portion of
intensity (via the pulse duration), the direct effects alone this variance, rather than reflecting true variations of the
cannot fully explain the influence of central wavelength focus position in the plasma source. Note that one could, in
shifts as captured by the model. principle, remove this measurement error of the focus
Nevertheless, even such small spectral changes can have a position by subtracting c ðλ − λ̄Þ from z , assuming that
substantial indirect influence on the LPA process through diag focthe spectrum is accurately measured. Nevertheless, a much
chromatic effects. A prime example are chromatic aberra- better approach would be to decrease the chromaticity of the
tions introduced by transmissive optical elements in thebeam diagnostics in future experiments. In summary, the analysis
path, which effectively cause a shift in focus position,
¼ ð − Þ shows that even though variations in the central laserΔzfoc csys λ λ̄ , given a change in central wavelength wavelength have a low direct effect on the acceleration
λ relative to the average λ̄. For our setup, we esti- dynamics, they can cause indirect chromatic effects and
mate csys ≈ 100 μm=nm. introduce measurement errors.
This focus shift should be accounted for by the wave The most significant parameters after focus position and
front sensor, as it measures the average focus position of all central laser wavelength were the higher-order wave front
wavelengths. However, the dependence between laser aberrations: vertical astigmatism, horizontal coma, and
wavelength, focus position, and electron properties indi- primary spherical, see Fig. 4. Figure 6 shows that the
cates the presence of another correlation between those electron properties have a similar relationship with these
properties. This can be attributed to additional chromaticity parameters as with focus position, namely, along these
cdiag ≈ 215 10 μm=nm that is accumulated in the diag- dimensions there is a minimum of energy spread. This
nostic beam path toward the wave front sensor (compare suggests that the higher wave front aberrations mainly affect
101302-9
F. BROGREN et al. PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
the electron spectra through the intensity in the injection xray xray
region, hence through the process of beam loading. The
foc
remaining wave front aberrations, oblique astigmatism and
vertical coma, had low variance compared to vertical
astigmatism and horizontal coma, which explains their
low importance, see Fig. 4.
To gain more knowledge about the influence of wave foc
front aberrations on the LPA, it is necessary to develop
datasets with larger variance in the wave front and analyze
alongside more realistic 3D PIC simulations than presented
in this work. Nevertheless, already at this level of accuracy,
it would be possible to use the model surface, partly shown
in Fig. 6, for tuning the electron quality by adjusting the
wave front.
C. Laser influence on x-ray generation xray xray max
In the preceding sections, we concentrated on the laser- FIG. 7. (a) Shows a waterfall plot of x-ray spectra, consisting of
plasma interaction and its immediate impact on the quality 40,000 shots, ordered and binned by electron median energy, and
and stability of the electron beam. Now, we shift our focus normalized by themean intensity of the peak signal of the 0th order
to the implications these dynamics can have on the diffraction. Black line marks the analytic estimation for the x-ray
generation of secondary radiation. wavelength with emission angle θ ¼ 0. (b) Shows the spectra in
To build a consistent picture, we begin with a straight- (a) at three different electron energies marked by red, black, and
forward analysis of the raw experimental data. Figure 7(a) blue dashed lines in (a). Subplot (c) and (d) show the distributions
shows the integrated x-ray spectrum sorted by the median of median energy and electron energy spread to central wavelength
and x-ray bandwidth. The explanatorymodel in (a), (c), (d) (dashed
electron energy (a moving mean has been applied along this green) is extracted by scanning focus position, zfoc, while fixing alldimension to reduce noise). The spectra exhibit a relatively other input to theworking point. The focus position is ranging from
broad bandwidth, with a clear dependence of the central zfoc ¼ ½−330; 330 ≈ ½−2σz; 2σz μm. Note the difference in
wavelength λxray on electron energy E. We can compare this strength of correlation between (c) and (d), which is a consequence
observation with the theoretical relationship between elec- of the underlying dependence on zfoc.
tron energy and x-ray wavelength, which is described by
the undulator equation
  deviation of the measured x-ray spectrum as output
ð Þ ¼ λu þ K
2
λ 2 2 parameters.xray γ; θ 1 þ γ θ ; (4)
2γ2 2 The spectral shape of the x-ray radiation is further linked
to the energy spread of the electron beam. Under optimal
where λxray is the radiation wavelength, λu the undulator beam loading conditions, where the energy spread is
period, K the undulator parameter, γ the relativistic factor, minimized (ΔE ≈ 2 MeV), broadening effects associated
and θ the observation angle. The LUX beamline has K ¼ with variations in γ are negligible. However, when the
0.29 and λu ¼ 5 mm. This relation is shown as a black solid plasma wave is nonideally loaded, deviations in the beam
line in Fig. 7(a), assuming on-axis emission (θ ¼ 0). energy occur, accompanied by an increased energy spread.
It is important to note that this theoretical curve intersects This leads to additional spectral broadening that is par-
the slope of the measured spectra toward shorter wave- ticularly visible toward shorter wavelengths.
lengths, rather than coinciding with their peak. This The behavior is further analyzed in Fig. 7(b), which
discrepancy is explained by the combined effects of off- compares normalized and shifted spectra for three beam
axis observation angles, captured within the spectrometer’s loading conditions: underloaded (red), overloaded (blue),
angular acceptance (∼6 mrad), and the natural divergence and optimally loaded (black). As the energy spread
of the electron beam as it propagates through the undulator. decreases, the lower end of the spectrum steepens, revealing
Both factors cause a broadening of the radiated spectrum a direct correlation between electron beam quality and x-ray
toward longer wavelengths due to the quadratic dependence characteristics.
on θ in the governing Eq. (4). Additionally, transverse The power of our explanatory model lies in its ability to
variations in the undulator parameter K can lead to further directly link these observed x-ray spectral effects to specific
broadening of the spectrum given the finite size of the laser properties. In order to isolate the contribution of
electron beam. Note that the predictions of the explanatory individual laser parameters, we can systematically vary
model [green dashed line in Fig. 7(a)] intrinsically account them independently and analyze the corresponding model
for these effects, since we use the mean and standard output. Through this approach, we identified that the laser
101302-10
xray xray
DATA-DRIVEN MODELING OF A LASER-PLASMA … PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
focus position is particularly decisive in shaping the
spectral characteristics of the emitted x-rays.
In Figs. 7(c) and 7(d), we present the measured energy
and relative bandwidth of the x-ray and electron spectrum,
with the explanatory model as a function of the focus
position, zfoc. Remarkably, scanning the focus over a range
of 660 μm, reproduces the strong correlation between the
electron median energy E and the central x-ray wavelength
λxray. As discussed previously, the shift in wavelength
originates from changes in the electron energy, which
are themselves induced by variations in beam loading.
The same physical principles are evident when comparing
the energy spread of the electron beam with the x-ray
bandwidth. As the focus is moved toward optimal charge
injection, both the electron energy spread and the x-ray
bandwidth decrease. However, moving the focus beyond the
optimal point reverses this trend, causing both quantities to
rise again. Furthermore, there is a difference between x-rays
generated by laser pulses withmore upstream or downstream
focus position, despite having the same relative energy
bandwidth. This could be explained by the corresponding
change in electron energy, which alters the beam transport.
While this effect is hidden in the rawdata, themodel is able to
elucidate the process by extracting the relationship between
the focus position and the properties of electrons and x-rays.
Beyond enhancing our understanding of the underlying
physics, these insights have practical implications. For
FIG. 8. (a) Shows the prediction for three models with different
instance, the parametrized function could serve as a tuning input variables (M1, M2, M3) and the measured spectrum (black
curve for the x-ray source. Simply adjusting the position of dashed). M1 uses the laser properties shown in Fig. 4, M2 uses
a lens within the final beam expander telescope of the laser the same laser properties and data from two beam position
system would offer a straightforward method to control the monitors, and M3 uses all aforementioned properties and the
x-ray output. With further refinement and by incorporating measured electron spectra. (b) Shows 50 consecutive measured x-
additional parameters, such a parametrized model could ray spectra and (c) the corresponding prediction of model M3.
even enable independent control over multiple beam (d) Shows the residual between measured and predicted spectrum
properties. In the long term, the approach may also form of model M3.
the basis for an active feedback system to maintain stable
x-ray output, ensuring consistent performance in future
applications. providing precise data without interfering with the exper-
imental setup.
For this predictive task, we compare three models, each
D. Virtual x-ray diagnostics using a different set of input variables. Model M1 uses only
The previous sections emphasized the usefulness of the laser properties, consistent with the input parameters used
explanatory model to gain insight into the physics. in the explanatory model. Model M2 adds noninvasive
However, in x-ray applications, such as time-resolved electron diagnostics, specifically data from beam position
pump-probe experiments [47], photoelectron spectroscopy monitors (BPM), which capture the bunch position and
[48], and absorption spectroscopy [49], the primary objec- charge. Importantly, this model remains fully noninvasive,
tive would not be to understand the underlying processes, as it leaves both the electron beam and the x-ray beam
but simply to obtain accurate knowledge of the radiated intact. This makes M2 potentially relevant for pump-probe
spectrum. This becomes particularly important when the experiments where both beams are used. Finally, model M3
source exhibits high variance, as discussed in [22]. incorporates measurements from the destructive electron
A crucial challenge in these applications is that standard spectrometer. This provides additional detailed information
diagnostics for characterizing x-ray sources are often about the electron spectrum. M3 is suitable for scenarios
destructive, making them unsuitable for real-time use. where only the x-ray beam is of interest.
This is where virtual diagnostics could play a transforma- Figure 8(a) illustrates the predicted spectrum for a
tive role. By leveraging predictive models, they can provide selected shot using each model. While a virtual diagnostic
noninvasive and real-time estimates of the x-ray spectrum, based solely on laser data alone already enables a decent
101302-11
F. BROGREN et al. PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
prediction of the x-ray spectrum (M1), it is outperformed to address this gap by guiding the development of more
by models that include direct information about the electron stable LPA systems and enabling model-based control and
beam (M2 and M3), as confirmed by the R2 scores for feedback strategies.
predicting the central wavelength, 0.61 (M1), 0.93 (M2), Looking ahead, the increase in high-performance com-
and 0.98 (M3). puting resources and experimental data—fueled by the
Remarkably, merely adding the measurements of the upcoming era of high-repetition-rate lasers—will only
beam charge and position to the input in M2 increases the increase the relevance of data-driven modeling in the field.
prediction accuracy by over 50% compared to M1. This Leveraging these developments will be critical to trans-
gain reflects the fact that these measurements provide direct forming laser-plasma accelerators from proof-of-concept
information on key beam properties: the charge encodes systems into practical, application-ready machines.
aspects of the energy spectrum via beam loading, while the
position reflects the trajectory through the beamline, which ACKNOWLEDGMENTS
not only determines how the beam interacts with the
undulator field but also carries information about the beam This work was supported by the Swedish Research
energy. The importance of this additional information is Council Grant No. 2023-06998. Computational resources
even more pronounced in predicting the x-ray flux (inte- were provided by the Maxwell HPC cluster operated at
grated first order signal), with R2 scores of 0.33 (M1), 0.91 Deutsches Elektronen-Synchrotron DESY, Hamburg,
(M2), and 0.97 (M3). The limited accuracy of M1 indicates Germany.
that the laser parameters alone may not sufficiently capture
the initial beam pointing and resulting trajectory, which in DATA AVAILABILITY
turn provide indirect information on possible losses during The data that support the findings of this article are not
transport. publicly available upon publication because it is not
Despite relying only on noninvasive diagnostics, M2 technically feasible and/or the cost of preparing, depositing,
achieves outstanding accuracy and is highly suitable for and hosting the data would be prohibitive within the terms
applications where destructive electron beam diagnostics of this research project. The data are available from the
are impractical. Still, additionally incorporating the full authors upon reasonable request.
electron spectrum in model M3 delivers the best perfor-
mance, as further demonstrated by the close agreement
between the measured and predicted spectra in Figs. 8(b)
and 8(c). A comprehensive list of model performance with [1] T. Tajima and J. M. Dawson, Laser electron accelerator,
all combinations of input parameters can be found in Phys. Rev. Lett. 43, 267 (1979).
Supplemental Material [36]. [2] E. Esarey, C. B. Schroeder, and W. P. Leemans, Physics of
laser-driven plasma-based electron accelerators, Rev. Mod.
V. CONCLUSION Phys. 81, 1229 (2009).
[3] A. Picksley, J. Stackhouse, C. Benedetti, K. Nakamura,
In summary, we have shown that data-driven models, H. E. Tsai, R. Li, B. Miao, J. E. Shrock, E. Rockafellow,
when carefully designed, can capture the complex interplay H. M. Milchberg, C. B. Schroeder, J. van Tilborg, E.
between laser parameters and electron dynamics in a laser- Esarey, C. G. R. Geddes, and A. J. Gonsalves, Matched
plasma accelerator. Central to our approach is a model guiding and controlled injection in dark-current-free, 10-
architecture that leverages an orthogonal-distance-based GeV-class, channel-guided laser-plasma accelerators,
Phys. Rev. Lett. 133, 255001 (2024).
training objective to effectively minimize the impact of [4] O. Lundh, J. Lim, C. Rechatin, L. Ammoura, A. Ben-
measurement errors in input parameters. The model effec- Ismaïl, X. Davoine, G. Gallot, J.-P. Goddet, E. Lefebvre, V.
tively uncovers the dependencies between laser parameters Malka, and J. Faure, Few femtosecond, few kiloampere
and electrons, eventually enabling tracing variations in the electron bunch produced by a laser–plasma accelerator,
resulting x-ray spectra back to specific changes in the drive Nat. Phys. 7, 219 (2011).
laser. Moreover, our predictive model precisely forecasts [5] Z. Guo, S. Liu, B. Zhou, J. Liu, H. Wang, Y. Pi, X. Wang,
the x-ray spectrum, offering a noninvasive diagnostic that Y. Mo, B. Guo, J. Hua, Y. Wan, and W. Lu, Preclinical
can provide prospective users with shot-to-shot parameters tumor control with a laser-accelerated high-energy electron
of the delivered radiation. The same model framework radiotherapy prototype, Nat. Commun. 16, 1895 (2025).
could also serve as the basis for active feedback loops that [6] J. Wenz, S. Schleede, K. Khrennikov, M. Bech, P. Thibault,
directly stabilize the secondary radiation during operation. M. Heigoldt, F. Pfeiffer, and S. Karsch, Quantitative x-rayphase-contrast microtomography from a compact laser-
Recent advances have highlighted the potential of laser- driven betatron source, Nat. Commun. 6, 7568 (2015).
plasma accelerators as compact drivers of secondary [7] W. Wang, K. Feng, L. Ke, C. Yu, Y. Xu, R. Qi, Y. Chen, Z.
radiation sources. However, their stability, reliability, and Qin, Z. Zhang, M. Fang et al., Free-electron lasing at 27
tunability still fall short of what is routinely achieved with nanometres based on a laser wakefield accelerator, Nature
conventional accelerator technologies. Our approach helps (London) 595, 516 (2021).
101302-12
DATA-DRIVEN MODELING OF A LASER-PLASMA … PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
[8] M. Labat et al., Seeded free-electron laser driven by a [22] A. Sanchez-Gonzalez et al., Accurate prediction of x-ray
compact laser-plasma accelerator, Nat. Photonics 17, 150 pulse properties from a free-electron laser using machine
(2023). learning, Nat. Commun. 8, 15461 (2017).
[9] S. K. Barber, F. Kohrell, C. E. Doss, K. Jensen, C. Berger, [23] C. Emma, A. Edelen, M. Hogan, B. O’Shea, G. White, and
F. Isono, Z. Eisentraut, S. Schröder, A. J. Gonsalves, K. V. Yakimenko, Machine learning-based longitudinal phase
Nakamura, G. R. Plateau, R. A. van Mourik, M. Gracia- space prediction of particle accelerators, Phys. Rev. Accel.
Linares, L. Labun, B. M. Hegelich, S. V. Milton, C. G. R. Beams 21, 112802 (2018).
Geddes, J. Osterhoff, C. B. Schroeder, E. H. Esarey, and J. [24] G. James, D. Witten, T. Hastie, R. Tibshirani, and J. Taylor,
van Tilborg, Greater than 1000-fold gain in a free-electron Statistical learning, inAn Introduction toStatisticalLearning:
laser driven by a laser-plasma accelerator with high With Applications in Python (Springer International Publish-
reliability, Phys. Rev. Lett. 135, 055001 (2025). ing, Cham, 2023), p. 17, 10.1007/978-3-031-38747-0.
[10] H. T. Kim, V. B. Pathak, K. Hong Pae, A. Lifschitz, F. [25] G. Shmueli, To explain or to predict?, Stat. Sci. 25, 289
Sylla, J. H. Shin, C. Hojbota, S. K. Lee, J. H. Sung, H. W. (2010).
Lee, E. Guillaume, C. Thaury, K. Nakajima, J. Vieira, L. O. [26] P. Mehta, M. Bukov, C. Wang, A. G. R. Day, C.
Silva, V. Malka, and C. H. Nam, Stable multi-GeVelectron Richardson, C. K. Fisher, and D. J. Schwab, A high-bias,
accelerator driven by waveform-controlled PW laser low-variance introduction to machine learning for physi-
pulses, Sci. Rep. 7, 10203 (2017). cists, Phys. Rep. 810, 1 (2019).
[11] M. Kirchen, S. Jalas, P. Messner, P. Winkler, T. Eichner, L. [27] R. J. Carroll, D. Ruppert, L. A. Stefanski, and C. M.
Hübner, T. Hülsenbusch, L. Jeppe, T. Parikh, M. Schnepp, Crainiceanu, Measurement Error in Nonlinear Models,
and A. R. Maier, Optimal beam loading in a laser-plasma 2nd ed. (Chapman and Hall/CRC, London, 2006), p. 1,
accelerator, Phys. Rev. Lett. 126, 174801 (2021). 10.1201/9781420010138.
[12] A. R. Maier, N. M. Delbos, T. Eichner, L. Hübner, S. Jalas, [28] J. A. Hutcheon, A. Chiolero, and J. A. Hanley, Random
L. Jeppe, S. W. Jolly, M. Kirchen, V. Leroux, P. Messner, measurement error and regression dilution bias, Br. Med. J.
M. Schnepp, M. Trunk, P. A. Walker, C. Werle, and P. 340, c2289 (2010).
Winkler, Decoding sources of energy variability in a laser- [29] R. J. Adcock, A problem in least squares, Analyst 5, 53
plasma accelerator, Phys. Rev. X 10, 031039 (2020). (1878).
[13] L. Rovige, J. Huijts, I. Andriyash, A. Vernier, V. Tomkus, [30] C. H. Kummell, Reduction of observation equations which
V. Girdauskas, G. Raciukaitis, J. Dudutis, V. Stankevic, P. contain more than one observed quantity, Analyst 6, 97
Gecys, M. Ouille, Z. Cheng, R. Lopez-Martens, and J. (1879).
Faure, Demonstration of stable long-term operation of a [31] L. Leng, T. Zhang, L. Kleinman, and W. Zhu, Ordinary
kilohertz laser-plasma accelerator, Phys. Rev. Accel. least square regression, orthogonal regression, geometric
Beams 23, 093401 (2020). mean regression and their applications in aerosol science, J.
[14] A. Döpp, C. Eberle, S. Howard, F. Irshad, J. Lin, and M. Phys. Conf. Ser. 78, 012084 (2007).
Streeter, Data-driven science and machine learning meth- [32] J. A. Hausman,W. K. Newey, H. Ichimura, and J. L. Powell,
ods in laser–plasma physics, High Power Laser Sci. Eng. Identification and estimation of polynomial errors-in-
11, e55 (2023). variables models, J. Econom. 50, 273 (1991).
[15] M. Streeter et al., Laser wakefield accelerator modelling [33] S. Schennach, Estimation of nonlinear models with meas-
with variational neural networks, High Power Laser Sci. urement error, Econometrica 72, 33 (2004).
Eng. 11, e9 (2023). [34] Z. Hu, Z. T. Ke, and J. S. Liu, Measurement error models:
[16] S. J. D. Dann et al., Laser wakefield acceleration with From nonparametric methods to deep neural networks,
active feedback at 5 Hz, Phys. Rev. Accel. Beams 22, Stat. Sci. 37, 473 (2022).
041303 (2019). [35] K. Hornik, Approximation capabilities of multilayer feed-
[17] A.Hanuka,C. Emma,T.Maxwell, A. S. Fisher, B. Jacobson, forward networks, Neural Netw. 4, 251 (1991).
M. J. Hogan, and Z. Huang, Accurate and confident pre- [36] See Supplemental Material at http://link.aps.org/
diction of electron beam longitudinal properties using supplemental/10.1103/dnnw-lfn5 for details on choice of
spectral virtual diagnostics, Sci. Rep. 11, 2945 (2021). model parameters σx; σy and weight regularization in the
[18] S. Jalas, M. Kirchen, P. Messner, P. Winkler, L. Hübner, J. orthogonal distance model as well as derivation of effects
Dirkwinkel, M. Schnepp, R. Lehe, and A. R. Maier, of chromaticity induced measurement errors and additional
Bayesian optimization of a laser-plasma accelerator, Phys. information on the predictive model.
Rev. Lett. 126, 104801 (2021). [37] K. P.Murphy,Machine learning:Aprobabilistic perspective,
[19] R. J. Shalloo et al., Automation and control of laser in Adaptive Computation and Machine Learning Series
wakefield accelerators using Bayesian optimization, Nat. (MIT Press, Cambridge, MA, 2012), pp. 225–227,
Commun. 11, 6355 (2020). 252.
[20] J. Lin, Q. Qian, J. Murphy, A. Hsu, A. Hero, Y. Ma, [38] S. C. Hora, Aleatory and epistemic uncertainty in proba-
A. G. R. Thomas, and K. Krushelnick, Beyond optimiza- bility elicitation with an example from hazardous waste
tion—supervised learning applications in relativistic laser- management, Reliab. Eng. Syst. Saf. 54, 217 (1996).
plasma experiments, Phys. Plasmas 28, 083102 (2021). [39] M. H. Shaker and E. Hüllermeier, Aleatoric and epistemic
[21] O. Convery, L. Smith, Y. Gal, and A. Hanuka, Uncertainty uncertainty with random forests, in Advances in Intelligent
quantification for virtual diagnostic of particle accelerators, Data Analysis XVIII (Springer International Publishing,
Phys. Rev. Accel. Beams 24, 074602 (2021). 2020), p. 444, 10.1007/978-3-030-44584-3_35.
101302-13
F. BROGREN et al. PHYS. REV. ACCEL. BEAMS 28, 101302 (2025)
[40] B. Efron, Bootstrap methods: Another look at the jackknife, [45] M. Santarsiero, D. Aiello, R. Borghi, and S. Vicalvi,
Ann. Stat. 7, 1 (1979). Focusing of axially symmetric flattened Gaussian beams,
[41] J. Franke and M. H. Neumann, Bootstrapping neural net- J. Mod. Opt. 44, 633 (1997).
works, Neural Comput. 12, 1929 (2000). [46] M. Tzoufras, W. Lu, F. S. Tsung, C. Huang, W. B. Mori, T.
[42] C. Molnar, 8.5 permutation feature importance, in Katsouleas, J. Vieira, R. A. Fonseca, and L. O. Silva, Beam
Interpretable Machine Learning (Leanpub, 2020), loading in the nonlinear regime of plasma-based accel-
https://christophm.github.io/interpretable-ml-book/feature- eration, Phys. Rev. Lett. 101, 145002 (2008).
importance.html. [47] A. Stolow, A. E. Bragg, and D. M. Neumark, Femtosecond
[43] Ensemble-PyTorch, PyTorch (software), version 0.9.1, time-resolved photoelectron spectroscopy, Chem. Rev.
https://github.com/TorchEnsemble-Community/Ensemble- 104, 1719 (2004).
Pytorch. [48] J. Chastain and R. C. King, Jr., Handbook of X-Ray
[44] R. Lehe, M. Kirchen, I. A. Andriyash, B. B. Godfrey, and Photoelectron Spectroscopy, edited by J. Chastain and
J.-L. Vay, A spectral, quasi-cylindrical and dispersion-free R. C. King, Jr(Perkin-Elmer Corporation, Eden Prairie,
particle-in-cell algorithm, Comput. Phys. Commun. 203, MN, 1992), Vol. 40, p. 10.
66 (2016). [49] J. Yano and V. K. Yachandra, X-ray absorption spectros-
copy, Photosynth. Res. 102, 241 (2009).
101302-14
Paper B
F. Brogren, C. Olofsson, J. Magnusson, A. Gonoskov.
"π-PIC: a framework for modular particle-in-cell developments and
simulations". arXiv preprint: 10, 48550.

π-PIC: a framework for modular particle-in-cell developments
and simulations
Frida Brogren, Christoffer Olofsson, Joel Magnusson, and Arkady Gonoskov
Department of Physics, University of Gothenburg, SE-41296 Gothenburg, Sweden
November 13, 2025
Abstract
Recently proposed modifications of the standard particle-in-cell (PIC) method resolve long-standing
limitations such as exact preservation of physically conserved quantities and unbiased ensemble down-
sampling. Such advances pave the way for next-generation PIC codes capable of using lower resolution
and fewer particles per cell, enabling interactive studies on personal computers and facilitating large-
scale parameter scans on supercomputers. Here, we present a Python-controlled framework designed
to stimulate the dissemination and adoption of novel PIC developments by providing unified interfaces
for accommodation, cross-testing, and comparison of PIC solvers and extensions written in Python
or low-level languages like C++ and Fortran. To demonstrate flexibility of proposed interfaces we
present and test implementations of several PIC solvers, as well as of extensions that is capable of
managing QED processes, moving-window, and tight focusing of laser pulses.
1 Introduction
In recent years, impressive progress has been made in developing advanced particle-in-cell (PIC) schemes
that overcome, or mitigate, fundamental limitations of this method. Significant effort has been dedicated
to creating faster, more versatile, exact and user-friendly PIC codes. Some of the most well-used codes
includes (but are not restricted to) PIConGPU [1], Smilei [2], Epoch, Osiris [3], PICADOR [4] and the
codes included in the BLAST project (e.g. WarpX [5], WakeT, HiPACE and FBPIC [6]). Codes, like
Smilei, Osiris, Epoch and WarpX are targeting a wide user base with a large variety of functionality
suitable for different types of interactions, for example ionization, Breit-Wheeler pair creation and ra-
diation reactions, ionization, Schwinger process and collisions. Other codes, including FBPIC, WakeT,
and HiPACE, aim to enhance speed through different strategies, such as employing simplified geometries,
adopting quasi-static methods, and utilizing boosted frame simulations.
Despite these advances, as the PIC community grows and the applicability of PIC codes enlarges the
need for testing and incorporating new physics and algorithms persists. At the same time, as these codes
evolve, several with focus on optimized execution, modifying them becomes more complex and expensive
from a developers perspective. In other areas of computational physics and applied mathematics, such
challenges have been alleviated by standard interfaces that enable testing and comparison of different
methods within a single design pattern (see, e.g. OpenFOAM [7], ASE [8], LAMMPS [9], BLAS [10],
LAPACK [11], FFTW [12], COIN-OR [13] and Scikit-learn [14]).
Certainly, the developers of almost every well-established PIC code strive to cover such functionality,
but for the wider PIC community, such interface remains to be established. Three distinct difficulties
can be identified. The first is related to compatibility of physical units, design patterns and interfaces
that vary across codes and communities. In particular, performance is often prioritized over flexibility
and simplicity, while these aspects appear to be crucial for cross-community efforts. The second difficulty
is centralization of developments: even in the case of open-source development using GitHub or other
platforms, in practice any contribution of a significantly new scheme requires setting up a collaboration
with the lead developers, which can slow down testing unestablished novel ideas. The third difficulty is
the necessity of contributors to have profound knowledge of the programming language chosen for the
project.
In this paper we report on our experience of developing an interface that appears to overcome these
obstacles. Compatibility is addressed by a framework that supports a variety of solvers and extensions,
1
which can be implemented in C++, Fortran or any other language that provide callable functions and,
thereby, does not compromise efficiency. Multithreading can be used directly within this setup, while
GPU and multi-node execution require additional development, which, however does not appear to face
any fundamental obstacles. The approach is decentralized: new solvers and extensions can be developed
completely independently and coupled within Python. The capabilities of the devised interface are
shown by considering a number of solvers, practical extensions and Python scripts that combine them
for performing validation tests.
The paper is arranged as follows: in section 2, we give a brief overview of selected PIC algorithms, nu-
merical inaccuracies connected to the methods and future prospects. Section 3 considers architecture and
utilization of the presented framework called π-PIC. Section 4 presents extensions for beam-plasma inter-
action and in section 5.1 the development and testing of an electrostatic solver. To further demonstrate
opportunities for incorporating non-conventional solvers in section section 5.2 we describe an explicit,
relativistic PIC solver, called emc2, that preserves energy and has improved properties with respect to
momentum conservation. In section 5.3, we show benchmark of a energy-conserving solver [15] imple-
mented in π-PIC against Smilei [2] in the context of laser-wakefield acceleration.
2 Limitations and prospects for the PIC method
The challenges of simulating plasma typically revolve around reducing numerical artifacts or increasing
computational efficiency. This section begins with a brief overview of the fundamentals of the PIC
algorithm, then turns to a discussion of various solvers, their role in mitigating or damping numerical
artifacts and instabilities, and the implications these choices have for computational performance.
The PIC algorithm is normally divided into four steps: (1) particle push, where the particles are
advanced using the Lorentz equation, (2) charge deposition, where the particle position and momentum
are used to calculate current and densities, (3) field advancement, where the electromagnetic field is
advanced and (4) field gathering, where the field is interpolated back to the particle position. Important
variations of this scheme includes the type of charge deposition, choice of field solver, macro-particle
shape function and particle push.
Methods are typically grouped into relativistic or non-relativistic and electrostatic or electromagnetic.
Deciding whether to use relativistic or non-relativistic and electrostatic or electromagnetic codes largely
depends on the plasma state being simulated and the governing physics. Additionally, methods can be
classified based on their implementation as implicit, explicit, spectral, global and local. Implicit methods
evolve some physical state by solving a system of equations, usually in the form of a matrix equation,
an operation which often is done iteratively. In contrast, explicit methods advance the state using only
the current state and past states of the system. On the other hand, spectral methods are based on an
explicit field advancement in Fourier space representation. Global methods employs a global calculation
to advance the state while local methods only uses the state of the current and neighboring nodes.
Spectral methods are inherently global due to the Fourier transform, while explicit and implicit methods
can be either local or global. More generally, PIC algorithms can use global or local and implicit or
explicit routines, using coordinate or spectral representation in different parts of the algorithm.
Due to their parallelization capability and fast execution, explicit local methods typically achieve
higher computational speeds than other approaches. Moreover, the rapid growth of parallel computing
resources has driven a broad preference for explicit methods. However, generally, these methods are less
stable compared to implicit (or spectral) methods which can allow larger time step and better conservation
of important quantities like energy and charge. Typically, numerical artifacts stem from inaccuracies in
numerical schemes which changes the dynamics of the system and leads to an accumulation of errors
manifesting as a violation in the laws of physics. In the following we examine three examples of this,
highlighting the trade-off between accuracy and computational speed: (1) numerical dispersion, (2) charge
accumulation, and (3) momentum and energy non-conservation.
Numerical dispersion emerges as a result of the discretization of electromagnetic waves, causing the
discretized waves to travel at velocities that differ from the speed of light. A popular choice of Maxwell
solver is Finite-Difference Time-Domain (FDTD) [16] which shows direction-dependent numerical disper-
sion. Suggestion for a modified FDTD schemes which mitigates numerical dispersion completely in some
directions and partly in others came in 1999 [17] and since then more have followed [18–21]. Setting aside
FDTD maxwell-solver, numerical dispersion can be entirely avoided by using a spectral solver [22] which
solves the spatial part of maxwells equations in Fourier space and advances the field in the time-domain.
However, this is done at the expense of doing a global Fourier transform operation. For this reason,
finite-difference methods was long preferred over spectral solvers. Yet, in 2013 a parallelized implementa-
2
tion using domain decomposition was presented in [23] and has since then been popularized and further
developed in codes such as WarpX [5] and FBPIC [6].
Apart from low numerical dispersion other desirable properties of PIC codes are charge conservation.
Charge conservation is equivalent to fulfilling the continuity equation. Ampères-Maxwells law is used to
advance the field in time, however, Gauss’s law is not included, as it does not depict the evolution of
the state. Consequently, one tactic to ensure charge conservation is by correcting the fields retroactively
using Gauss law [22,24,25]. Another option is to adopt a so called current deposition scheme that ensures
current conservation locally [26–30]. Again, the second options has some benefits as it is a local operation
while the correction of Gauss law must be carried out globally.
Furthermore, the PIC scheme in its standard implementation does not conserve energy nor momentum.
This may result in numerical heating, in particular if the Debye length is not resolved aliased modes
will appear as thermal fluctuations [31–33]. One method to address this is employing smoother shape
functions, acting as a low-pass filter for small wavelength oscillations [22], essentially reducing spatial
invariance thereby achieving better momentum conservation. Another way to limit numerical heating
is to employ energy conserving schemes. In general these are constructed using an implicit solver [34].
π-PIC however comes equipped with an explicit energy conserving (ec and ec2) solver [15]. Another
promising approach for achieving energy conservation has been recently proposed by Rickertson and Hu
in [35] and demonstrated for the case of non-relativistic PIC simulations.
The preceding discussion indicates a connection between the accumulation of numerical artifacts and
the failure of numerical methods to comply with conservation laws. Unsurprisingly, a discretization of
a continuous theory with a set of symmetries and associated conservation laws, does not automatically
inherent these symmetries. However, in the past decade a new research area for the design of PIC codes
which inherently respect conservation laws has emerged. So called structure-preserving codes [36] uses
geometric integrators to achieved exact charge conservation by building algorithms based on principles
of gauge invariance [37, 38]. Similar techniques can also provide energy conservation [39]. We consider
this line of research a good prospect for further development of the PIC algorithm.
To summarize, a diverse range of PIC implementations is currently available. To further advance the
field, it is necessary to develop a unified platform where different code variations can be tested, compared
and studied to reveal their preferable domains of use. The objective of π-PIC is to support a wide
range of diverse solvers, while ensuring they remain compatible with extensions that enable additional
functionalities, such as ionization injection, radiation reaction and strong-field QED [40, 41], as well as
ensemble downsampling [42,43].
3 The π-PIC framework
To accommodate maximum freedom for the user, the π-PIC framework offers three interaction levels that
can be used to modify the code: (1) the Python interface, (2) extensions and (3) solvers.
The first, entry level, offers functionality for initializing the simulation with arbitrary fields and
particle densities as well as advancing the system state, see fig. 1. Additionally, it provides runtime read
capabilities for particle phase-space, and both runtime read and modification access of fields. The PIC
algorithm itself is written in C++ and the linking between Python and C++ is realized using pybind11
[44]. Custom functions, defined in the Python interface used to preform actions on the simulation state,
are compiled from Python into C++ code using the Numba decorator @cfunc(). The decorator indicates
parts of the code that are precompiled to avoid a potential loss of efficiency in performance-critical parts.
There are three types of predefined decorator functions available in π-PIC: add particle callback,
field loop callback and particle loop callback. Each of these callbacks accepts a predefined
set of arguments (such as particle position, weight, and electric field), which are described in the example
script (section 4).
The final two arguments of every π-PIC callback, data double and data int, correspond to C++
pointers referencing a float64 array and an int32, respectively. These can be used to pass memory
addresses for storing or retrieving data (see fig. 1). For instance, one might supply the address of
a Numpy array to which the electric field should be written. This mechanism enables dynamic data
exchange between Python and the C++ backend across PIC update steps.
The second interaction level is the extensions, which are added in the Python interface and then
attached to each execution in connection to each PIC update. Extensions can be written in Python
or in C++ and provides a handle for further modification of particle state as well as field state. This
is particularly useful for developing, for example, new ionization routines or incorporating collisional
3
and strong-field effects. The development and application of extensions are further discussed and ex-
emplified in section 4. The extensions are added to the advance routine in the Python interface by the
add handler() call, see fig. 1.(L:8). This places the defined handlers in a list in the handlerManager
which is then called during each advance call of the simulation. A step-by-step tutorial on how to develop
and employ an extension can be found at [45].
The final interaction level is dedicated to the employment of new solvers. In π-PIC solvers are
composed of two defining classes: A pic solver and a field solver that specifies the advancement
of particle states and fields respectively. The field solver is a mandatory member of the pic solver
class, reflecting the fact that field and particle updates are generally coupled and must be advanced
consistently. This composition enforces a clear ownership hierarchy: the pic solver manages the
field solver. Developers may still implement or replace individual components (e.g. particle pusher)
while reusing an existing field solver or pic solver. Furthermore, this class structure allows for a
large degree of freedom in terms of the implementation of the PIC update, for example, part of the field
update can be done concurrently with the particle update, as is the case for the energy conserving (EC)
solver [15] implemented in π-PIC.
To maintain compatibility with the rest of the framework, there is a number of mandatory and optional
class methods that every field solver and pic solver implements. These are defined in the header
file named interfaces.h. Provided that the new solver follows the same structure as the template, it
will be compatible with all previously developed extensions.
With the flexibility granted by this framework, the only fundamental limitation imposed on the
developer is that particles are stored in the CellInterface on a cell-by-cell basis. This ensures that
extensions can access particles efficiently. Fields, in contrast, are accessed by extension through the
field solver, allowing developers to define and update them according to their needs. This enables
support for staggered grids and non-Cartesian geometries. In section 5 we will further describe the
implementation of solvers by the example of an electrostatic solver, as well as testing the accuracy of
solvers in simulating Laser Wakefield acceleration.
4 Extensions
In this section, we describe the role of extensions in π-PIC, along with some of the extensions, useful
for laser-plasma interaction, that have already been implemented. Extensions provide functionality for
modifications of particles and fields which can be applied on top of the chosen PIC solver.
To access particles and fields during the advance call extensions has two defining methods: a
handler, for modifying, adding and removing particles, and a fieldHandler for modifying fields.
The input arguments to the handler handler is repacked into a cellInterface for easy access to cell
data, such as fields, cell dimension, cell node position, number of particles and particle buffer size. At
the advance call, the ensemble loops over the existing cells in a multi-threaded process, applying the
action of handler. In handler particles can be removed by setting their weight to zero and added by
writing to the particle buffer. Additionally, the particle momentum can be arbitrarily modified, however,
the position should be changed with care, since particle migration to other cells can cause the particle
buffer to overflow. The fieldHandler has input arguments, such as field, position and index. As for
the handler the fieldHandler method is applied at the advance call to all points in the field grid. It
is, however, the fieldSolver which defines how the loop is executed, enabling different types of grid
geometries and staggering.
Here we present two extensions which utilizes both the handler and fieldHandler: an imple-
mentation of absorbing boundaries with replacement of particles (section 4.1) and a moving window
(section 4.2). Additionally, we present an extension for mapping a spatially extended pulse to its focus
(section 4.3).
4.1 Absorbing boundaries
Periodic and absorbing (or open) boundary conditions are arguably the most common choice in plasma
simulations. Generally, implementing periodic boundary conditions is straightforward, as it mainly in-
volves linking the ends of the simulation box. Free boundary conditions are more challenging, as it
theoretically implies placing the boundary at infinity or restricting the type of electro-magnetic waves
allowed at the boundary - namely no incoming waves. For this reason, open boundary conditions are most
commonly emulated by introducing an absorbing medium at the boundary. A popular implementation is
Perfectly Matched Layers (PML) [46], where, so called, complex coordinate stretching is used to emulate
4
Figure 1: To the left, a simplified main script for simulation with π-PIC is shown, a full version can be
seen in [45]. Line 1 initiates the simulation instance. Line 2-5 adds particles and fields to the simulation.
Line 6-7 defines functions which loops over fields and particles for saving or modification. All predefined
decorator functions in the Python interface has a data double and data int as arguments. This
allows the user to pass an address to a Numpy array which can be used to record the particle or field
state, see line 10-11. Line 8 adds the extension to the advance call. During the advance call (line 9)
the extensions and solver are applied to each cell and associated particles in a multi-threaded process.
Line 9 advances the simulation state, executing the PIC solver and field solver and added extensions.
The use omp key word is used to activate multi-threading. The existence of multi-threaded processes is
denoted as double arrow in the right scheme. Line 9 and 10 loops through the particles and field, passing
an address to the double arrays particle dd and field dd, respectively, for saving the simulation
state. The different instances of the code in the scheme to the right are color coded to match the Python
commands that activates them.
an absorbing space within a few grid points. Here, we take on a minimalistic approach known as masking,
where incoming electromagnetic fields are gradually depleted by iterative attenuation [47].
For our implementation of free boundary condition the field is attenuated at every iteration by
Ẽ(x) = R(x)E(x), B̃(x) = R(x)B(x) (1)
where R(x) is the masking function and Ẽ, B̃ are the updated fields. To avoid absorption changing with
the length of timestep, we require that the cumulative attenuation after n time steps ∆t, to be equal
to the attenuation obtained using a single time step of size n∆t, that is R(x,∆t)n = R(x, n∆t). This
condition is full-filled by the choice
R(x,∆t) = e∆tr(x). (2)
Moreover, to minimize reflection we demand that r(x ) = 0 and dr(xb)b dx = 0, where xb is the start of the
boundary. Additionally we require that limx→x R(x) → 0 which gives the condition limmax x→x r(x) →max
−∞. A function that satisfies these conditions is
[ ]
r(xr) = α cos(πxr/2)−
1
, (3)
cos(πxr/2)
5
where x = x−xbr x −x is the relative, normalized position at the boundary such that xr ∈ [0, 1] and α is amax b
shape parameter controlling the steepness of the absorption, see fig. 2.(b). We assess the absorption by
letting a Gaussian beam with wavelength λ = 1 cm pulse duration ts = 5 fs and equal transverse extent,
propagate at an angle of 45◦ towards the boundary as shown in fig. 2.
Figure 2: Panel (a) shows the initial field of the simulation: a Gaussian pulse is propagating at 45◦
angle towards the boundary where the masking function R(y) is shown in orange color scale. Panel (b)
shows the difference in absorption function with different shape factors α. The red curves corresponds to
boundaries with maximal absorption as shown in panel (c). For these simulations the boundary size was
set to four wavelengths. Panel (d) shows the dependence of the absorption on boundary size. The shape
factor corresponding to ∆tα = 1 was used for these simulations. The absorption was calculated as the
decrease in total electromagnetic energy before and after interaction with the boundary.
The disadvantage of this approach, compared to other implementations, such as PMLs, is that it
requires a large damping region. Here, we find that a damping region > 4 wavelengths achieves sufficient
absorption.
To fully emulate open boundaries in plasma simulations particles needs to be added and removed to
maintain a constant density and velocity distribution at the boundary. Similar to fields, this must be
done gradually and smoothly to avoid introducing artifacts and high-frequency noise at the boundary.
We propose removing particles at the same rate as fields, while simultaneously adding particles with
velocity distributions determined by temperature.
Since the primary reason for this work was to showcase the use of extensions to modify fields and
particles, no comparative study with other methods has been made. These comparisons, as well as the
possible implementation of more advanced methods, are left for future work. The code implementation
can be found at [45].
4.2 Moving window
For simulation involving, for example, laser-plasma interaction the computational demands can typically
be greatly reduced by only simulating a region of space that moves with the process of interest. For
6
example, this is the case when simulating a laser pulse or a relativistic bunch of electrons. To achieve
this a moving window is used.
The available solvers in π-PIC all have spectral field solvers, consequently the default boundary
condition is periodic and a propagating field will reappear at the opposite side of the simulation box if
it passes through the boundary. This makes the implementation of the moving window rather simple.
For a moving window with speed v we need to remove the fields and all plasma perturbation in a region
which moves around the simulation box with the same speed. The full implementation and use of the
moving window is documented in [45].
4.3 Mapping of tightly focused laser pulses
Tightly focused laser pulses are central to studies of high-intensity laser–matter interactions [48,49], fun-
damental quantum systems [50] and various areas of attosecond physics [51]. While analytical methods
exist [52–56] to describe the propagation of their associated electromagnetic fields, numerical approaches
provide a more direct and flexible alternative, free from the inaccuracies introduced by simplifying assump-
tions. Such computations are particularly relevant for vacuum electron acceleration [57], and radiation
generation driven by laser–electron dynamics, where interactions are highly sensitive to the structure
of the focal electromagnetic field [58–60]. However, numerical techniques based on direct integration of
Green’s functions or FDTD methods can be computationally intensive and prone to numerical dispersion.
Spectral field solvers on the other hand, as implemented in π-PIC, leverage fast Fourier transforms to
evolve the electromagnetic field. This enables dispersion-free computation of the focal field of tightly
focused laser pulses, given a known far-field profile, within a single time step. Despite this advantage, the
simulation box ΩB must be sufficiently large to encompass the initial laser pulse profile. Consequently,
excessive computational effort is spent evolving the simulation in regions where no significant dynamics
occur, even though the physical interaction of interest is many times confined to a much smaller box
Ωb ⊆ ΩB (see fig. 3). To avoid this computational overhead, the periodic boundary conditions imposed
by spectral solvers can be leveraged. In what follows, we further develop the approach proposed in [61]:
since the divergence of previously focused radiation defined in a space with cyclic geometry results in
overlaying and summing of electromagnetic radiation, one can time reverse this process and obtain a
focused radiation by time advancement of summed, overlaid radiation before focusing.
Indeed, any point r = (x, y, z) inside a given simulation box is equivalent to a point outside the
computational domain according to the identification
r ∼ r+ L, L = (Lx, Ly, Lz), (4)
where Lx, Ly and Lz are the simulation box lengths along each coordinate axis, respectively. In other
words, any two coordinates that differ by a multiple of the box lengths, represent the same physical
point. Moreover, eq. (4) implies that electromagnetic fields initialized in ΩB admit an equivalent repre-
sentation within Ωb. This suggests that the initial pulse can be mapped directly onto Ωb, after which
the electromagnetic field can be evolved over a single time step, thereby avoiding the need to include the
full extent of the far-field. To perform such a mapping, a set of coordinates enclosing the initial pulse
profile, relative to the simulation box Ωb must be defined. One approach is to assume that the laser pulse
initially occupies a spherical segment with thickness ℓ, aperture angle θmax, and radius R0. The latter is
interpreted as the distance from the pulse centroid to the location of the focus. Then, letting i, j, k ∈ Z
denote integer steps that constitutes the box Ωb, the electric (and magnetic) field at any location in Ωb
can be written as
h∑=h+ h∑=h
+
h∑=h+x x y y z z
E (i∆x, j∆y, k∆z) = E ([hxnx + i] ∆x, [hyny + j] ∆y, [hznz + k]∆z) , (5)
h − − −x=hx hy=hy hz=hz
where ∆x,∆y,∆z and nx, ny, nz are the spatial resolution and number of cells along each coordinate
axes of Ωb, respectively. Here, the integers hm with m ∈ {x, y, z} are referred to as hypercells, which are
periodic replicas of Ωb arranged to fully enclose the initial pulse profile. The hypercell bound for each
coordinate axis is given by ⌊ ⌋
± Lm/2± (R0 + ℓ/2)hm = . (6)Lm
An example of Ωb is depicted in fig. 3(a) as the shaded quadratic region, whereas the three bottom rows
of unshaded squares demarcate the required hypercells to accurately capture the pulse profile. To avoid
7
Figure 3: Electric field strength of the initial laser pulse after initialization within (a) ΩB and (c) Ωb,
along with the corresponding field states after a time step R0/c shown in (b) and (d), respectively. In
(a), the shaded gray square indicates the extent of Ωb, while the surrounding unshaded squares represent
equivalent coordinates within ΩB . Panel (e) shows the average runtime across 15 runs for the initialization
and propagation of the electromagnetic fields as a function of R0.
superimposing negligible electromagnetic fields, arising from e.g. empty cells, only the coordinates that
satisfy
| ℓR0 − r| ≤ , θ = arccos(r̂ · n̂) ≤ θmax, (7)
2
are included. Here, r = |r|, r̂ = r/r, and n̂ is the unit vector along the direction of R0. The specific
implementation of the periodic mapping, including an optimized procedure, is provided in [45], and is
demonstrated next.
To demonstrate the periodic mapping, we perform 2D simulations for a linearly polarized focused
pulse with central wavelength λ. We take the axis of propagation to be along y, using a dimensionless
far-field pulse profile u analogous to that given in [61]:
u∥(r −R0)
u(r) = u⊥(θ), (8)
r
where ( ) ( ) ( )
2πξ πξ ℓ ℓ
u∥(ξ) = sin cos
2 Π ξ,− , , (9)
λ ℓ 2 2
is the longitudinal profile and Π (ξ, ξmin, ξmax) is the rectangular function that equals unity whenever
ξ ∈ [ξmin, ξmax] and zero otherwise. The transverse shape is given by
( ) ( ( ))
− ε 2 π θ − θ
ε ( )
max + 2 − ε εu⊥(θ) = Π θ, 0, θmax + cos Π θ, θmax , θmax + , (10)
2 2ε 2 2
where ε = 0.1 is an angle that allows for a smooth transition of eq. (10) around θmax. As an example, we
initialize the pulse using θ ◦max = 59 , ℓ = 2λ, R0 = 16λ, and advance the electromagnetic field described by
eqs. (8) and (10) within a quadratic box ΩB with side lengths 2(R0+ℓ) using the field loop callback
in π-PIC. Then, a separate pulse initialization is done using the periodic mapping implementation (re-
ferred to as periodic in fig. 3) within a smaller domain Ωb with side lengths 4λ. Both configurations have
8
a spatial resolution of 16 cells per wavelength and are evolved over a single time step R0/c, yielding the
initial and final electromagnetic field states as shown in fig. 3(a)–(d). The relative root-mean-square error
for the electromagnetic energy density computed within the focal region, was found to be 1.43 × 10−4,
demonstrating the accuracy of the periodic mapping implementation. To assess the efficiency of periodic
mapping, we benchmark the runtime using 8 Intel(R) Core(TM) i7-10700K CPUs for initializing the field
and performing a single time step R0/c as a function of R0, while maintaining the same spatial resolution.
This time, the box ΩB is set to have dimensions Lx = R0+ ℓ/2+2λ and Ly = 2R0 sin (ε+ θmax) so as to
precisely enclose Ωb and the pulse profile u. The resulting benchmark, presented in fig. 3(e), demonstrates
a significant speed-up for the periodic mapping algorithm across all values of R0.
5 Solvers
In π-PIC, the solver class is an object containing all code relevant to the PIC update. This functionality is,
as mentioned in section 3, realized by two defining classes pic solver and field solver. Both classes
has a number of mandatory and optional function with which different functionality can be implemented.
The field solver has two mandatory methods: fieldLoop and customFieldLoop. The first
implement a loop over the whole grid and the second a loop over a subset of coordinates. These methods
are used by the Python interface and extensions to access data, it is therefore important that they are
implemented by the developer since different solvers can feature different field grids.
The pic solver class defines a single mandatory method, advance. Although its implementa-
tion is left to the developer, the recommended approach is to base it on the advance singleLoop
method in ensemble.h, since this routine handles additionally calls the extensions in a thread-safe
and efficient manner. Due to the template-based structure of the code, and in order to avoid run-
time compilation overhead, the developer is required to implement this function explicitly. The method
advance singleLoop will call the optional methods of pic solver: preStep(), preLoop(),
startSubLoop(), processParticle(), endSubLoop(), postLoop() and postStep() and,
thereby, defines the execution order as below. Calls related to extensions, defined by the handlerManager,
are highlighted in blue.
1: preStep()
2: for n in numberOfIterations:
3: apply fieldHandlers()
4: preLoop()
5: for cell in CellInterface:
6: apply actOnCellHandlers()
7: for particle type in cell:
8: startSubLoop()
9: apply particleHandlers()
10: for particle in particles of this type in cell:
11: processParticle()
12: endSubLoop()
13: postLoop()
14: postStep()
By organizing the PIC update into the optional methods of pic solver, a wide range of PIC
algorithms can be implemented while maintaining compatibility with extensions. As an example, we
next present an electrostatic solver (section 5.1) that demonstrates how a solver can be constructed using
this framework. The energy-conserving solver in π-PIC [15] can be complemented with an approximate
momentum conservation scheme. This addition is described in section 5.2. We proceed by simulating
laser wakefield acceleration (section 5.3) using two additional solvers of π-PIC: the Fourier-Boris (FB)
solver – a spectral electromagnetic solver with Boris pusher – and the energy conserving (EC) solver [15].
As a benchmark, we compare the results with the same simulation performed by the PIC code Smilei [2]
using the M4 solver [62].
5.1 Example: A electrostatic solver
To demonstrate the implementation of solvers a one-dimensional electrostatic solver was implemented
in π-PIC. The solver advances the plasma state according to the electrostatic Vlasov equation and the
9
Ampere-Maxwell equation assuming zero magnetic field
∂f ∂f qEx ∂f
+ vx + = 0, (11)
∂t ∂x m ∂vx
∂Ex
+ 4πjx = 0. (12)
∂t
The particle push corresponding to this system is ∂ v = qExt m , which can be discretized using a forward
finite difference scheme with time staggering
∆tqEj(xj )
vj+1/2
x p
p = v
j−1/2
p + , (13)m
xj+1 = xj +∆tvj+1/2p p p (14)
where upper index denotes time step indices and lower particle index. The time staggering is realized
− ∆tqE01/2 (x0)
through initializing v half a timestep backwards in time as v = v0 − x pp p 2m .
The electric field is advanced as
Ej+1x,i = E
j
x,i −∆t4πjj+1x,i (15)
where upper indices denotes time indices and lower space indices. The current is interpolated to the grid
from the particle position and momentum as
∑N
jj+1x,i = wqS(x
j+1 − x )vj+1/2p i p (16)
p
where N is the number of particles, S is the shape function (Cloud-In-Cell), w is the macro-particle
weight and the lower index denotes particle index. The electric field is interpolated back to the particle
position similarly
∑Ng
Ejx(x) = S(xi − x)Ejx,i (17)
i
where Ng are the number of adjacent grid nodes to x.
The relevant code for pic solver and field solver can be found at [45], however, in short this
is how the methods in pic solver was used to execute the algorithm: Apart from evolving the system
1: preStep() - Pull back particle velocities ∆t/2 to implement staggered grid.
2: preLoop() - Advance Ex and set current to 0 everywhere.
3: processParticle() - Interpolate field to particle position, advance particle and deposit current.
4: postStep() - Advance particle velocities ∆t/2 to remove staggering.
according to eqs. (11) and (12), to accurately model the electrostatic system, the discretization must also
satisfy }
∇ ·E = 4πρ 21D ∂ ϕ
∇×E = 0 ↔ E = −∇ −−→ 4πρ+ = 0.ϕ ∂x2
This requires the system to be initialize in a state consistent with these equations. As the simulation
progresses, this relation will be approximately maintained since the continuity equation (derivable from
eq. (11) by integrating with respect to velocity) together with eq. (12) imposes
( ) ( )
∂ 1 ∂ ∂ 1 ∂2ϕ
ρ− Ex = ρ+ = 0. (18)
∂t 4π ∂x ∂t 4π ∂x2
However, this restricts the plasma systems which can be realistically simulated using periodic boundary
conditions to systems with zero spatial average charge and current. Nevertheless, with this simple set
of equations is possible to simulate quasi-electrostatic phenomenons such as plasma oscillations and the
early stages of charge-current neutral two-stream instability as demonstrated in fig. 4.
10
Figure 4: A two-stream instability by initializing opposite and equal currents flowing in a constant
background field. The instability was simulated using two solvers: (left) electrostatic solver and (right)
a electromagnetic solver (spectral field solver with Boris particle pusher). The density perturbation in
phase space, n, is normalized to the initial particle density n0. The electrostatic solver resolves the
excitation of the two stream instability, but as the instability grows - the electrostatic simulation starts
to gain total energy and eventually becomes invalid.
5.2 Energy-conservative and momentum semi-conservative solver
In this section, we describe an explicit, relativistic, second-order PIC solver that exactly preserves energy
and is improved in relation to momentum conservation. The idea is based on splitting the deviation
from exact preservation of these physically conserved quantities into three fundamental contributions
and eliminating two of them. In doing this analysis, we indicate the path towards eliminating the last
contribution as well. That is why this solver represents an essential step towards a possible way for
simultaneous conservation of both energy and momentum.
Our approach is based on a modification of the PIC method proposed and used to achieve energy
conservation in [15]. This modification exploits a particle-wise splitting of the Maxwell-Vlasov set of
equations such that the current in Maxwell’s equations and the electric component of the Lorentz force
are accounted for within a single subsystem for each particle. This means that the state of each particle
is coupled with the state of electromagnetic field in the nearby nodes, resulting in the evolution equation
that can be solved with enforced energy conservation. Here, we extend this approach to also enforce the
exact conservation of momentum related to the magnetic component of the Lorentz force.
We start by considering the energy and momentum exchange between the electromagnetic field and
a single charge. For this auxiliary analysis, we assume that the charge is distributed within a rigid body,
the center of which is at rb. By considering the effect of all the terms in the Maxwell-Vlasov set of
equations, we separate the part that is responsible for the energy and momentum exchange (CGS units
11
are used):
∂
E (r) = −4πJ (r) , (19)
∂t
c
J (r) = ρ (r− rb) pb, (20)
∫ γ ( )
∂ 1 c
pb = ρ(r− rb) E (r) + pb ×B (r) dr, (21)
∂t mbc V γb
where B, E are magn√etic and electric field vectors; J is the charge current; c is the speed of light, mb,
ρ(r−rb), pb and γ = 1 + p2b are the mass, charge density, momentum (given in units ofmbc) and gamma
factor of the rigid body in consideration, while Vb denotes the part of the space occupied by the body.
The first equation describes the account of current in Maxwell’s equations, the second is the definition of
current and the third describes cumulative effect of the Lorentz force on the change of momentum. To
see the origin of energy conservation in this system, we consider time derivative of the total energy W
being a sum of electromagnetic energy and the kinetic energy of the body in question:
(∫ )
∂ ∂ E2(r) +B2(r) !
W = dr +m 2bc γ = 0. (22)
∂t ∂t 8π
Since the magnetic field is not altered by any of eqs. (19) to (21) and electric field is modified only within
Vb, we can write ∫
! ∂ 1 · ∂E m
2
bc · ∂pb0 = W = E dr + pb . (23)
∂t V 4π ∂t γ ∂tb
Substituting the electric component of the Lorentz force to the second term (the magnetic component
does no work) and using eqs. (19) and (20) to express ∂E/∂t in the first term, we obtain:
∫ (∫ )
∂ − c cW = E · pbρ(r− rb)dr + pb · Eρ (r− rb) dr = 0. (24)
∂t V γ γb Vb
This indicates that the system given by eqs. (19) to (21) describes the changed of field and particle states
in accordance with energy conservation, while the electric component of Lorentz force is the one that
is responsible for this. To see the role of magnetic component of the Lorenz force, we can consider it
separately: ∫
∂ (B) 1 − cpb = ρ(r rb) pb ×Bdr. (25)∂t mbc V γb
Using eqs. (19) and (20), we obtain (note that B is not changed by eqs. (19) to (21)):
∫ ∫
∂ (B) × ∂ 1pb = J Bdr = − E×Bdr. (26)∂t V ∂t 4πb
On the right hand side, to within a factor of c, we have time derivative of the Poynting vector (c/4π)E×B,
which describes energy flux of electromagnetic field. Since a photon carrying energy ℏω carries momentum
ℏk = (ℏω/c)k/k, we can interpret this as an equation that describes the transfer of momentum carried
by electromagnetic waves to pb and vice versa. The remaining part of momentum conservation concerns
the exchange of momentum between the body and electrostatic field.
In such a way we can observe three essential channels for the deviation from the exact conservation
of energy and momentum. The first one concerns energy conservation and can be eliminated by coupling
the update of particle momentum due to the electric component of the Lorentz force with the update
of electric field due to current induced by this particle (this has been done in [15]). The second channel
concerns momentum exchange between photons and the particle in question (eq. (26)). Analogously, this
channel can be eliminated by coupling the update of electric field with the update of particle momentum,
but now due to both electric and magnetic components of the Lorentz force. Finally, the last channel
concerns momentum exchange between the particle and electrostatic part of the field. Eliminating this
last channel requires adding to the subsystem the account for the charge deposition that must be coupled
with current in a way that exactly complies with the charge continuity equation. Unfortunately, exact
local agreement with the charge continuity equation is inconsistent with the Fourier solver (approximate
agreement can be enhanced by the use of higher-order particle shapes). Nevertheless, it is interesting
to consider the use of Esirkepov charge conservation scheme [63] in combination with the consecrative
12
update of electromagnetic field within the framework of finite element exterior calculus [64,65]. Leaving
this consideration for future work, we here describe the elimination of the second channel that can be
seen as an essential step for this approach.
To eliminate the first and the second channels of deviations from the energy-momentum conservations
we need to find exact solution of the following system:
∂p q q
= E+ p×B, (27)
∂t mc mcγ
∂E
= −4πJ, (28)
∂t
where q and m are the charge and mass of the particle, whereas J, E and B are the weighted current,
electric and magnetic field, respectively. The weighting is defined by weight coefficients ci (given the use
of Fourier method we here assume collocated grids for all quantities and use subscript to indicate grid
node):
∑
E = ciEi, (29)
∑
B = ciBi, (30)
∑ qc p
J = ci , (31)
Vg γ
where Vg is the volume of a grid cell. Taking time derivative of eq. (27) and substituting E/∂t from
eq. (28), we obtain:
( )
∂2 −η q ∂p = p+ p ×B, (32)
∂t2 ∑γ mcγ ∂t
4πq2 c2
η = i . (33)
mVg
We use the approach of effective mass, considering γ fixed under this update and using correction later
(see [15]). Introducing rescaled time, electric and magnetic field:
√
η
τ = t, (34)
γ
q
b = √ B, (35)
mc ηγ
√
γ q
e = E (36)
η mc
we obtain
p̈ = −p+ ṗ× b, (37)
ṗ = e+ p× b, (38)
where we use dot to denote time derivative with respect to rescaled time τ . We now need to solve eq. (37)
with initial conditions given by particle’s initial p and ṗ defined via eq. (38) for the initial electric field.
The value of ṗ after one time step defines the updated value of the electric field.
To solve eq. (37), we introduce an orthogonal basis with one axis being parallel to b. Let us denote
the component of b along this axis by p3 and the two other components by p1 and p2. Then for p3 we
have harmonic oscillator equation:
p̈3 + p3 = 0, (39)
the exact solution of which can be expressed using sine and cosine function computation for the given
time step. For the other two components we obtain:
p̈1 = −p1 + ṗ2b, (40)
p̈2 = −p2 − ṗ1b. (41)
Introducing four-component vector
T
s = (ṗ1, ṗ2, p1, p2) , (42)
13
we can transform these set of equation into matrix form
 
0 b −1 0
−b 0 0 −1ṡ = As, A =   . (43)1 0 0 0 
0 1 0 0
Following standard procedure we transform this matrix equation into the one with diagonal matrix by
transforming the basis. Under such transformation the diagonal elements have the form:
i ( ( ) )1/2
w = ±√ b2 2 1/21,2,3,4 + 2± b b + 4 . (44)
2
The exact solution can then be expressed using two sine and cosine function computations. This proce-
dure is implemented in π-PIC for the second-order solver that is available under the name emc2. The
implementation follows the outline patterns [45].
5.3 Benchmark: Laser wake-field acceleration
PIC codes are widely used for the simulation of laser wakefield acceleration (LWFA). Since its proposal
[66], the acceleration technique has attracted much attention because it can sustain accelerating gradients
that far exceed those produced by conventional accelerators [67], thereby enabling compact sources of
relativistic particles.
In this section, we will benchmark simulations of LWFA using the energy-conserving solver [15] and
the Fourier-Boris (FB) solver of π-PIC, comparing it against the well-established PIC-code Smilei [2].
In particular, we will focus on performance at low-resolution. The FB solver is a standard spectral
solver with a Boris particle-pusher. For simulations with Smilei we use the M4 maxwell solver [62] since
it exhibits reduced numerical dispersion, which slows down the propagation of electro-magnetic waves,
compared to the standard FDTD solver. In contrast the energy-conserving and FB solver in π-PIC uses
as spectral solver, which inherently avoids this problem.
We begin by simulating the one-dimensional wakefield excitation that occurs when a laser interacts
with a plasma density upramp.. We employ a laser with wavelength λ = 1 µm and a Gaussian time-
envelope spanning 30 laser cycles at full width half maximum (FWHM). The laser envelope has a peak
amplitude a0 = e|E|/(mecω) = 4, where e denotes the electron charge, E the electric field amplitude, me
the electron mass, c the speed of light and ω the laser frequency. The plasma density ramps from 0 to
4 · 1018 cm−3 in 210λ = 1024 µm. The simulations have a space step ∆x = λ/16, time step ∆t = λ/(64c)
and 8 particles per cell.
The 1D simulations, shown in fig. 5 shows near excellent agreement except for a minor difference in
the build up of the longitudinal electric field, see fig. 5.c. However, the EC and FB solver shows absolute
agreement.
14
Figure 5: Longitudinal electric field (a, c) and plasma density (b, d) for 1D LWFA simulations using
π-PIC with EC and FB solver as well as Smilei with M4 solver. The black rectangles in (a) and (b)
indicate the enlarged area shown in (c) and (d) respectively.
Having shown satisfactory agreement in 1D we move on to full 3D simulations of LWFA, where
additionally injection and subsequent acceleration of particles is simulated. This is achieved via density
downramp injection, using a plasma profile that decreases to one-third of the peak density (ne = 4 ·
1018, cm−3) at the end of the upramp, see fig. 6(d). The simulation parameters for the laser are identical
to the 1D case, with the addition of laser focusing, characterized by a focus at the position of the
downramp with a spot size of 40λ. While the results in fig. 6 are qualitatively similar, the on-axis
accelerating fields vary in amplitude and bubble-length between π-PIC and Smilei, despite the wave
fronts being identical. Moreover, there is a noticeable difference in the amount of injected charge and the
resulting phase-space distribution of the accelerated beam.
The observed differences in particle injection suggest that the discrepancy between the solvers arises
from their distinct shape functions. Smilei employs a 3-point stencil, whereas π-PIC uses a CIC shape
function corresponding to a 1-point stencil. In addition, Smilei enforces exact charge conservation through
Esirkepov’s deposition scheme [27]. In contrast, the solvers in π-PIC do not implement a similar deposition
scheme; however, the FBs solver applies a Poisson-based charge correction. The agreement in plasma
response between the EC and FB solvers indicates that charge non-conservation does not account for the
differences observed between the π-PIC and Smilei solvers, pointing instead to the shape function as the
most probable source of the discrepancy. Since the FB and energy-conserving solvers are once again in
complete agreement, we have excluded FB from the subsequent simulations.
15
Figure 6: Panel (a-c) shows the on-axis density at three different positions in the plasma marked in panel
(d). The upper half of (a-c) shows the π-PIC simulation and the lower the Smilei simulation, additionally
the on-axis electric field is shown in for π-PIC using the EC solver (full black) and FB solver (dashed
orange) as well as the M4 solver of Smilei (green dashed). Panel (e-f) shows the phase-space at the
longitudinal position corresponding to subplot (c).
To look further into these differences, we perform several simulations with different spatial and tem-
poral resolutions, as shown in fig. 7. The timestep was ∆t = ∆x/(c4) and ∆y = ∆x = 4λ. At high
resolution (∆x = λ/32), Smilei and π-PIC show strong agreement at the front of the wake, but a quali-
tative difference emerges at the rear end of the wake. Additionally, Smilei converges to a solution faster
with resolution, as the electric field for ∆x = λ/32 and ∆x = λ/16 is practically identical – which is not
the case for π-PIC.
For lower resolution there is a decrease in impact of the laser for both codes, shown at the front of
the wake. This is to be expected, since the ponderomotive force which drives the wake, is an averaged
force and by under-sampling the oscillations of the field the particles experience a lower effective push.
For π-PIC there is also an effective decrease of plasma wavelength for lower resolutions. Smilei, on the
other hand, shows a decrease in propagation speed of the wake which becomes progressively worse with
resolution. At a resolution of nx/λ = 4 (red line) the head of the wake is almost out of window.
In conclusion, we find that π-PIC retains a higher degree of accuracy at low resolution but the
convergence of the code with higher resolution is weak compared to Smilei.
16
Figure 7: Shows the longitudinal electric field at the plasma upramp corresponding to position (a) in
fig. 6 for simulations of π-PIC (upper) and Smilei (lower) for different resolutions.
6 Conclusion
PIC methods provide powerful tools for studying plasmas across a wide range of environments and scales.
The expansion of computing power makes previously unfeasible simulations achievable. Despite signif-
icant progress and ongoing research, challenges remain regarding numerical stability and conservation
properties. To accommodate further research, we have developed a flexible PIC framework, which en-
ables users to modify and extend the PIC scheme. Its modular structure supports parallel development
on multiple levels, such as solver and extension implementation, with the aim of making it a valuable
resource for the broader plasma community.
17
References
[1] H. Burau, R. Widera, W. Hönig, G. Juckeland, A. Debus, T. Kluge, U. Schramm, T. E. Cowan,
R. Sauerbrey, and M. Bussmann, “PIConGPU: A fully relativistic particle-in-cell code for a GPU
cluster,” IEEE Transactions on Plasma Science, vol. 38, no. 10, pp. 2831–2839, 2010.
[2] J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plot-
nikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech, “Smilei : A collaborative, open-source, multi-
purpose particle-in-cell code for plasma simulation,” Comput. Phys. Commun., vol. 222, pp. 351–373,
2018.
[3] R. A. Fonseca, L. O. Silva, F. S. Tsung, V. K. Decyk, W. Lu, C. Ren, W. B. Mori, S. Deng, S. Lee,
T. Katsouleas, and J. C. Adam, “OSIRIS: A three-dimensional, fully relativistic particle in cell code
for modeling plasma based accelerators,” in Computational Science — ICCS 2002 (P. M. A. Sloot,
A. G. Hoekstra, C. J. K. Tan, and J. J. Dongarra, eds.), pp. 342–351, Springer, 2002.
[4] S. Bastrakov, R. Donchenko, A. Gonoskov, E. Efimenko, A. Malyshev, I. Meyerov, and I. Surmin,
“Particle-in-cell plasma simulation on heterogeneous cluster systems,” Journal of Computational
Science, vol. 3, no. 6, pp. 474–479, 2012.
[5] J. L. Vay, A. Almgren, J. Bell, L. Ge, D. P. Grote, M. Hogan, O. Kononenko, R. Lehe, A. Myers,
C. Ng, J. Park, R. Ryne, O. Shapoval, M. Thévenet, and W. Zhang, “Warp-x: A new exascale
computing platform for beam–plasma simulations,” Nuclear Instruments and Methods in Physics
Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 909,
pp. 476–479, 2018.
[6] R. Lehe, M. Kirchen, I. A. Andriyash, B. B. Godfrey, and J.-L. Vay, “A spectral, quasi-cylindrical
and dispersion-free particle-in-cell algorithm,” Comput. Phys. Commun., vol. 203, pp. 66–82, 2016.
[7] OpenFOAM Foundation, “OpenFOAM: The Open Source CFD Toolbox.” https://www.
openfoam.com/. Accessed: 2025-08-20.
[8] Atomic Simulation Environment (ASE), “ASE: Atomic Simulation Environment.” https://wiki.
fysik.dtu.dk/ase/. Accessed: 2025-08-20.
[9] LAMMPS Development Team, “LAMMPS Molecular Dynamics Simulator.” https://www.
lammps.org/. Accessed: 2025-08-20.
[10] BLASWorking Group, “Basic Linear Algebra Subprograms (BLAS).” http://www.netlib.org/
blas/. Accessed: 2025-08-20.
[11] LAPACK Project, “LAPACK: Linear Algebra PACKage.” http://www.netlib.org/lapack/.
Accessed: 2025-08-20.
[12] Frigo, Matteo and Johnson, Steven G., “FFTW: Fastest Fourier Transform in the West.” http:
//www.fftw.org/. Accessed: 2025-08-20.
[13] COIN-OR Foundation, “COIN-OR: Computational Infrastructure for Operations Research.”
https://www.coin-or.org/. Accessed: 2025-08-20.
[14] Pedregosa, Fabian and Varoquaux, Gaël and Gramfort, Alexandre and Michel, Vincent and others,
“Scikit-learn: Machine Learning in Python.” https://scikit-learn.org/. Accessed: 2025-
08-20.
[15] A. Gonoskov, “Explicit energy-conserving modification of relativistic pic method,” J. Comput. Phys.,
vol. 502, p. 112820, Apr. 2024.
[16] K. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in
isotropic media,” IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp. 302–307,
1996.
[17] A. Pukhov, “Three-dimensional electromagnetic relativistic particle-in-cell code VLPL (virtual laser
plasma lab),” Journal of Plasma Physics, vol. 61, no. 3, pp. 425–433, 1999.
18
[18] M. Kärkkäinen, E. Gjonaj, T. Lau, and T. Weiland, “Low-dispersion wake field calculation tools,”
Proc. International Computational Accelerator Physics Conference, 01 2006.
[19] J. Cole, “A high-accuracy realization of the yee algorithm using non-standard finite differences,”
IEEE Transactions on Microwave Theory and Techniques, vol. 45, no. 6, pp. 991–996, 1997.
[20] R. Lehe, A. Lifschitz, C. Thaury, V. Malka, and X. Davoine, “Numerical growth of emittance in
simulations of laser-wakefield acceleration,” Physical Review Special Topics - Accelerators and Beams,
vol. 16, no. 2, p. 021301, 2013.
[21] B. M. Cowan, D. L. Bruhwiler, J. R. Cary, E. Cormier-Michel, and C. G. R. Geddes, “Generalized
algorithm for control of numerical dispersion in explicit time-domain electromagnetic simulations,”
Physical Review Special Topics - Accelerators and Beams, vol. 16, no. 4, p. 041303, 2013.
[22] C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulation. CRC Press, 2018.
[23] J.-L. Vay, I. Haber, and B. B. Godfrey, “A domain decomposition method for pseudo-spectral elec-
tromagnetic simulations of plasmas,” J. Comput. Phys., vol. 243, pp. 260–268, 2013.
[24] B. Marder, “A method for incorporating gauss’ law into electromagnetic PIC codes,” J. Comput.
Phys., vol. 68, no. 1, pp. 48–55, 1987.
[25] C. D. Munz, P. Omnes, R. Schneider, E. Sonnendrücker, and U. Voß, “Divergence correction tech-
niques for maxwell solvers based on a hyperbolic model,” J. Comput. Phys., vol. 161, no. 2, pp. 484–
511, 2000.
[26] J. Villasenor and O. Buneman, “Rigorous charge conservation for local electromagnetic field solvers,”
Comput. Phys. Commun., vol. 69, no. 2, pp. 306–316, 1992.
[27] T. Z. Esirkepov, “Exact charge conservation scheme for particle-in-cell simulation with an arbitrary
form-factor,” Comput. Phys. Commun., vol. 135, no. 2, pp. 144–153, 2001.
[28] T. Umeda, Y. Omura, T. Tominaga, and H. Matsumoto, “A new charge conservation method in
electromagnetic particle-in-cell simulations,” Comput. Phys. Commun., vol. 156, no. 1, pp. 73–85,
2003.
[29] X. Kong, M. C. Huang, C. Ren, and V. K. Decyk, “Particle-in-cell simulations with charge-conserving
current deposition on graphic processing units,” J. Comput. Phys., vol. 230, no. 4, pp. 1676–1685,
2011.
[30] K. Steiniger, R. Widera, S. Bastrakov, M. Bussmann, S. Chandrasekaran, B. Hernandez, K. Holsap-
ple, A. Huebl, G. Juckeland, J. Kelling, M. Leinhauser, R. Pausch, D. Rogers, U. Schramm, J. Young,
and A. Debus, “EZ: An efficient, charge conserving current deposition algorithm for electromagnetic
particle-in-cell simulations,” Comput. Phys. Commun., vol. 291, p. 108849, 2023.
[31] A. Pukhov, “Particle-in-cell codes for plasma-based particle acceleration,” CERN Yellow Reports,
vol. 1, pp. 181–181, 2016.
[32] A. T. Powis and I. D. Kaganovich, “Accuracy of the explicit energy-conserving particle-in-cell method
for under-resolved simulations of capacitively coupled plasma discharges,” Phys. Plasmas, vol. 31,
no. 2, p. 023901, 2024.
[33] A. B. Langdon, “Effects of the spatial grid in simulation plasmas,” J. Comput. Phys., vol. 6, no. 2,
pp. 247–267, 1970.
[34] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles, ch. 5-5, p. 149. CRC
Press, 2021.
[35] L. F. Ricketson and J. Hu, “An explicit, energy-conserving particle-in-cell scheme,” J. Comput.
Phys., vol. 537, p. 114098, Sept. 2025.
[36] J. Xiao, H. Qin, and J. Liu, “Structure-preserving geometric particle-in-cell methods for vlasov-
maxwell systems,” Plasma Science and Technology, vol. 20, no. 11, p. 110501, 2018.
19
[37] J. Squire, H. Qin, and W. M. Tang, “Geometric integration of the vlasov-maxwell system with a
variational particle-in-cell scheme,” Phys. Plasmas, vol. 19, no. 8, p. 084501, 2012.
[38] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker, “GEMPIC: geometric electromagnetic
particle-in-cell methods,” Journal of Plasma Physics, vol. 83, no. 4, p. 905830401, 2017.
[39] K. Kormann and E. Sonnendrücker, “Energy-conserving time propagation for a structure-preserving
particle-in-cell vlasov–maxwell solver,” J. Comput. Phys., vol. 425, p. 109890, 2021.
[40] A. Gonoskov, S. Bastrakov, E. Efimenko, A. Ilderton, M. Marklund, I. Meyerov, A. Muraviev,
A. Sergeev, I. Surmin, and E. Wallin, “Extended particle-in-cell schemes for physics in ultrastrong
laser fields: Review and developments,” Physical review E, vol. 92, no. 2, p. 023305, 2015.
[41] V. Volokitin, J. Magnusson, A. Bashinov, E. Efimenko, A. Muraviev, and I. Meyerov, “Optimized
event generator for strong-field qed simulations within the hi-chi framework,” Journal of Computa-
tional Science, vol. 74, p. 102170, 2023.
[42] A. Gonoskov, “Agnostic conservative down-sampling for optimizing statistical representations and
PIC simulations,” Comput. Phys. Commun., vol. 271, p. 108200, Feb. 2022.
[43] A. Muraviev, A. Bashinov, E. Efimenko, V. Volokitin, I. Meyerov, and A. Gonoskov, “Strategies for
particle resampling in PIC simulations,” Comput. Phys. Commun., vol. 262, p. 107826, May 2021.
[44] W. Jakob, “pybind11: Seamless interoperability between C++ and Python.” https://pybind11.
readthedocs.io/en/stable/. Accessed: 2025-09-30.
[45] hi chi, “π-pic: A python-controlled interactive particle-in-cell code.” https://github.com/
hi-chi/pipic/tree/pipic_tutorial, 2023. [Code].
[46] J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput.
Phys., vol. 114, no. 2, pp. 185–200, 1994.
[47] T. Tajima and Y. C. Lee, “Absorbing boundary condition and budden turning point technique for
electromagnetic plasma simulations,” J. Comput. Phys., vol. 42, no. 2, pp. 406–412, 1981.
[48] G. A. Mourou, T. Tajima, and S. V. Bulanov, “Optics in the relativistic regime,” Reviews of Modern
Physics, vol. 78, p. 309–371, Apr. 2006.
[49] M. Marklund and P. K. Shukla, “Nonlinear collective effects in photon-photon and photon-plasma
interactions,” Reviews of Modern Physics, vol. 78, p. 591–640, May 2006.
[50] A. Di Piazza, C. Müller, K. Z. Hatsagortsyan, and C. H. Keitel, “Extremely high-intensity laser
interactions with fundamental quantum systems,” Reviews of Modern Physics, vol. 84, p. 1177–1228,
Aug. 2012.
[51] F. Krausz and M. Ivanov, “Attosecond physics,” Reviews of Modern Physics, vol. 81, p. 163–234,
Feb. 2009.
[52] L. W. Davis, “Theory of electromagnetic beams,” Physical Review A, vol. 19, p. 1177–1179, Mar.
1979.
[53] M. Couture and P.-A. Belanger, “From gaussian beam to complex-source-point spherical wave,”
Physical Review A, vol. 24, p. 355–359, July 1981.
[54] A. M. Fedotov, K. Y. Korolev, and M. V. Legkov, “Exact analytical expression for the electromag-
netic field in a focused laser beam or pulse,” in ICONO 2007: Physics of Intense and Superintense
Laser Fields; Attosecond Pulses; Quantum and Atomic Optics; and Engineering of Quantum Infor-
mation (M. V. Fedorov, W. Sandner, E. Giacobino, S. Kilin, S. Kulik, A. Sergienko, A. Bandrauk,
and A. M. Sergeev, eds.), vol. 6726, p. 672613, SPIE, June 2007.
[55] S. M. Sepke and D. P. Umstadter, “Analytical solutions for the electromagnetic fields of tightly
focused laser beams of arbitrary pulse length,” Optics Letters, vol. 31, p. 2589, Aug. 2006.
[56] Y. Salamin, “Fields of a gaussian beam beyond the paraxial approximation,” Applied Physics B,
vol. 86, p. 319–326, Sept. 2006.
20
[57] K. I. Popov, V. Y. Bychenkov, W. Rozmus, and R. D. Sydora, “Electron vacuum acceleration by a
tightly focused laser pulse,” Phys. Plasmas, vol. 15, Jan. 2008.
[58] D. E. Cardenas, T. M. Ostermayr, L. Di Lucchio, L. Hofmann, M. F. Kling, P. Gibbon, J. Schreiber,
and L. Veisz, “Sub-cycle dynamics in relativistic nanoplasma acceleration,” Scientific Reports, vol. 9,
May 2019.
[59] A. A. Gonoskov, A. V. Korzhimanov, A. V. Kim, M. Marklund, and A. M. Sergeev, “Ultrarelativistic
nanoplasmonics as a route towards extreme-intensity attosecond pulses,” Physical Review E, vol. 84,
Oct. 2011.
[60] C. Harvey, M. Marklund, and A. R. Holkundkar, “Focusing effects in laser-electron thomson scat-
tering,” Physical Review Accelerators and Beams, vol. 19, Sept. 2016.
[61] E. Panova, V. Volokitin, E. Efimenko, J. Ferri, T. Blackburn, M. Marklund, A. Muschet, A. De An-
dres Gonzalez, P. Fischer, L. Veisz, et al., “Optimized computation of tight focusing of short pulses
using mapping to periodic space,” Applied Sciences, vol. 11, no. 3, p. 956, 2021.
[62] Y. Lu, P. Kilian, F. Guo, H. Li, and E. Liang, “Time-step dependent force interpolation scheme for
suppressing numerical cherenkov instability in relativistic particle-in-cell simulations,” J. Comput.
Phys., vol. 413, p. 109388, 2020.
[63] T. Esirkepov, “Exact charge conservation scheme for particle-in-cell simulation with an arbitrary
form-factor,” Comput. Phys. Commun., vol. 135, pp. 144–153, Apr. 2001.
[64] D. N. Arnold, R. S. Falk, and R. Winther, “Finite element exterior calculus, homological techniques,
and applications,” Acta Numerica, vol. 15, p. 1–155, May 2006.
[65] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker, “Gempic: geometric electromagnetic
particle-in-cell methods,” Journal of Plasma Physics, vol. 83, July 2017.
[66] T. Tajima and J. M. Dawson, “Laser electron accelerator,” Physical Review Letters, vol. 43, no. 4,
pp. 267–270, 1979.
[67] E. Esarey, C. B. Schroeder, and W. P. Leemans, “Physics of laser-driven plasma-based electron
accelerators,” Reviews of Modern Physics, vol. 81, no. 3, pp. 1229–1285, 2009.
21

91
References
[1] E. Esarey, C. B. Schroeder, and W. P. Leemans, “Physics of laser-driven
plasma-based electron accelerators,” Rev. Mod. Phys., 81, no. 3, pp. 1229–
1285, Aug. 27, 2009. DOI: 10.1103/RevModPhys.81.1229.
[2] A. Hasegawa, “Introduction to plasma instabilities,” Plasma Insta-
bilities and Nonlinear Effects, A. Hasegawa, Ed., Berlin, Heidelberg:
Springer, 1975, pp. 1–26, ISBN: 978-3-642-65980-5. DOI: 10.1007/978-3-
642-65980-5_1.
[3] T. Tajima and J. M. Dawson, “Laser electron accelerator,” Physical
Review Letters, 43, no. 4, pp. 267–270, Jul. 23, 1979. DOI: 10 . 1103/
PhysRevLett.43.267.
[4] M. Galletti et al., “Prospects for free-electron lasers powered by plasma-
wakefield-accelerated beams,” Nat. Photon., 18, no. 8, pp. 780–791, Aug.
2024. DOI: 10.1038/s41566-024-01474-3.
[5] S. F. Martins, R. A. Fonseca, J. Vieira, L. O. Silva, W. Lu, and W. B. Mori,
“Modeling laser wakefield accelerator experiments with ultrafast particle-
in-cell simulations in boosted frames,” Phys. Plasmas, 17, no. 5, p. 056 705,
Mar. 29, 2010. DOI: 10.1063/1.3358139.
[6] J. L. Vay et al., “Modeling of 10 GeV-1 TeV laser-plasma accelerators us-
ing lorentz boosted simulations,” Phys. Plasmas, 18, no. 12, p. 123 103,
Dec. 13, 2011. DOI: 10.1063/1.3663841.
[7] B. B. Godfrey, “Numerical cherenkov instabilities in electromagnetic
particle codes,” J. Comput. Phys., 15, no. 4, pp. 504–521, Aug. 1, 1974.
DOI: 10.1016/0021-9991(74)90076-X.
[8] D. A. Gurnett, Introduction to plasma physics : with space and labo-
ratory applications, in collab. with Internet Archive. Cambridge, UK ;
New York : Cambridge University Press, 2005, 470 pp., ISBN: 978-0-521-
36483-6 978-0-521-36730-1.
[9] A. A. Vlasov, “THE VIBRATIONAL PROPERTIES OF AN ELECTRON
GAS,” Sov. Phys. Usp., 10, no. 6, p. 721, Jun. 30, 1968. DOI: 10.1070/
PU1968v010n06ABEH003709.
[10] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion. Boston,
MA: Springer US, 1984, ISBN: 978-1-4419-3201-3 978-1-4757-5595-4.
DOI: 10.1007/978-1-4757-5595-4.
92 RESEARCH PAPERS
[11] D. B. Melrose, “Kinetic plasma physics,” Plasma Astrophysics, J. G. Kirk,
D. B. Melrose, E. R. Priest, A. O. Benz, and T. J.-L. Courvoisier, Eds.,
Berlin, Heidelberg: Springer, 1994, pp. 113–223, ISBN: 978-3-540-31627-
5. DOI: 10.1007/3-540-31627-2_2.
[12] A. Picksley et al., “Matched guiding and controlled injection in dark-
current-free, 10-GeV-class, channel-guided laser-plasma accelerators,”
Phys. Rev. Lett., 133, no. 25, p. 255 001, Dec. 18, 2024. DOI: 10.1103/
PhysRevLett.133.255001.
[13] A. J. Gonsalves et al., “Petawatt laser guiding and electron beam accel-
eration to 8 GeV in a laser-heated capillary discharge waveguide,” Phys.
Rev. Lett., 122, no. 8, p. 084 801, Feb. 25, 2019. DOI: 10.1103/PhysRevLett.
122.084801.
[14] E. Esarey, A. Ting, P. Sprangle, D. Umstadter, and X. Liu, “Nonlinear
analysis of relativistic harmonic generation by intense lasers in plas-
mas,” IEEE Trans. Plasma Sci., 21, no. 1, pp. 95–104, Feb. 1993. DOI:
10.1109/27.221107.
[15] B. A. Shadwick, C. B. Schroeder, and E. Esarey, “Nonlinear laser energy
depletion in laser-plasma accelerators,” Phys. Plasmas, 16, no. 5, May 1,
2009. DOI: 10.1063/1.3124185.
[16] S. Bulanov, N. Naumova, F. Pegoraro, and J. Sakai, “Particle injection
into the wave acceleration phase due to nonlinear wake wave break-
ing,” Phys. Rev. E, 58, no. 5, R5257–R5260, Nov. 1, 1998. DOI: 10.1103/
PhysRevE.58.R5257.
[17] M. Chen, E. Esarey, C. B. Schroeder, C. G. R. Geddes, and W. P. Leemans,
“Theory of ionization-induced trapping in laser-plasma accelerators,”
Phys. Plasmas, 19, no. 3, p. 033 101, Mar. 2012. DOI: 10.1063/1.3689922.
[18] T. C. Katsouleas, S. Wilks, P. Chen, J. M. Dawson, and J. J. Su, “Beam
loading in plasma accelerators,” Part. Accel., 22, pp. 81–99, 1987.
[19] M. Kirchen et al., “Optimal beam loading in a laser-plasma accelera-
tor,” Phys. Rev. Lett., 126, no. 17, p. 174 801, Apr. 2021. DOI: 10.1103/
PhysRevLett.126.174801.
[20] K. Jensen et al., “Improved laser-plasma accelerator stability via high-
bandwidth longitudinal focal position stabilization of a 100 TW-class
laser system,” Phys. Rev. Accel. Beams, 28, no. 9, p. 092 802, Sep. 26,
2025. DOI: 10.1103/j5md-sckl.
93
[21] A. R. Maier et al., “Decoding sources of energy variability in a laser-
plasma accelerator,” Phys. Rev. X, 10, no. 3, p. 031 039, Aug. 18, 2020.
DOI: 10.1103/PhysRevX.10.031039.
[22] L. Rovige et al., “Demonstration of stable long-term operation of a kilo-
hertz laser-plasma accelerator,” Phys. Rev. Accel. Beams, 23, p. 093 401,
9 Sep. 2020. DOI: 10.1103/PhysRevAccelBeams.23.093401.
[23] A. Sanchez-Gonzalez et al., “Accurate prediction of x-ray pulse proper-
ties from a free-electron laser using machine learning,” Nat. Commun,
8, no. 1, p. 15 461, Jun. 5, 2017. DOI: 10.1038/ncomms15461.
[24] R. J. Adcock, “A problem in least squares,” The Analyst, 5, no. 2, pp. 53–
54, 1878. DOI: 10.2307/2635758.
[25] C. H. Kummell, “Reduction of observation equations which contain
more than one observed quantity,” The Analyst, 6, no. 4, pp. 97–105,
1879, ISSN: 0741-7918. DOI: 10.2307/2635646.
[26] J. A. Hausman, W. K. Newey, H. Ichimura, and J. L. Powell, “Identi-
fication and estimation of polynomial errors-in-variables models,”
J. Econom., 50, no. 3, pp. 273–295, Dec. 1, 1991. DOI: 10.1016/0304-
4076(91)90022-6.
[27] S. Schennach, “Estimation of nonlinear models with measurement
error,” Econometrica, 72, p. 33, 2004. DOI: https://doi.org/10.1111/j.
1468-0262.2004.00477.x.
[28] L. Leng, T. Zhang, L. Kleinman, and W. Zhu, “Ordinary least square
regression, orthogonal regression, geometric mean regression and
their applications in aerosol science,” J. Phys. Conf. Ser., 78, p. 012 084,
2007. DOI: https://doi.org/10.1088/1742-6596/78/1/012084.
[29] Z. H., Z. T. K., and J. S. L., “Measurement Error Models: From Non-
parametric Methods to Deep Neural Networks,” Statistical Science, 37,
p. 473, 2022. DOI: https://doi.org/10.1214/21-STS834.
[30] C. Molnar, “8.5 permutation feature importance,” Interpretable Ma-
chine Learning. Leanpub, 2020. [Online]. Available: https://christophm.
github.io/interpretable-ml-book/feature-importance.html.
[31] C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simu-
lation. Boca Raton: CRC Press, Oct. 8, 2018, 504 pp., ISBN: 978-1-315-
27504-8. DOI: 10.1201/9781315275048.
[32] M. Pohl, M. Hoshino, and J. Niemiec, “PIC simulation methods for
cosmic radiation and plasma instabilities,” Prog. Part. Nucl. Phys., 111,
p. 103 751, Mar. 1, 2020. DOI: 10.1016/j.ppnp.2019.103751.
94 RESEARCH PAPERS
[33] A. Gonoskov et al., “Extended particle-in-cell schemes for physics in
ultrastrong laser fields: Review and developments,” Phys. Rev. E, 92,
no. 2, p. 023 305, Aug. 18, 2015. DOI: 10.1103/PhysRevE.92.023305.
[34] J.-L. Vay, “Noninvariance of space- and time-scale ranges under a
lorentz transformation and the implications for the study of relativistic
interactions,” Phys. Rev. Lett., 98, no. 13, p. 130 405, Mar. 30, 2007. DOI:
10.1103/PhysRevLett.98.130405.
[35] S. F. Martins, R. A. Fonseca, L. O. Silva, W. Lu, and W. B. Mori, “Numerical
simulations of laser wakefield accelerators in optimal lorentz frames,”
Comput. Phys. Commun., 181, no. 5, pp. 869–875, May 1, 2010. DOI:
10.1016/j.cpc.2009.12.023.
[36] J. L. Vay, C. G. R. Geddes, E. Cormier-Michel, and D. P. Grote, “Nu-
merical methods for instability mitigation in the modeling of laser
wakefield accelerators in a lorentz-boosted frame,” J. Comput. Phys.,
230, no. 15, pp. 5908–5929, Jul. 1, 2011. DOI: 10.1016/j.jcp.2011.04.003.
[37] P. Yu et al., “Elimination of the numerical cerenkov instability for spec-
tral EM-PIC codes,” Comput. Phys. Commun., 192, pp. 32–47, Jul. 1,
2015. DOI: 10.1016/j.cpc.2015.02.018.
[38] P. Yu et al., “Enabling lorentz boosted frame particle-in-cell simulations
of laser wakefield acceleration in quasi-3d geometry,” J. Comput. Phys.,
316, pp. 747–759, Jul. 1, 2016. DOI: 10.1016/j.jcp.2016.04.014.
[39] B. B. Godfrey, J.-L. Vay, and I. Haber, “Numerical stability analysis of
the pseudo-spectral analytical time-domain PIC algorithm,” J. Comput.
Phys., 258, pp. 689–704, Feb. 1, 2014. DOI: 10.1016/j.jcp.2013.10.053.
[40] M. Kirchen, R. Lehe, S. Jalas, O. Shapoval, J.-L. Vay, and A. R. Maier,
“Scalable spectral solver in galilean coordinates for eliminating the nu-
merical cherenkov instability in particle-in-cell simulations of stream-
ing plasmas,” Phys. Rev. E, 102, no. 1, p. 013 202, Jul. 7, 2020. DOI: 10.
1103/PhysRevE.102.013202.
[41] I. Haber, R. Lee, H. Klein, and J. Boris, “Advances in electromagnetic
simulation techniques,” Proc. Sixth Conf. Num. Sim. Plasmas, pp. 46–
48, 1973.
[42] “Proceedings of the conference on the numerical simulation of plas-
mas (4th) held at the naval research laboratory, washington, d.c. on
2, 3 november 1970. ”[Online]. Available: https://apps.dtic.mil/sti/
citations/ADA023511.
95
[43] A. Gonoskov, “Explicit energy-conserving modification of relativistic
PIC method,” J. Comput. Phys., 502, p. 112 820, Apr. 1, 2024. DOI: 10.
1016/j.jcp.2024.112820.
[44] M. Palmroth et al., “Vlasov methods in space physics and astrophysics,”
Living Reviews in Computational Astrophysics, 4, no. 1, pp. 1–54, Dec. 1,
2018. DOI: 10.1007/s41115-018-0003-2.
[45] A. Pukhov, “Particle-in-cell codes for plasma-based particle acceler-
ation,” CERN Yellow Reports, 181 Pages, Feb. 16, 2016. DOI: 10.5170/
CERN-2016-001.181.
[46] A. T. Powis and I. D. Kaganovich, “Accuracy of the explicit energy-
conserving particle-in-cell method for under-resolved simulations
of capacitively coupled plasma discharges,” Phys. Plasmas, 31, no. 2,
p. 023 901, Feb. 1, 2024. DOI: 10.1063/5.0174168.
[47] A. Langdon, “Effects of the spatial grid in simulation plasmas,” J. Com-
put. Phys., 6, no. 2, pp. 247–267, Oct. 1970. DOI: 10.1016/0021-9991(70)
90024-0.
[48] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Parti-
cles. Boca Raton: CRC Press, Mar. 24, 2021, 540 pp., ISBN: 978-0-367-
80693-4. DOI: 10.1201/9780367806934.
[49] B. McMillan, “Is it necessary to resolve the debye length in standard or
δ f PIC codes?” Phys. Plasmas, 27, no. 5, 2020. DOI: 10.1063/1.5139957.
[50] J. U. Brackbill, “On energy and momentum conservation in particle-
in-cell plasma simulation,” J. Comput. Phys., 317, pp. 405–427, Jul. 15,
2016. DOI: 10.1016/j.jcp.2016.04.050.
[51] D. C. Barnes and L. Chacón, “Finite spatial-grid effects in energy-
conserving particle-in-cell algorithms,” Comput. Phys. Commun., 258,
p. 107 560, Jan. 1, 2021. DOI: 10.1016/j.cpc.2020.107560.
[52] H. R. Lewis, “Energy-conserving numerical approximations for vlasov
plasmas,” J. Comput. Phys., 6, no. 1, pp. 136–141, Aug. 1, 1970. DOI:
10.1016/0021-9991(70)90012-4.
[53] A. B. Langdon, ““energy-conserving” plasma simulation algorithms,” J.
Comput. Phys., 12, no. 2, pp. 247–268, Feb. 1, 1973. DOI: 10.1016/S0021-
9991(73)80014-2.
[54] S. Markidis and G. Lapenta, “The energy conserving particle-in-cell
method,” J. Comput. Phys., 230, no. 18, pp. 7037–7052, Aug. 1, 2011.
DOI: 10.1016/j.jcp.2011.05.033.
96 RESEARCH PAPERS
[55] J. R. Angus et al., “An implicit particle code with exact energy and
charge conservation for electromagnetic studies of dense plasmas,”
J. Comput. Phys., 491, p. 112 383, Oct. 2023. DOI: 10.1016/j.jcp.2023.
112383.
[56] G. Lapenta, “Exactly energy conserving semi-implicit particle in cell
formulation,” J. Comput. Phys., 334, pp. 349–366, Apr. 1, 2017. DOI:
10.1016/j.jcp.2017.01.002.
[57] J. Derouillat et al., “Smilei : A collaborative, open-source, multi-purpose
particle-in-cell code for plasma simulation,” Comput. Phys. Commun.,
222, pp. 351–373, Jan. 1, 2018. DOI: 10.1016/j.cpc.2017.09.024.
[58] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker, “GEM-
PIC: Geometric electromagnetic particle-in-cell methods,” Journal
of Plasma Physics, 83, no. 4, p. 905 830 401, Aug. 2017. DOI: 10.1017/
S002237781700040X.
[59] J. Villasenor and O. Buneman, “Rigorous charge conservation for lo-
cal electromagnetic field solvers,” Comput. Phys. Commun., 69, no. 2,
pp. 306–316, Mar. 1, 1992. DOI: 10.1016/0010-4655(92)90169-Y.
[60] T. Z. Esirkepov, “Exact charge conservation scheme for particle-in-cell
simulation with an arbitrary form-factor,” Comput. Phys. Commun.,
135, no. 2, pp. 144–153, Apr. 1, 2001. DOI: 10.1016/S0010- 4655(00)
00228-9.
[61] T. Umeda, Y. Omura, T. Tominaga, and H. Matsumoto, “A new charge
conservation method in electromagnetic particle-in-cell simulations,”
Comput. Phys. Commun., 156, no. 1, pp. 73–85, Dec. 1, 2003. DOI: 10.
1016/S0010-4655(03)00437-5.
[62] X. Kong, M. C. Huang, C. Ren, and V. K. Decyk, “Particle-in-cell simula-
tions with charge-conserving current deposition on graphic processing
units,” J. Comput. Phys., 230, no. 4, pp. 1676–1685, Feb. 20, 2011. DOI:
10.1016/j.jcp.2010.11.032.
[63] K. Steiniger et al., “EZ: An efficient, charge conserving current deposi-
tion algorithm for electromagnetic particle-in-cell simulations,” Com-
put. Phys. Commun., 291, p. 108 849, Oct. 1, 2023. DOI: 10.1016/j.cpc.
2023.108849.
[64] B. Marder, “A method for incorporating gauss’ law into electromagnetic
PIC codes,” J. Comput. Phys., 68, no. 1, pp. 48–55, Jan. 1, 1987. DOI:
10.1016/0021-9991(87)90043-X.
97
[65] C. .-. Munz, P. Omnes, R. Schneider, E. Sonnendrücker, and U. Voß,
“Divergence correction techniques for maxwell solvers based on a
hyperbolic model,” J. Comput. Phys., 161, no. 2, pp. 484–511, Jul. 1,
2000. DOI: 10.1006/jcph.2000.6507.
[66] A. Blinne, D. Schinkel, S. Kuschel, N. Elkina, S. G. Rykovanov, and
M. Zepf, “A systematic approach to numerical dispersion in maxwell
solvers,” Comput. Phys. Commun., 224, pp. 273–281, Mar. 1, 2018. DOI:
10.1016/j.cpc.2017.10.010.
[67] K. Yee, “Numerical solution of initial boundary value problems in-
volving maxwell’s equations in isotropic media,” IEEE Transactions
on Antennas and Propagation, 14, no. 3, pp. 302–307, May 1966. DOI:
10.1109/TAP.1966.1138693.
[68] A. Pukhov, “Three-dimensional electromagnetic relativistic particle-
in-cell code VLPL (virtual laser plasma lab),” Journal of Plasma Physics,
61, no. 3, pp. 425–433, Apr. 1999. DOI: 10.1017/S0022377899007515.
[69] M. Kärkkäinen, E. Gjonaj, T. Lau, and T. Weiland, “Low-dispersion wake
field calculation tools,” Proc. International Computational Accelerator
Physics Conference, Jan. 2006.
[70] J. Cole, “A high-accuracy realization of the yee algorithm using non-
standard finite differences,” IEEE Transactions on Microwave Theory
and Techniques, 45, no. 6, pp. 991–996, Jun. 1997. DOI: 10.1109/22.
588615.
[71] R. Lehe, A. Lifschitz, C. Thaury, V. Malka, and X. Davoine, “Numerical
growth of emittance in simulations of laser-wakefield acceleration,”
Phys. Rev. Accel. Beams, 16, no. 2, p. 021 301, Feb. 28, 2013. DOI: 10.
1103/PhysRevSTAB.16.021301.
[72] B. M. Cowan, D. L. Bruhwiler, J. R. Cary, E. Cormier-Michel, and C. G. R.
Geddes, “Generalized algorithm for control of numerical dispersion in
explicit time-domain electromagnetic simulations,” Phys. Rev. Accel.
Beams, 16, no. 4, p. 041 303, Apr. 4, 2013. DOI: 10.1103/PhysRevSTAB.
16.041303.
[73] X. Xu et al., “Numerical instability due to relativistic plasma drift in
EM-PIC simulations,” Comput. Phys. Commun., 184, no. 11, pp. 2503–
2514, Nov. 1, 2013. DOI: 10.1016/j.cpc.2013.07.003.
[74] J.-P. Berenger, “A perfectly matched layer for the absorption of elec-
tromagnetic waves,” J. Comput. Phys., 114, no. 2, pp. 185–200, Oct. 1,
1994. DOI: 10.1006/jcph.1994.1159.
98 RESEARCH PAPERS
[75] T. Tajima and Y. C. Lee, “Absorbing boundary condition and budden
turning point technique for electromagnetic plasma simulations,” J.
Comput. Phys., 42, no. 2, pp. 406–412, Aug. 1, 1981. DOI: 10.1016/0021-
9991(81)90254-0.
[76] R. Lehe, M. Kirchen, B. B. Godfrey, A. R. Maier, and J.-L. Vay, “Elimina-
tion of numerical cherenkov instability in flowing-plasma particle-in-
cell simulations by using galilean coordinates,” Phys. Rev. E, 94, no. 5,
p. 053 305, Nov. 14, 2016. DOI: 10.1103/PhysRevE.94.053305.
[77] M. Kirchen et al., “Stable discrete representation of relativistically drift-
ing plasmas,” Phys. Plasmas, 23, no. 10, p. 100 704, Oct. 10, 2016. DOI:
10.1063/1.4964770.
[78] O. Shapoval, E. Zoni, R. Lehe, M. Thévenet, and J.-L. Vay, “Pseudospec-
tral particle-in-cell formulation with arbitrary charge and current-
density time dependencies for the modeling of relativistic plasmas,”
Phys. Rev. E, 110, no. 2, p. 025 206, Aug. 23, 2024. DOI: 10.1103/PhysRevE.
110.025206.
APPENDIX
A Derivation of numerical dispersion relation
In this section the algebraic manipulations to go from eq. (4.25) to eq. (4.30)
is shown as well as the calculations for eq. (4.31). Rewriting the right hand
side of eq. (4.25) and evaluating the derivative
‚ Œ

2 ω ∂ Ek +εk l mεmnp pl kn Ep/γmω ω p p γ ∂ p i
=
′
k ω −pj k ′ j /γm p⃗=(p0,0,0)

1 v 2E 1 ” —
= δ 0 x′ 0i + Ei + ′ εi 0mεmnpγkn v0Ep/ω +γ(ω′− v 2 ′0kx ) c γ(ω − v0kx )
 
v 2 Ep 
+δ 0i 0 (−2ω′+ v k ′0 x )γEk + v
′
0kkγ(Ek +εk 0mεmnp k v )
 ,
c 2 n 0 ω 
v⃗=(v0,0,0)
where the last term is without summation on i . Considering only two dimen-
sions with Ez = 0 and extracting i = 0, 1 corresponding to x , y
i = 0 : (ω2− c 2k 2)Ex + c 2kx (kx Ex +ky Ey ) =
 
ω2pω ” —
= ω′/γ2 + v 2 ′ ′
γ(ω′− v k ′ )2 0
ky ky /ω Ex + v0ky (1− v0kx /ω)Ey
0 x
i = 1 : (ω2− c 2k 2)E 2y + c ky (kx Ex +ky Ey ) =
ω2  pω v0
γ(ω′− v k ′
Ey + (ky E −k) ω x x
Ey ) .
0 x
100 DERIVATION OF NUMERICAL DISPERSION RELATION
Gathering all Ex and Ey terms on either side
– 2 ‚ 2 ′ Œ™ωpω v0 ky ky
i = 0 : E 2 2 2 2 2 ′ 2x ω − c k + c kx − ω /γ + =γ(ω′− v0k ′ 2x ) ω
– 2 2 ′ ™ ωpωv0 ky v k
= E 2y −c kx k
0 x
y + 1− →γ(ω′− v0k ′x )2 ω
– 2 ™ω
E ω2− c 2k 2 + c 2 2 p ′ 2 2 ′x kx − ′ (ωω /γ + v0 kγ(ω′− v k )2 y
ky ) =
0 x
–
ω2 v k ′
™
2 p 0− y= Ey c kx ky + (ω− v k ) ,γ(ω′− v0k ′x )2
0 x
– ™
ω2p v0ky
i = 1 : Ex c
2ky kx − =γ(ω′− v ′0kx )
– ™
ω2
= E −ω2 + c 2k 2− c 2k 2 py y + γ(ω′− v k ′
(ω− v
) 0
kx ) .
0 x
To summarize, the components in the matix equation E E⃗ = 0 becomes
–
ω2
™
ε = ω2− c 2 211 ky −
p
(ωω′/γ2 + v 2k k ′ )
γ(ω′− v k ′ )2 0 y y0 x
–
ω2 v k ′
™
p 0 y
ε12 = c
2kx ky − (ω− v k )γ(ω′− v k ′ )2 0 x0 x
–
ω2
™
p v0ky
ε21 = c
2ky kx − γ(ω′− v ′0kx )
– 2 ™ω
ε = ω2− c 2k 2 − p22 x γ(ω′− v k ′
(ω− v0k) x
)
0 x
101
The dispersion relation is given by ε11ε22−ε12ε21 = 0. Gathering terms with
factors (ω2p/γ)
0, (ω2 /γ)1p , (ω
2
p/γ)
2 we get
(ω2 /γ)0 2 2 2 2p :ω (ω − c k )
2 1 1
€ Š
(ω /γ) : − (ω2− c 2k 2 )(ω− v k ) + c 2k 2p ω′− v k ′ y 0 x y
kx v0 +
0 x
1 € Š
′ − (ωω
′/γ2 + v 2k k ′ )(ω2− c 2k 2) + (ω− v k )v c 2k k ′ k
(ω′− v k )2 0 y y x 0 x 0 y y x0 x
2 2 1 ω
′ω
(ωp/γ) : (ω− v k )(ω′− v0k ′ )3x γ2
0 x
Multiplying all lines with (ω′− v k ′ 20 x ) /ω
2
(ω2 /γ)0 :(ω′− v k ′ )2(ω2− c 2k 2p 0 x )
ω′− v k ′ € Š
( 0ω2p/γ)
1 : x − (ω2− c 2k 2y )(ω− v0k ) + c
2k 2
ω2 x y
kx v0 +
1 € Š
− (ωω′/γ2 + v 20 ky k
′
y )(ω
2− c 2k 2x ) + (ω− v0k
2 ′
ω2 x
)v0c ky ky kx =
ω′
=− (ω2− c 2k 2)− (ω′− v0k ′x )(ω− vωγ2 0
kx )+
ω′ c
2k 2
− c 2 yk 2y + (ω
′− v k ′0 x )+ωγ2 ω
k ′ 2y v0c kx ky v
2k ′y k
+ (ω− 0 yv0k )− (ω2x − c 2k 2)ω2 ω2 x
ω′ ω− v k
( 0 xω2p/γ)
2 :
ωγ2 ω′− v ′0kx
The red colored terms can be gathered as
ω4ω′p ω− v k ω
2
p ω
2ω′
0 x − (ω′− v k ′ p′ 0 x )(ω− v0kx )− (ω
2− c 2k 2)
γ4ω ω′− v k γ γ30 x ω
+ (ω′− v k ′ )20 x (ω
2− c 2k 2) =
ω2 ′ ω2€ Š€ Š
= (ω′− v k ′ p
ω
)2− ω2 2 2 p
ω− v0k− c k − x0 x γ3 ω γ ω′− v k ′0 x
102 DERIVATION OF NUMERICAL DISPERSION RELATION
The remaining terms are
ω2
 ′ c 2k 2
g (ω,ω′
p ω
, k⃗ , k⃗ ′) = − c 2k 2 y+ (ω′− v k ′
γ ωγ2 y ω 0 x
)+
k ′ v c 2 2 ′

y 0 kx ky v0 ky ky
+ (ω−v k)− (ω
2− c 2k
2
0 x x ) =ω2 ω2

ω2 ′ 2

p ω € 1 Š2 2 v c k= c ky 1−
0 y
+ (k ′y kx −k
′ 2
x ky )− v0 k k
′
γ ω γ2 ω y y
=
2  ωp €2 ω
′ v c 2Š
− ′ 0
ky
= v0 ky ky ky + (k
′
y kx −k
′
x ky ) .γ ω ω