# Modeling of a Liquid Leaf Target TNSA Experiment using Particle-In-Cell
Simulations and Deep Learning

B. Schmitz<EMAIL_ADDRESS>Technische Universität Darmstadt
Institut für Teilchenbeschleunigung und Elektromagnetische Felder (TEMF),
Schlossgartenstr.8, 64289 Darmstadt, Germany  D. Kreuter  University of
Cambridge
Department of Applied Mathematics and Theoretical Physics (DAMTP), Centre for
Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
Technische Universität Darmstadt
Institut für Teilchenbeschleunigung und Elektromagnetische Felder (TEMF),
Schlossgartenstr.8, 64289 Darmstadt, Germany  O. Boine-Frankenheim  Technische
Universität Darmstadt
Institut für Teilchenbeschleunigung und Elektromagnetische Felder (TEMF),
Schlossgartenstr.8, 64289 Darmstadt, Germany  GSI Helmholtzzentrum für
Schwerionenforschung GmbH, Planckstr. 1, 64291 Darmstadt, Germany

(August 28, 2024)

###### Abstract

Liquid leaf targets show promise as high repetition rate targets for laser-
based ion acceleration using the Target Normal Sheath Acceleration (TNSA)
mechanism and are currently under development. In this work, we discuss the
effects of different ion species and investigate how they can be leveraged for
use as a possible laser-driven neutron source. To aid in this research, we
develop a surrogate model for liquid leaf target laser-ion acceleration
experiments, based on artificial neural networks. The model is trained using
data from Particle-In-Cell (PIC) simulations. The fast inference speed of our
deep learning model allows us to optimize experimental parameters for maximum
ion energy and laser-energy conversion efficiency. An analysis of parameter
influence on our model output, using Sobol’ and PAWN indices, provides deeper
insights into the laser-plasma system.

TNSA, Deep Learning, PIC, Liquid Leaf, Lorentz Boost, Multi-Species, Numerical
Optimization

††preprint: APS/123-QED

##  I Introduction

Laser-accelerated ions have great potential for various applications, such as
compact medical accelerators [1, 2, 3, 4], neutron sources [1, 5, 6, 7] or as
injectors for conventional accelerators[8]. These applications require a high
repetition rate to overcome the drawback of the exponential energy
distribution, typical for Target Normal Sheath Acceleration. However,
conventional solid-state targets cannot achieve high repetition rates due to
engineering difficulties and target supply[9, chapter 4.2]. For this reason,
different targets such as gas [10] or liquid-based targets [11, 12] are
currently being developed.

In this work, we investigate a liquid leaf target [13] currently under
development at TU Darmstadt. This target is a major step towards achieving
reproducible, high repetition rate ion bunches from laser-plasma interactions,
which is necessary for any kind of application. This new system allows for the
operation of a repetitive target with arbitrary H2O/D2O ratios for the first
time, which we investigate in this work.

In particular, we aim to train a surrogate model for a liquid leaf target in a
target normal sheath acceleration (TNSA) experiment to understand the
characteristics of the liquid leaf and its composition, predict ideal
operating points and understand how multiple ion species interact with each
other. The first aim for the target at TU Darmstadt is the creation of a
viable compact neutron source, which requires proton energies larger than the
production threshold of neutrons ($>$1.7\text{\,}\mathrm{MeV}$$)[6].

Previous attempts at modeling laser-plasma acceleration experiments have been
made [14, 15, 16, 17]. With our contributions, presented in this work, we
expand on the state-of-the-art by taking more parameters of the experimental
setup into account and providing a surrogate for the theory of intricate
effects a multi-species target can have on the energy spectrum of the
accelerated ions.

Huebl et al. [10] have found a strong influence of the mixture ratio of
multiple species on the resulting energy spectra in hydrogen-deuterium
targets. In subsection II.1, we expand on their idea and provide indicators
for the importance of this effect.

Furthermore, while attempts at modeling a laser-plasma acceleration experiment
using neural networks have been made [14], we extend previous work with a vast
number of Particle-In-Cell simulations to train our models. The chosen
approach via deep learning also ensures that expansion (transfer learning) of
the model with experimental data is possible. We demonstrate our surrogate
model’s high performance and utility by optimizing an example laser-plasma
acceleration experiment, leveraging non-trivial relationships between the
experimental parameters not yet understood by theory (see section III).

##  II Plasma target models

The following section details our contributions related to the considered
multi-species target experiment as well as the creation of our simulation
datasets and the training of our surrogate model.

For this work, we carried out various PIC simulations to generate our
datasets. The bulk of the simulations was computed on the Virgo High-
Performance Computing cluster [18] at GSI Helmholtzzentrum, Darmstadt. For
these simulations, we used the Smilei [19] PIC code. We determined the
resulting surrogate by training an artificial neural network from the
simulated data.

###  II.1 Multi-species target considerations

We are considering a liquid leaf target which consists of multiple different
atom species. These ion species differ in their charge-to-mass ratio
$q_{i}/m_{i}$. and an ion can have up to $Z$ different charge states. Taking
water, for example, one can have up to 8 ionization states of oxygen and an
additional one for the hydrogen component.

Water occurs naturally with different isotopes of hydrogen. Taking into
account regular water (H2O) and heavy water (D2O), an additional degree of
freedom—the mixtures between the two—must be considered. Since several ion
species are present, this can be denoted as a $n$-species plasma, where $n$ is
the number of ion states in the plasma.

The final non-relativistic kinetic energy of species $i$ accelerated in a
constant electric field $E_{0}$ scales as

$E_{\text{kin}}\propto E_{0}^{2}\frac{q_{i}^{2}}{m_{i}}\;.$ (1)

Therefore, species with a higher $q_{i}^{2}/m_{i}$ ratio will gain more
energy. Provided that the initial densities are similar, species with higher
$q_{i}^{2}/m_{i}$ will deplete most of the available field energy. This energy
is then split between different particle species and limits the acceleration
efficiency of a single species.

Several species interact with each other, leading to a deformation of the
particle spectrum. Faster particles take electrons from the sheath and screen
the acceleration field for the following heavier particles. These heavier
particles are then accelerated in the screened field and hence have less
kinetic energy per nucleon and a lower velocity than expected from the
assumption above. Mid-energy lighter particles are accelerated by the
following heavier particle front due to the Coulomb force, getting compressed
in the momentum space. This compression causes plateaus and quasi-
monoenergetic features to form.

This effect is described in detail and analytically calculated for the
asymptotic case for two particle species (deuterium and hydrogen gas) by Huebl
et al. [10]. We applied their solutions for 2 species, since we assumed a
fully ionized plasma in our simulations to reduce the degrees of freedom
inside the plasma.

The compression effect on the lighter particle spectrum is visualized in
Figure 1 using a PIC simulation for regular H2O.

Figure 1: Particle energy spectra of hydrogen and oxygen after TNSA PIC
simulation of liquid leaf water target. The dotted lines are the corresponding
Mora [20] fits for the displayed spectrum. Large deviations from the spectrum
($30\text{\,}\mathrm{MeV}$ and up) can be explained by the multi-species
effect. The simulation setup is described in subsection II.2. The investigated
features are sharp and their shape varies. One dimensional simulations have a
sharper profile, while higher dimensional ones and real-life experiments are
smoother [10, 21].

Deviations from the ideal Mora[20, 10] can be seen.

We can make two observations: Firstly, the lower energy part of the spectrum,
in this case until half of the maximum energy, is coarser than the
corresponding higher energy part of the spectrum. There is also a peak at
around $10\text{\,}\mathrm{MeV}$ in both spectra which makes it possible to
compare the spectra against each other. These peaks are shifted by the same
amount as the oxygen cut-off is shifted from the hydrogen plateau, which can
be assumed to be a correlation due to the particle interaction in this energy
range.
Secondly, there is a plateau in the hydrogen spectrum starting at around
$30\text{\,}\mathrm{MeV}$. This plateau and the corresponding dip before it
deviate fairly strongly from the established Mora theory for TNSA (indicated
by the dotted lines). This drop/increase combination is explainable by the
previously introduced multi-species effect and we want to investigate,
predict, and leverage this behavior.

If we can describe and predict this effect, we can optimize our ion beam for
specific applications. To do this, we need to find a surrogate model for the
full spectrum problem. Further insights gained from considering multiple
species become evident in subsubsection II.2.6. Fully ionized oxygen has the
same $q_{i}/m_{i}$ ratio as deuterium, for example, which reduces the
efficiency of deuterium acceleration. In this work, we only modeled the proton
part of the spectrum because our data is ambiguous for the Oxygen/Deuterium
combination part. However, expanding the dataset to include a sweep of the
oxygen’s charge number would resolve this ambiguity and yield clearer modeling
results for deuterium.

###  II.2 Particle-In-Cell Simulations Setup

The simulations reflect a real experiment in reduced dimensions. To sample a
larger parameter space in a reasonable time, we reduced the dimensions of the
simulation to 1.5D. This means simulating one space and three momentum
components. The fields are also sampled in three dimensions.

We further applied an additional method to account for angle dependency by
applying a transverse Lorentz boost to the system. Details on both the Lorentz
boost and the method itself can be found in Appendix C. A sketch of the full
setup is displayed in Figure 2.

Figure 2: Overview of the simulation setup. Green marks the plasma target. The
lighter green areas indicate the pre-plasma and the skirt implemented. The
laser, indicated by the red arrow, hits the plasma under an angle $\Theta$ –
relative to the target normal. After the acceleration time, the momenta of the
accelerated particles, given in blue, are registered. For the liquid leaf
target, $d_{1}$ is assumed to be equal to $d_{3}$.

####  II.2.1 Plasma target

The target in the conducted simulations models a liquid leaf target under
development at TU Darmstadt Institute of Nuclear Physics and which is similar
to the work by George et al. [13]. The liquid leaf target’s width is some cm,
while the typical irradiation size of a laser is in the order of
$\mathrm{\SIUnitSymbolMicro m}$. We assume, that the surface roughness is
negligible and the plasma surface is therefore considered to be planar.

When the target is only dependent on one coordinate, it can be described fully
by its particle density profile. Thus, the simulation only allows movement
along the $x$ coordinate and is independent of $y$ and $z$. We also assume
that the plasma is expanded when the main pulse hits the target, and that the
pre-plasma and skirt follow an exponential profile. We chose the scale length
for the exponential profile as $0.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ so
as to be longer than a comparable setup with a cryogenic (i.e. less
evaporative) jet target [22] while still ensuring a well-defined plasma
border. The exponential profile thus takes the shape

$\displaystyle n_{\text{exp}}(x)$
$\displaystyle=\frac{n_{0}}{1+\exp(-(x-x_{\text{front}})/l_{\text{s}})}$ (2)

where $l_{\text{s}}=$0.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$x_{\text{front}}$ is the location of the target front. The skirt has
identical functional shape for the backside of the target.

Since a liquid leaf target evaporates, we superimposed the typical vapor
density distribution for a liquid leaf target, given by

$\displaystyle n_{\text{LLT}}(r)$ $\displaystyle\approx
n\left(r_{\text{jet}}\right)\left(\frac{r_{\text{jet}}}{r}\right)^{2}\frac{L_{\text{L}}}{\sqrt{r^{2}+L_{\text{L}}^{2}}}\;,$
(3)

where $n\left(r_{\text{jet}}\right)$ is the water vapor density at the liquid
jet surface, $L_{\text{L}}\approx$3\text{\,}\mathrm{cm}$$ is the liquid jet
length, and $r_{\text{jet}}$ is the liquid jet radius [23]. Note that the
second term in the above expression has been squared as we expect a faster
drop-off of the liquid leaf density in our proposed experimental setup.

The assumed particle densities are
$n_{0}=$6.68\text{\times}{10}^{28}\text{\,}{\mathrm{m}}^{-3}$$ and
$n\left(r_{\text{jet}}\right)=$1.62\text{\times}{10}^{23}\text{\,}{\mathrm{m}}^{-3}$$
stemming from the liquid density of water and the density estimated at the
saturation vapor pressure at $0\text{\,}\mathrm{\SIUnitSymbolCelsius}$ [23].
We also introduced a cut-off of the profile
$4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ before and after the target, which
washes out by approximately $1.4\text{\,}\mathrm{\SIUnitSymbolMicro m}$ by the
time the laser hits the target. This cut-off is only introduced to optimize
the simulation’s performance.

We chose to investigate multi-species effects resulting from a combination of
different ion species inside the target. We simulated regular water, heavy
water, and a potential mixture between the two. This mixture is indicated by
the mixture parameter listed in Table 1, which we varied in discrete steps.
The simulation thus consists of up to four species: electrons (e-), hydrogen
(H+), deuterium (D+), and oxygen (On+). As the ionization of oxygen is of
importance to the model, this was varied as well. All particle species follow
the same distribution function defined above.

The ion species are initialized cold while the electrons received an initial
temperature of $30\text{\,}\mathrm{keV}$ to simulate interaction with a pre-
pulse. We used Smilei’s defaults for particle initialization, including no
ionization or radiation model [24]. The length of one cell is the Debye length
at the initial electron temperature, around $5\text{\,}\mathrm{nm}$. For the
time resolution, a CFL number of 0.98 was used. The interpolation order of the
particle shape functions is set to four and the particle per cell count for
each species is 800.

####  II.2.2 Laserpulse

In a 1.5D simulation, the laser pulse is given by its time profile only. We
assumed a Gaussian time profile, using Smilei’s tgaussian profile with the
following shape:

$\displaystyle I_{\text{envelope}}(t)$
$\displaystyle=\begin{cases}\exp\left(\frac{-(t-\tau_{\text{L}})^{2}}{(\tau_{\text{L}}/2)^{2}/\ln(2)}\right)&\text{if
}t\leq 2\tau_{\text{L}}\\\ 0&\text{otherwise}\end{cases}\;,$ (4)

where $\tau_{\text{L}}$ is the main laser pulse duration. In this work, we
deal with lasers that have a pulse duration
$\tau_{\text{L}}<$1\text{\,}\mathrm{ps}$$ and an $a_{0}>1$. The laser energy
$E_{\text{L}}$, pulse length $\tau_{\text{L}}$, polarization, incident angle
$\theta_{\text{L}}$, wavelength $\lambda_{\text{L}}$ and the target thickness
$d_{\text{T}}$ are variable and are uniformly sampled from the defined
intervals in Table 1. Our thought process in choosing exactly these parameters
was that we needed to cover the full system, which required the use of 9
parameters. These parameters were chosen based on two different, sometimes
contradictory paradigms: one was to allow the similarity equations to take
full effect, while the other was to enable experimental validation of the
model (see also Appendix A and Appendix B)

Table 1: Table of the physical parameters that were used for sampling of the input files to the 1.5D PIC simulations. Mixture defines the percentage of hydrogen substituted by deuterium. No |  | Attribute | Sign | Range | Units
---|---|---|---|---|---
1 | Laser | Energy | $E_{\text{L}}$ | [0.001, 50] | $\mathrm{J}$
2 | Laser | Focus-FWHM | FWHM | [2,20] | $\mathrm{\SIUnitSymbolMicro m}$
3 | Laser | Pulse length | $\tau_{\text{L}}$ | [15, 150] | $\mathrm{fs}$
4 | Laser | Polarization |  | {s, p} |
5 | Laser | Incidence angle | $\theta_{\text{L}}$ | [0, 85] | $\mathrm{\SIUnitSymbolDegree}$
6 | Laser | Wavelength | $\lambda_{\text{L}}$ | [550, 1100] | $\mathrm{nm}$
7 | Target | Thickness | $d_{\text{T}}$ | [0.6, 3] | $\mathrm{\SIUnitSymbolMicro m}$
8 | Target | Mixture | Mix | [0, 100] | $\mathrm{\char 37\relax}$
9 | Target | Oxygen Charge | $Z_{\text{eff}}$ | {7, 8} |

####  II.2.3 Simulation Output Quantities

The diagnostics recorded are the particles’ $x$-coordinate, all components of
the momentum $\vec{p}$, and the macro-particle weight $w$ at the acceleration
time

$\displaystyle t_{\text{acc}}=\tau_{\text{L}}+d_{\text{T}}/c_{\text{s}}\;,$
(5)

where $c_{\text{s}}$ is the ion-acoustic velocity. Lécz [25] has found that
this is a suitable acceleration time after which an isothermal plasma
expansion model no longer holds. From these recorded values we reconstruct the
energy spectrum of the particles in the lab frame by using Eq. (34). Since all
energy spectra have an individual shape and cut-off energy, the spectra were
each normalized to the energy range $[0,1]$, counted into 100 bins, and stored
as a list together with their respective cut-off energies. In order to keep
the numbers more practical, we took the logarithm. An entry for the results of
a simulation thus has the following shape:
$\\{\ln\left(\frac{\mathrm{d}n}{\mathrm{d}E}\right)_{\text{Bin
1}},\ldots,\ln\left(\frac{\mathrm{d}n}{\mathrm{d}E}\right)_{\text{Bin
100}},E_{\text{max}}\\}$. Exponentiating and re-scaling by $E_{\text{max}}$
restores the original energy spectrum accordingly. This same recording scheme
is used for all four species for all simulations.

We chose that the parameters in Table 1 are uniformly sampled with exception
of the laser energy $E_{\text{L}}$ which we sampled following a square root
scale and the mixture was varied in discrete steps. This type of sampling
results in significantly more simulations with low $a_{0}$ than with high
$a_{0}$. To deal with this we forced additional simulations onto dedicated
intervals of $a_{0}$. Although the laser focus-FWHM is technically not
relevant in the 1D case we sampled it nonetheless such that together with the
sampled laser energy and pulse length the correct $a_{0}$ was written in the
input file. This also ensures comparability with higher-order simulations and
experimental data.

####  II.2.4 Simulation statistics

We used the setup described above to create a dataset of simulations for our
subsequent surrogate model. All parameters were stochastically sampled and
their combination can be thought of as a sparse grid. The Virgo cluster
employs the Simple Linux Utility for Resource Management (SLURM) [26] to
schedule incoming jobs where up to $10\,000$ jobs can be added to the queue
simultaneously. The jobs were queued using a script to sample a certain number
of parameter combinations and then start a simulation job for each of them.
The number of simulations varies between the different species. There were
$508\,200$ simulations for hydrogen and $762\,426$ simulations for deuterium,
resulting in a total of $1\,270\,626$ simulations. However, the precise number
of simulations is not crucial, as long as the number of simulations is in a
similar order of magnitude, the results should be comparable. The reduced
model, which utilizes only the pure H2O data without D2O component was trained
on a subset of the full dataset with $68\,973$ entries accordingly.

####  II.2.5 Limitations of 1.5D PIC

We used 1.5D simulations as mentioned earlier. These low-dimensional
simulations do have some drawbacks. While they, together with our introduced
transversal Lorentz boost method (Appendix C), are capable of describing
several effects, some are not possible. The main limitation is created by the
expansion of the plasma behind the target. In one spatial dimension, no
transversal drift of the particles is possible, therefore also no decay of
space charge effects exists. The expansion continues until infinity if it is
not stopped. Even though we introduced an effective acceleration time
$t_{\text{acc}}$, this problem persists. Since we keep both setup and method
constant, the relative behavior of the cut-off energies can still be taken
into account, but the absolute value is overestimated. This overestimation is
predictable and when applied makes the models directly comparable. Lécz et al.
[27] have shown that the acceleration time $t_{\text{acc}}$ cuts off the
spectrum, such that it is a good approximation of 2D simulations. The
simulations have been verified with experiments as well, which have shown that
the bias can be mitigated. Furthermore, Sinigardi et al. [28] have shown
further scalings between 2D and 3D cutoff energies. Taking both arguments into
account we can deduce, that there is a constant scaling factor from 1D to
real-world experiments and also to 3D simulations. Similarly, because of the
lack of transversal particle movement, we cannot evaluate divergence opening
angles in a 1D simulation.

####  II.2.6 Data Discussion by Example

We display an example of the spatial distribution from the simulations in
Figure 3. An example of the energy spectrum is already displayed in Figure 1.

Firstly, in this simulation, we assumed that the target consists of regular
water and is fully ionized by the implied laser pre-pulse. Thus, the three
species (e-, H+, and O8+) are initialized with a density ratio of $10:2:1$
respectively such that overall neutrality is conserved. We display the
species’ positions at $t=t_{\text{acc}}$ in Figure 3.

Figure 3:  Example PIC simulation of water leaf target TNSA experiment. The
plot shows the particle distribution at the previously proposed acceleration
time $t_{\text{acc}}$ (Eq. (5)).

In this simulation the laser incidence angle is
$0\text{\,}\mathrm{\SIUnitSymbolDegree}$, the target thickness is
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and the dimensionless laser
amplitude is $a_{0}=20$. We observe the two ion species, H+ and O8+ at
$t=t_{\text{acc}}$. The figures show that the species have different positions
at the measured time, which means that the species are accelerated separately
by the sheath field.

The ion front position at the acceleration time for different species varies
due to the different charge and mass values as mentioned in subsection II.1.
Calculating the expected variation, following the relation from Huebl et al.,
for only fully ionized oxygen and hydrogen present, yields a scaling factor of
$x_{F}^{O^{8+}}/x_{F}^{H}\approx 0.68$. The corresponding factor from Figure 3
is $x_{F,\text{sim}}^{O^{8+}}/x_{F,\text{sim}}^{H}\approx$0.67(1)$$, where the
uncertainty results from the binning.

We can see that the general TNSA mechanism is still applicable. Although the
dynamics of the different particle species with each other are more complex,
as we will see later, the general behavior appears to follow classical TNSA
theory. This is supported by the kinetic energy spectra of the ion species
after acceleration, an example of which is shown in Figure 1. The figure shows
the energy spectra of hydrogen and oxygen ions, along with Mora’s predicted
ideal curve.

###  II.3 Deep Learning Application

We have to correlate the different simulations to each other and find
relations and interpolations to allow for an optimization of the full setup.
We decided to use a neural network approach with fully connected feedforward
topologies and built them in Keras [29] running inside of Tensorflow 2 [30].
For hyperparameter tuning, we used the Keras Tuner module [31].

####  II.3.1 Model Training

To predict a particle spectrum, two models are needed. The spectrum model
continuously maps $\\{\text{[physical
parameters]}\\}=\\{E,\text{mix},E_{\text{L}},r_{\text{L}},\tau_{\text{L}},\text{s/p-pol.},\theta_{\text{L}},\lambda_{\text{L}},d_{\text{T}}\\}$
onto $\ln\left(\frac{\mathrm{d}n}{\mathrm{d}E}(E)\right)$ while a second
cutoff model only predicts the maximum energy (i.e. when to cut off the
continuous spectrum from the first model). We trained a reduced model pair,
not taking deuterons into account for regular H2O, and a full model pair
containing different ratios between H2O and D2O. The dedicated features of the
PIC simulation can be seen better with the reduced model. We assume that this
is a result of the lower number of input dimensions and therefore of the
deviating degrees of generalization.

We essentially think of the energy spectrum as the graph of a continuous
function $f$, the first model maps $\\{x,\text{[system parameters]}\\}$ to
$f(x)$ while the second model predicts the point $x$ at which the graph gets
cut off. Details about the training parameters and the procedure are given in
Appendix E. The reduced spectrum model has 6 hidden layers $(x\rightarrow
320\rightarrow 288\rightarrow 288\rightarrow 256\rightarrow 256\rightarrow
320\rightarrow 1)$ while the full spectrum model has 11 hidden layers with 460
neurons each. The cutoff models both have 8 hidden layers $(x\rightarrow
320\rightarrow 284\rightarrow 288\rightarrow 512\rightarrow 32\rightarrow
480\rightarrow 512\rightarrow 32\rightarrow 1)$. It is worth noting that the
input dimension of the reduced spectrum model is one less than the full
spectrum model since it does not include the mix parameter. All networks were
fully-connected architectures with ReLU activations on their hidden layers. We
will now briefly discuss and evaluate the trained models:

##### Reduced Model Pair:

The precision of the cutoff models, which attempt to map $\\{\text{[physical
parameters]}\\}$ onto $E_{\text{max}}$ can be estimated rather easily. For the
reduced problem, the model achieved a mean squared error of
$8.93\text{\,}{\mathrm{MeV}}^{2}$ on validation data (confer to appendix E),
meaning the average error on the prediction of the hydrogen spectrum’s maximum
energy is projected to be around $\pm 2.99\text{\,}\mathrm{MeV}$.

To more intuitively evaluate the reduced spectrum model’s predicting
capabilities and potential shortcomings, ten simulations with equal parameters
(except for random seed) were computed such that their hydrogen ion energy
spectra could be compared to the predicted spectrum of the model. A plot of
all the spectra is shown in Figure 4.

Figure 4: Reconstructed hydrogen ion energy spectra of ten simulations
differing only in their random seed. The energy spectrum prediction by the
trained neural network model is indicated with a red dashed line, while the
average of the simulations is indicated by the blue dotted line. The curve is
obtained from the reduced continuous model and is cut off at the maximum
energy determined by the maximum energy model.

The overall agreement of the model with the simulations is evident. The
maximum energy predicted by the cutoff model falls centrally between the
maximum energies of the ten simulations, only differing from the simulation
average by $0.2\text{\,}\mathrm{MeV}$. Looking at more intricate features of
the simulation spectra, however, it is clear that the model possibly
generalized slightly too much. At around $10\text{\,}\mathrm{MeV}$ a dip,
possibly due to multi-species effects, can be observed in most of the
simulations and yet is barely present in the model prediction. Generally, the
fluctuations in the simulation spectra are greatly reduced in the neural
network predicted spectrum. A reason for this is likely the sheer vastness of
differing spectra the model was trained on. Since the parameter space for the
training simulations was so large, the model had to generalize to many very
different output spectra.

##### Full Model Pair:

The full model pair were trained exactly the same as the reduced model pair
but with an additional parameter and a larger dataset. The full cutoff model
converged with a mean squared error of $7.25\text{\,}{\mathrm{MeV}}^{2}$ on
validation data resulting in a prediction error of $\pm
2.7\text{\,}\mathrm{MeV}$ for the maximum energy of the hydrogen spectrum
(confer to appendix E).

Again, as given above, the sensitivity of the spectrum model is more
complicated to estimate.

Figure 5: Model comparison for hydrogen spectra with reference PIC simulation.
Dashed lines give the result for the full model, the number indicates the
value for the mixture parameter. The PIC reference (for mixture = 0) is
displayed with a solid line and the reduced model with a dotted line.

In Figure 5, we can see that both numerical models, the full and the reduced
model do deviate from one another slightly, even if the mixture is set to
zero. This is expected behavior since there is a statistical variation in the
training of neural networks. Important to note is the deviation in the spectra
for different mixture ratios. An influence of the mixture parameter on the
spectrum is visible and it can be used to tune the spectrum. The behavior for
the full spectrum model is the same as the one for the reduced spectrum model
presented in Figure 4, the model is generalizing to a specific degree and has
an uncertainty of few MeV for the cut-off energy.

####  II.3.2 Model Efficiency

Calling the models in a Python code environment is similar to calling any
other function and takes around $20\text{\,}\mathrm{ms}$ on a personal laptop.
This time is in stark contrast to the four hours on 16 CPUs taken to run a
similar 1D PIC simulation on the HPC cluster. To put this in perspective, we
can run inference on the models roughly $720\,000$ times, while one PIC
simulation calculates.

We made other attempts at fitting the regression problem using various kernel
combinations and Gaussian Process Regression [32], however, they never
produced energy spectrum predictions that even came close to the neural
network prediction seen in Figure 4, usually being off from simulations by
orders of magnitude. As expected, the adaptability of modern unparameterized
machine learning methods such as neural network models stands out from other
regressors.

##  III Application of the Model

With a trained surrogate in hand, we were able to take advantage of the model
to perform a numerical optimization of an experiment as well as evaluate our
models’ interpretability.

###  III.1 Optimization of Parameters for Laser Plasma Interactions

In this section, we optimize a TNSA experiment with a water leaf target. We
aim at an ideal set of laser and target parameters and apply the previously
obtained reduced machine learning model pair.

We chose some base parameters that have a proven repetition rate of at least
$1\text{\,}\mathrm{Hz}$: Ti:Sa lasers with a central wavelength of
$800\text{\,}\mathrm{nm}$ and p-polarized laser light. Exemplary systems would
be the VEGA-3 laser at the Centro De Laseres Pulsados in Salamanca, Spain
(CLPU) [33] or DRACO laser at the Helmholt-Zentrum Dresden-Rossendorf [34].
Following the procedure in this section, the model can also be applied to any
other system, if its parameters are inside the minimal and maximal physical
parameters of our model (see Table 1). If the system’s parameters are not
included, our model could be expanded by retraining with additional data,
using transfer learning [35], or other modern domain adaptation methods [36].
The assumed initial parameters of the laser system are stated in Table 2.

Table 2: Table of the physical parameters to be optimized for the laser system. Both initial and optimized values are shown. Rows in bold remained fixed during optimization. The dimensionless laser amplitude $a_{0}$ also remained fixed during optimization to encourage the convergence towards non-trivial parameter combinations. $\eta_{\text{conv}}$ is a measure for energy conversion efficiency (see Eq. (6)), normalized to the initial parameter case. No | Attribute | Initial | Optimized | Optimized | Units
---|---|---|---|---|---
|  |  | ($E_{\text{max}}$) | ($E\text{-conversion}$) |
1 | Laser energy | 30 | 6.6 | 1.4 | $\mathrm{J}$
2 | Focus-FWHM | 20 | 4.2 | 2 | $\mathrm{\SIUnitSymbolMicro m}$
3 | Pulse length | 30 | 149.9 | 137.6 | $\mathrm{fs}$
4 | Polarization | p | p | p |
5 | Incidence angle | 12.2 | 32.2 | 29.3 | $\mathrm{\SIUnitSymbolDegree}$
6 | Wavelength | 800 | 800 | 800 | nm
7 | Thickness | 2 | 3.0 | 3.0 | $\mathrm{\SIUnitSymbolMicro m}$
| $E_{\text{max}}$ | 13.8 | 51.5 | 51.2 | $\mathrm{MeV}$
| $\eta_{\text{conv}}$ | 1.0 | 7.8 | 41.3 |

In this section, we investigate two different optimization goals. The first
goal is to find the maximum cut-off energy, while the second goal is to
maximize the laser energy deposition into the plasma.

As mentioned we assumed polarization and central laser wavelength as fixed but
otherwise allowed all parameters to change, as long as they stayed in the
given physical constraints. Since the obvious solution to maximizing output
energy is to maximize input energy, the optimizations were computed under the
constraint of a constant dimensionless laser amplitude $a_{0}$. This ensures
optimization by exploiting complicated relationships between the physical
parameters of the system; a task that can only be feasibly solved with a
rapidly callable model.

We implemented the optimization utilizing the SciPy Python library [37] and
the Byrd-Omojokun algorithm [38] included in its scipy.optimize.minimize
routine. The Byrd-Omojokun algorithm allows us to include both, boundary
conditions according to Table 1 as well as the equality constraint of constant
$a_{0}$, to leverage the aforementioned non-trivialities of the system.

The optimized parameters are displayed in Table 2. The optimizer seems to have
taken advantage of incidence angle-dependent absorption effects such as
resonance absorption. Additionally, by dramatically increasing laser focus
while simultaneously decreasing laser power (energy over time) the maximum ion
energy could be optimized without changing the dimensionless laser amplitude
$a_{0}$. Overall, the optimizer was able to increase the maximum ion output
energy by a factor of roughly 4. The hydrogen energy spectra for these
optimized parameters as well as for the initial parameters are depicted in
Figure 6.

A more intricate measure of a TNSA experimental system is the laser-ion energy
conversion efficiency, i.e. the measure of how much of the laser’s input
energy gets transported into the accelerated particles. We thus consider the
optimization of the ratio of the total kinetic energy of the ions
$E_{\text{H}}$ to the laser pulse energy $E_{\text{L}}$:

$\displaystyle\operatorname*{arg\,max}_{x\in\\{\text{params}\\}}\frac{E_{\text{H}}(x)}{E_{\text{L}}}$
$\displaystyle=\operatorname*{arg\,max}_{x\in\\{\text{params}\\}}\frac{1}{E_{\text{L}}}\cdot\int_{0}^{E_{\text{max}}}\frac{\mathrm{d}N}{\mathrm{d}E}\cdot
E\,\mathrm{d}E\;,$ (6)

where $\frac{\mathrm{d}N}{\mathrm{d}E}(E,x)$ and $E_{\text{max}}$ are given by
the neural network models, and {params} is the set of all parameter
combinations within the ranges specified in Table 1. It is important to note
that the Smilei output gives $\frac{\mathrm{d}n}{\mathrm{d}E}$ which has to be
scaled by a unit volume $V$ to arrive at the expression needed. For further
explanation of how to arrive at the above integral term we refer to Appendix
D. Here, we also allow the variation of laser energy $E_{\text{L}}$,
increasing the complexity of the problem. The optimization described in Eq.
(6) was carried out by solving the numerical integral using the composite
trapezoidal rule and once again employing the Byrd-Omojokun algorithm.

(a) Initial parameter selection.

(b) Optimized parameter selection with respect to maximum ion energy.

Figure 6:  Energy spectra of H-ions for a TNSA water leaf target experiment
using both initial VEGA-3 as well as optimized parameters with respect to the
maximum ion energy. Predicted spectra by the neural network model and spectra
from a 1D PIC simulation are shown. The parameters are given in Table 2.

As seen in Table 2, despite having a slightly lower maximum ion energy than
the first optimization task, the calculated energy conversion efficiency is
more than five times greater. This gives a strong indication that laser
coupling into the target in a laser-plasma experiment depends on the physical
parameters of the system in a highly non-trivial way.

###  III.2 Sensitivity Analysis

Artificial neural networks are generally difficult to interpret, which is a
drawback we have to accept. Nevertheless, the importance of specific
parameters for a model can be evaluated. One way to quantify the impact of a
model’s input parameters on its output is to use variance-based global
sensitivity analysis, also known as the Sobol’ method. The corresponding
sensitivity metrics are known as Sobol’ indices[39, 40, 41]. The Sobol’
indices are calculated by Monte Carlo sampling of parameters and corresponding
model outputs. This method is used to apportion the variance of the output to
the inputs and their combinations. The number of evaluations of our model is
$N\times(2D+2)$, where $D$ is the number of input features and $N$ is the
number of samples drawn. $N$ is ideally selected as a power of 2, where we
selected $2^{18}=$262\,144$$ drawn samples.

We used the PAWN method [42] for a second sensitivity analysis to complement
the Sobol’ method due to its shortcomings for the higher order of the input
features. The PAWN method uses a different approach for when the variance
might not be a good measure for the outcome of a system. It utilizes the
traits of the Cumulative Distribution Functions with similar Monte Carlo
sampling as for the Sobol’ indices, giving a different approach to determine
the sensitivity of a model. A combination of these two methods was also
proposed by Baroni et al. [43].

#### Reduced Model

The reduced cutoff model has 7 input features which are mapped to 1 output
prediction for the maximal energy. Our results of the Sobol’ analysis for the
reduced H2O-only model are given in Figure 7.

(a) Total variation which explains the cut-off energy variation.

(b) Matrix of dependencies to explain the cut-off energy variation. The
diagonal gives first-order Sobol’ indices, while the lower gives the second-
order Sobol’ indices for the corresponding variables. The upper line is the
numerical value, the lower line gives the corresponding error.

Figure 7: Sobol’ sensitivity analysis results showing the influence of various
physical parameters on the cut-off energy of H-ions for a TNSA water leaf
target experiment, for the reduced model utilizing only the H2O data. Errors
are given in the $95\text{\,}\mathrm{\char 37\relax}$ confidence level.

The larger the value of a Sobol’ index is, the more influence the independent
parameter has on the result. The total Sobol’ indices, normally referred to as
$S_{T}$ give a measure of the total importance of the given features. The
total Sobol’ indices can neither describe how much of the variance is
attributed to which combination of parameters nor are they normalized for the
total expression. This is due to multiple counting of effects, e.g. if there
is a second-order contribution for $\Theta_{L}$ and $r_{L}$, then this
contribution is added to both of the values in the total representation. It
doubles the counting for the second order, triples for the third order, and so
on.

Due to this complication, the determination of higher-order dependencies makes
it necessary to display the first and second-order Sobol’ indices, as done in
Figure 7 (b). The values are displayed in a matrix, such that the interaction
between $(x_{i},y_{j})$ can be displayed. The first order Sobol’ indices are
shown on the main diagonal ($x_{i}=y_{j}$). It is evident from the plot, that
the sum does not add up to 1, leaving approximately $21\text{\,}\mathrm{\char
37\relax}$ of data variance unexplained. The consequence of this is, that even
higher order interactions are necessary to fully explain the variation in our
model. Full calculation of higher orders has been omitted as it was deemed
unfeasible due to the extreme computational cost for higher dimensions.

The results of the PAWN method are displayed in Figure 8.

Figure 8: PAWN indices for the reduced model as a measure for parameter
importance. Boxes consist of uncertainty value, minimum, median, maximum, and
upper uncertainty value. The numerical value given is the median.

PAWN can only give us a measure of the full importance of the individual
parameters. A subsequent division into main effects and higher order is not
possible.

Table 3: Importance ranking of the model parameters as calculated by the Sobol’ and PAWN methods. Importance | Sobol’ | PAWN
---|---|---
1 | $r_{L}$ | $\Theta_{L}$
2 | $\Theta_{L}$ | $r_{L}$
3 | $E_{L}$ | $E_{L}$
4 | $\tau_{L}$ | $\tau_{L}$
5 | $\lambda_{L}$ | $\lambda_{L}$
6 | $d_{T}$ | $d_{T}$
7 | Pol | Pol

The importance ranking from PAWN does not entirely match the order found by
the Sobol’ method but is rather close. Both are listed in Table 3. If not the
total, but the sum of first and second-order Sobol’ is taken, then the first
two features change places.

The sensitivity analyses thus suggest that higher-order interactions are
important in this model and a simple optimization (e.g. maximizing only one
quantity) is not sufficient. Our previously presented optimizations take this
implicitly into account. Furthermore, the incidence angle and the irradiation
area appear to be important. The angle $\Theta_{\text{L}}$’s high influence is
expected, considering laser-ion absorption mechanisms, and is faithfully
implemented into the 1D simulation space using the Lorentz boosted geometry
(see Appendix C). While the third quantity, the laser energy, directly scales
the laser’s dimensionless amplitude $a_{0}$, the irradiation radius
$r_{\text{L}}$’s influence is more difficult to understand. The irradiation
area is not directly represented in a 1.5D PIC simulation. However, since
$a_{0}$ is dependent on $r_{\text{L}}$ an indirect influence is included.

#### Full Model

Running the same analysis for the full cutoff model including the mixture
parameters yields the results given in Figure 9 and Figure 10.

(a) Total variation which explains the cut-off energy variation.

(b) Matrix of dependencies to explain the cut-off energy variation. The
diagonal gives first-order Sobol’ indices, while the lower gives the second-
order Sobol’ indices for the corresponding variables. The upper line is the
numerical value, the lower line gives the corresponding error.

Figure 9: Sobol’ sensitivity analysis results showing the influence of various
physical parameters on the cut-off energy of H-ions for a TNSA water leaf
target experiment, for the full model utilizing only the H2O data. Errors are
given in the $95\text{\,}\mathrm{\char 37\relax}$ confidence level. Figure 10:
PAWN indices for the full model as a measure for parameter importance. Boxes
consist of uncertainty value, minimum, median, maximum, and upper uncertainty
value. The numerical value given is the median.

As can be seen in the display of the data, the mixture has, according to the
Sobol’ analysis, minimal if not zero influence on the maximum energy of the
hydrogen component, while the PAWN analysis gives a higher influence.
Furthermore, the variance of the output can be better explained in this model,
than from the reduced model although the geometry did not change.

#### Sensitivity Analysis Discussion

We performed a sensitivity analysis on our models and were able to evaluate
the importance of the different parameters. We found evidence that the model
describing the laser-plasma system is highly non-linear. It should be noted
that a deep learning model approximates the physical system very well. It does
not, however, provide a closed-form solution for the underlying physics, which
would require further theoretical work.

The models apply regression for the simulated data and as such are able to
reproduce a mean curve for the data, which is for example displayed in Figure
6 or in Figure 4. The spectrum models take the energy bin value into account
and predict the continuum of accelerated ions. This means, that the dependency
on the bin’s energy would be included as well. An explainable analysis of
taking all energy bin values into account then becomes infeasible as each bin
would require a separate Sobol/PAWN analysis.

The cut-off energy of a TNSA spectrum is the main parameter investigated in
the literature, which is for example analyzed by Zimmer et al. [15]. We have
seen that we can get similar results to Zimmer et al. for the cut-off energy
dependencies. We see from the Sobol’ analysis, that several parameters are of
importance, therefore having only a single parameter to describe the cut-off
is not sufficient. Since second-order Sobol’ indices are not zero, we have to
take them into account as well.

Since neither model is explaining the cut-off variation close to
$100\text{\,}\mathrm{\char 37\relax}$, when only 1st and 2nd-order variations
are taken into account, we can conclude that the calculated models require
consideration of even higher order variations to describe an additional
$10\text{\,}\mathrm{\char 37\relax}$-$20\text{\,}\mathrm{\char 37\relax}$ of
the cut-off variation. Such a large reliance on higher-order interactions
implies that simple scaling models are unideal for optimizations since these
effects are not taken into account. A model capable of approximating highly
non-linear effects, such as our neural network models, should thus be
preferred.

The Sobol’ indices decompose the function into a unique space[39], this could
be used to construct a polynomial chaos expansion[44] polynomial from it. This
polynomial can describe the same amount of variation as indicated by first and
second-order Sobol’. It is therefore neither a complete representation of our
network-based models nor is it physically interpretable.

###  III.3 Interpretations

Two major observations from the numerical study are of interest for the
understanding of the modeled system.

The first observation is the deviation from the exponential Mora-like shape
towards the plateau-like features as presented in Figure 1. An explanation of
this effect is the particle-particle interaction inside the expanding plasma.
The driver of this effect is the higher-mass particle species (heavy ions),
which is accelerated later than the lower-mass particle species (protons).
Heavy ions, with their higher inertia co-propagate with lower energy protons
and interact via the Coulomb force. Due to the higher inertia, the protons are
pushed away from the heavy ions, being accelerated as a result. This effect is
especially highlighted by 1D PIC simulations since no transverse particle
movement is allowed. For higher order dimensions [10] or experimental data
[21] the effect is less dominant and the transitions are smoother. If the
particles are accelerated purely in the longitudinal direction, the
divergence, which is given by the quotient of transversal and longitudinal
momentum, can be reduced as well [45, Eq. 2].

The second observation is the high increase in energy absorption. The fraction
of the energy passed onto the protons increases by a factor of about 42. As
shown in the validation for the angular Lorentz boost scheme in Appendix C, a
large increase in the absorption efficiency is a result of the angle-dependent
resonance absorption. Indications for this are displayed in Figure 12 and
Figure 14. The optimization algorithm exploits this behavior directly and
therefore finds ideal angle values. However, there are at least two sides to
this coin. The goal of the approach was to describe the TNSA process in a
model which allows for the optimization of the output depending on the input.
Approaching this directly and analytically is not possible. The time
development of the governing Maxwell-Vlasov system, which already is a
simplification using the collision-free case, can not be solved in closed
form. No true relations for the cut-off energy, for example, have been derived
so far. To get as close to this ground truth as possible, and to become able
to extract it at a later point (with sufficient experimental data), a complex
numerical model must be used. In our case, we adopted an artificial neural
network approach. Artificial neural networks have desirable properties as they
have been shown to be universal function approximators [46]. However, the
explainability of such complex models has been a critical point in their
analysis for some time. The sensitivity analyses shown in subsection III.2
were used as a way to mitigate the complexity and gain some explainability of
the model. The Sobol’ indices method, or global variance-based method,
underlines, that the interaction of the different parameters is of importance.
First and second order can only explain 79% of the models’ variance (Figure 7
b). This means that higher-order dependencies of the input parameters are
necessary to explain a significant part (21%) of the variance. The model can
not explain which higher-order effect i.e. which combination of input
quantities, is exactly responsible. The model’s goal is to allow engineering
optimization of the TNSA process. As a result of this optimization, this
higher-order dependency was found.

##  IV Conclusion

In this study, we modeled and optimized a possible TNSA experiment using a
liquid leaf target by employing a combination of Particle-In-Cell simulations
and deep learning. In agreement with previous studies [10, 21], we have seen
that the accelerated spectra from a multi-species target behave untypical in
comparison to regular one-species TNSA which is described by Mora.

We developed surrogate models that replicate computationally costly PIC
simulations using a deep learning approach. Deep learning is well-suited for
optimizing complex systems. To take advantage of the trained models’ inference
speed, we used the Byrd-Omojokun algorithm to find an optimal parameter
configuration for the system. This yielded a set of parameters that resulted
in optimal maximum hydrogen energy (8 times greater than the initial
parameters) and a set of parameters that resulted in optimal laser energy
conversion efficiency (41 times greater than the initial parameters). We
verified these findings with additional PIC simulations.

We applied sensitivity analysis methods to evaluate the influence of the
different parameters and successfully identified the relevant ones. We showed
that such sensitivity analysis methods bear great potential for the
understanding and quantification of physical dependencies when a closed-form
solution is not known.

The data-based model that we developed can be extended in the future to
improve predictions and better understand the system. This can be achieved by
incorporating future experimental data for the liquid jet.

## Author Declarations

The authors have no conflicts of interest to disclose.

## Data availability statement

Codes and data are available on request.

## Funding Statement

This work was funded by HMWK through the LOEWE center “Nuclear Photonics.”

This work is also supported by the Graduate School CE within the Centre for
Computational Engineering at Technische Universität Darmstadt.

The results presented here are based on simulations, which were performed on
the Virgo HPC cluster at the GSI Helmholtzzentrum für Schwerionenforschung,
Darmstadt (Germany) in the frame of FAIR Phase-0.

###### Acknowledgements.

We would like to thank the Smilei team for providing valuable discussions. We
would also like to thank Ion Gabriel Ion and Dimitrios Loukrezis for the
helpful discussions on sensitivity analysis and the physical interpretation of
artificial neural network models.

##  Appendix A Units and Dimensionality

The simulations in this work were done using the particle-in-cell (PIC) method
[47]. In the following section, we discuss the units and dimensions of the
underlying Maxwell Vlasov system and extract a lower number of relevant
parameters which give valuable physical insight. This is important to
understand why 9 parameters were used in our model. The basis we construct in
this chapter can be represented by the physical parameters we sampled for the
simulation part of our study.

###  A.1 Basis Maxwell Vlasov System and Normalization

TNSA requires a high-intensity laser pulse to heat plasma electrons up to MeV
temperatures. We assume that the mean free path is larger than the target
thickness and the whole process can therefore be assumed as collision-free
[48, 49]. If the process is collision-free, it can be characterized by the
Maxwell-Vlasov system of partial differential equations.

$\displaystyle\nabla\cdot\vec{B}$ $\displaystyle=0$
$\displaystyle\nabla\cdot\vec{E}$
$\displaystyle=\frac{\varrho}{\varepsilon_{0}}$
$\displaystyle\nabla\times\vec{E}$
$\displaystyle=-\frac{\partial\vec{B}}{\partial t}$ (7)
$\displaystyle\nabla\times\vec{B}$
$\displaystyle=\mu_{0}\vec{j}+\mu_{0}\varepsilon_{0}\frac{\partial\vec{E}}{\partial
t}$ $\displaystyle 0$ $\displaystyle=\frac{\partial f_{\alpha}}{\partial
t}+\vec{v}_{\alpha}\cdot\nabla
f_{\alpha}+q_{\alpha}\left(\vec{E}+\vec{v}_{\alpha}\times\vec{B}\right)\cdot\frac{\partial
f_{\alpha}}{\partial\vec{p}}$

Solving these coupled equations efficiently with numerical methods makes it
important to simplify relations. A normalization towards reference quantities
is the first step:

$\displaystyle t^{\prime}$ $\displaystyle=\frac{t}{\tau}\quad\text{where
}\tau\text{ is the pulse length}$ (8) $\displaystyle r^{\prime}$
$\displaystyle=\frac{r}{L}\quad\text{where $L$ is the focus size on the
target}$ (9) $\displaystyle p^{\prime}$
$\displaystyle=\frac{p}{p_{0}}\quad\text{with}\quad
p_{0}=\frac{eE_{0}}{\omega_{\text{L}}}$ (10) $\displaystyle\vec{E}$
$\displaystyle=E_{0}\hat{\vec{E}}$ (11) $\displaystyle\vec{B}$
$\displaystyle=\frac{E_{0}}{c}\hat{\vec{B}}$ (12)

The charge distribution and the current can further be expressed by

$\displaystyle\varrho$
$\displaystyle=\int\sum_{\alpha}f_{\alpha}\operatorname{d}^{3}\\!\vec{p}\quad\text{and}$
(13) $\displaystyle\vec{j}$
$\displaystyle=\int\sum_{\alpha}f_{\alpha}\vec{v}_{\alpha}\operatorname{d}^{3}\\!\vec{p}\quad\text{with}\quad\vec{v}_{\alpha}=\frac{\vec{p}}{m_{\alpha}}\left(1+\frac{p^{2}}{m_{\alpha}^{2}c^{2}}\right)^{-1/2},$
(14)

where $f_{\alpha}$ denotes the charge density. These can be normalized to a
reference quantity as well by modifying $f_{\alpha}$ accordingly shifting all
dimensions into the new $n_{\alpha}$

$\hat{f}_{\alpha}=\frac{f_{\alpha}}{n_{\alpha}}$ (15)

###  A.2 Similitude Relations and Dimensional Reduction

A system of equation can be simplified, by applying the Buckingham $\Pi$
theorem[50]. This theorem allows us to take the dimensional quantities of a
problem into account and find underlying dimensionless quantities which
reflect the actual physical meaning.

If the boundary and initial conditions are similar, then fewer dimensionless
parameters than dimensional parameters can be found to fully represent this
equation system. This implies that the shape of the electromagnetic wave,
defined by $\vec{B}$ and $\vec{E}$, and the normalized charge density
$\hat{f}_{\alpha}$, have to be similar. Similar in this case means, that the
governing function is the same except for some parameters which themselves can
be derived using the Buckingham $\Pi$ Theorem as well.

To ensure similarity in this work, a Gaussian profile was assumed for the
electromagnetic wave, leaving the laser frequency $\omega_{\text{L}}$, the
pulse length $\tau_{\text{L}}$, and the corresponding electric peak field
$E_{0}$ as variable quantities. The initial plasma distribution is defined as
a homogeneous slab with particle density $n_{0}$ and a thickness
$d_{\text{T}}$ with exponential decaying pre-plasma and skirt. The exact
relations are given in subsection II.2.

Keeping these initial conditions fixed allows us to apply the Buckingham $\Pi$
Theorem to the Maxwell-Vlasov system of equations. This results in
dimensionless quantities $\Pi_{i}$ which are capable of describing all
dimensional quantities inside the equation system. The dimensional quantities
are given in Table 4.

Table 4: Overview of the dimensional quantities of the Maxwell-Vlasov EQS. Dimensions are listed in SI base dimensions. Buckingham $\Pi$ parameters are calculated by defining primary quantities which are used multiplicatively in each parameter. Quantity | Dimensions | Type
---|---|---
$\tau$ | T1 |
$L$ | L1 |
${q}/{m}$ | C1 T1 Mass-1 | Primary
$E_{0}$ | M1 L1 C-1 T-3 |
$\omega$ | T-1 | Primary
$\mu_{0}$ | M1 L1 T-2 C-2 | Primary
$\varepsilon_{0}$ | M-1 L-3 T4 C2 | Primary
$n_{\alpha}$ | C1 T1 L-3 |
T: Time, L: Length, C: Current, M: Mass

Using these dimensional quantities, the following $\Pi_{i}$ are determined:

$\displaystyle\Pi_{1}$ $\displaystyle=\omega_{\text{L}}\tau_{\text{L}}$ (16)
$\displaystyle\Pi_{2}$
$\displaystyle=\omega_{\text{L}}L\sqrt{\mu_{0}\varepsilon_{0}}=\frac{\omega_{\text{L}}L}{c}$
(17) $\displaystyle\Pi_{3}$
$\displaystyle=\frac{q_{\alpha}n_{\alpha}}{\varepsilon_{0}m_{\alpha}\omega_{\text{L}}^{2}}=\begin{cases}\dfrac{\hat{n}_{\text{e}}e^{2}}{\varepsilon_{0}m_{\text{e}}\omega_{\text{L}}^{2}}\quad\text{for
electrons}\\\\[17.22217pt]
\dfrac{\hat{n}_{\alpha}Z_{\alpha}^{2}e^{2}}{\varepsilon_{0}m_{\alpha}\omega_{\text{L}}^{2}}\quad\text{for
ions}\end{cases}$ (18) $\displaystyle\Pi_{4}$
$\displaystyle=\frac{q_{\alpha}E_{0}}{m_{\alpha}\omega_{\text{L}}}\sqrt{\mu_{0}\varepsilon_{0}}=\begin{cases}\dfrac{E_{0}}{\omega_{\text{L}}c}\dfrac{e}{m_{\text{e}}}\quad\text{for
electrons}\\\\[17.22217pt]
\dfrac{E_{0}}{\omega_{\text{L}}c}\dfrac{Z_{i}e}{m_{i}}\quad\text{for
ions}\end{cases}$ (19)

The theorem also states that the number of resulting dimensionless parameters
is lower than the number of dimensional parameters, reducing the complexity of
the model.

###  A.3 Interpretation of the dimensionless quantities

These quantities are sufficient to describe and condition the EQS from a
mathematical standpoint. From a physical standpoint, this also creates
valuable insight. $\Pi_{1}$ gives the number of $\vec{E}$-oscillations in the
laser pulse and $\Pi_{2}$ the irradiation size of the laser. $\Pi_{3}$
correlates the laser with the target since it is the ratio of the particle
density to the critical plasma density defined by the laser. $\Pi_{4}$
describes the particle dynamic inside the laser’s amplitude for each species.
For electrons, $\Pi_{4}$ is identical to the dimensionless quiver velocity
$a_{0}$. The meaning is equivalent for the different ion species.

Writing down the EQS and substituting the $\Pi$ parameters makes their
importance apparent:

$\displaystyle\tilde{\nabla}\cdot\hat{\vec{B}}$ $\displaystyle=0$ (20)
$\displaystyle\tilde{\nabla}\cdot\hat{\vec{E}}$
$\displaystyle=\underbrace{\int\sum_{\alpha}\frac{\Pi_{2}\cdot\Pi_{3\alpha}}{\Pi_{4\alpha}}\operatorname{d}^{3}\tilde{\vec{p}}}_{=0\text{
for }t=0}$ (21) $\displaystyle\tilde{\nabla}\times\hat{\vec{E}}$
$\displaystyle=-\frac{\Pi_{2}}{\Pi_{1}}\frac{\partial\hat{\vec{B}}}{\partial\tilde{t}}$
(22) $\displaystyle\tilde{\nabla}\times\hat{\vec{B}}$
$\displaystyle=\frac{\Pi_{2}}{\Pi_{1}}\frac{\partial\hat{\vec{E}}}{\partial\tilde{t}}+\int\sum_{\alpha}\frac{4\pi\Pi_{2}\Pi_{3\alpha}\tilde{\vec{p}}}{\sqrt{\tilde{p}^{2}\Pi_{4\alpha}^{2}+1}}\hat{f}_{\alpha}\operatorname{d}^{3}\tilde{\vec{p}}$
(23) $\displaystyle 0$
$\displaystyle=\frac{1}{\Pi_{1}}\frac{\partial\hat{f}_{\alpha}}{\partial\tilde{t}}+\frac{1}{\Pi_{2}}\frac{\Pi_{4\alpha}\tilde{\vec{p}}}{\sqrt{\tilde{p}^{2}\Pi_{4\alpha}^{2}+1}}\frac{\partial\hat{f}_{\alpha}}{\partial\tilde{\vec{r}}}$
$\displaystyle+\left({Z_{\alpha}}\hat{\vec{E}}+\frac{\Pi_{4\alpha}}{\sqrt{\tilde{p}^{2}\Pi_{4\alpha}^{2}+1}}\tilde{\vec{p}}\times\hat{\vec{B}}\right)\frac{\partial\hat{f}_{\alpha}}{\partial\tilde{p}}$
(24)

If these parameters are constant, then the equations are all the same and
therefore behave the same. This results in the same time development of the
system and therefore yields the same results. One can say that for constant
$\Pi_{i}$ areas of iso-dynamics exist, finally simplifying any model
approaches by reducing the dimensions to be examined. Models therefore only
need these 4 parameters to precisely determine a system.

###  A.4 Correlating dimensionless parameters and simulation input

The system, therefore, has a dedicated number of $\Pi$ Parameters which have
to be taken into account: $\Pi_{1}$ and $\Pi_{2}$ are laser-relevant
quantities and therefore particle species independent. $\Pi_{3}$ and $\Pi_{4}$
describe quantities of the particle species, therefore introducing a
multiplicity in the parameter, denoted by $\alpha$. In the case investigated
here, the multiplicity is 4: electrons, oxygen, hydrogen, and deuterium. The
system is initialized with the same spatial distribution function
$\hat{f}_{\alpha}$ for each species. To mimic ionization and mixture some
conditions apply:

$N_{O}=2\times N_{H/D}\quad\text{and}\quad N_{e}=Z^{\text{eff}}_{O}+Z_{H/D}$
(25)

Taking these assumptions into account resolves the multiplicity and the
corresponding $\Pi$s can be expressed with a multiplicative factor. The
construction, including the multiplicative factors and the needed parameters,
are given in Table 5. Red and green mark the relevant, varying parameters to
be taken into account.

Table 5: Construction of $\Pi$ parameters. Parameter | Definition
---|---
$\Pi_{1}$ | ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{L}}\tau_{\text{L}}$
$\Pi_{2}$ | ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{L}}L/c$
$\Pi_{3O}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}Z}^{\text{eff}}_{O}/m_{O}\times n_{O}/{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{L}}^{2}$
$\Pi_{3e}$ | $\Pi_{3O}\times\left({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}Z}^{\text{eff}}_{O}+2\right)\times\frac{m_{O}}{m_{e}}\cdot\frac{q_{e}}{e\cdot{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}Z}^{\text{eff}}_{O}}$
$\Pi_{3H}$ | $\Pi_{3O}\times 2\text{{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}mix}}\times\frac{m_{O}}{m_{H}}\cdot\frac{q_{H}}{e\cdot{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}Z}^{\text{eff}}_{O}}$
$\Pi_{3D}$ | $\Pi_{3O}\times 2\left(1-\text{{\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}mix}}\right)\times\frac{m_{O}}{m_{D}}\cdot\frac{q_{D}}{e\cdot{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}Z}^{\text{eff}}_{O}}$
$\Pi_{4e}$ |  $q_{e}/m_{e}\times E_{0}/{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\omega_{L}}c$ ($a_{0}$)
$\Pi_{4H}$ | $\Pi_{4e}\times\frac{m_{e}}{m_{H}}\cdot\frac{q_{H}}{q_{e}}$
$\Pi_{4D}$ | $\Pi_{4e}\times\frac{m_{e}}{m_{D}}\cdot\frac{q_{D}}{q_{e}}$
$\Pi_{4O}$ | $\Pi_{4e}\times\frac{m_{e}}{m_{O}}\cdot\frac{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}Z}^{\text{eff}}_{O}}{q_{e}}$

#### Mapping physical to dimensionless parameters

As stated in the main body of this work, several physical input quantities are
used. They are chosen based on keeping datasets consistent and comparable.
Therefore some parameters are sampled which do not exist in 1D. This also
ensures that the data can be broken down into the $\Pi_{i}$ with the following
relations. It is important to note that from the $\Pi_{4}$ possibilities only
the electron variant (equivalent to $a_{0}$) has to be passed to the PIC code.

$\displaystyle\Pi_{1}$ $\displaystyle=\frac{2\pi
c\tau_{\text{L}}}{\lambda_{\text{L}}}$ (26) $\displaystyle\Pi_{2}$
$\displaystyle=\frac{2\pi\text{FWHM}}{\lambda_{\text{L}}}$ (27)
$\displaystyle\Pi_{3}$
$\displaystyle=\frac{Z_{\text{eff}}^{2}e^{2}n_{O}\lambda_{\text{L}}^{2}}{\varepsilon_{0}4\pi^{2}m_{O}c^{2}}$
(28) $\displaystyle\Pi_{4}$
$\displaystyle=\left(\frac{E_{\text{L}}\lambda_{\text{L}}^{2}\sin\left(\Theta_{\text{L}}\right)}{\tau_{\text{L}}\pi\text{FWHM}^{2}\cdot
1.37\times 10^{18}}\right)^{2}$ (29)

This culminates in a needed dimensionality of 9 for the list of parameters: 4
$\Pi_{i}$, 2 parameters to deal with the ambiguity $Z^{\text{eff}}_{O}$ and
the mixture parameters, 1 parameter for the plasma slab $d_{T}$ (particle
density is fixed) and 2 for dealing with laser’s polarization: Selection
whether p / s linear polarization and to make a difference, variation of the
incidence angle $\Theta_{L}$.

Taking the mapping into account (Eqs. (26)–(29)). A proper physical sampling
includes:

  1. 1.

Ionization of oxygen

  2. 2.

Mixtures (deuterium vs hydrogen)

  3. 3.

Laser Polarization

  4. 4.

Laser Energy / Joule

  5. 5.

Laser Pulse Time / second

  6. 6.

Laser Irradiation Size / micron

  7. 7.

Laser Wavelength / meter

  8. 8.

Laser Incidence Angle / Degree (to the plasma normal)

  9. 9.

Plasma Slab Thickness

These 9 parameters have the same dimensionality as the parameter space
calculated by the $\Pi_{i}$ which is necessary since the construction of
dedicated quantities (especially the electric field of the laser) cannot be
determined easily and a composition of these parameters has to be taken into
account.

##  Appendix B Parameter Ranges

As mentioned in subsubsection II.2.2 two paradigms are relevant for the
selection of the parameter ranges.

The class of Petawatt laser systems we based our work on uses mainly linearly
polarized laser light and is capable of varying the incidence angle. Taking
this into account we get two possibilities for the laser’s polarization:
s-polarization and p-polarization. Similarly, we can get angles from
$0\text{\,}\mathrm{\SIUnitSymbolDegree}$ to less than
$90\text{\,}\mathrm{\SIUnitSymbolDegree}$. At
$90\text{\,}\mathrm{\SIUnitSymbolDegree}$ the laser is not hitting the target
and is traveling parallel to the plasma surface, we, therefore, chose to cut
the interval at $85\text{\,}\mathrm{\SIUnitSymbolDegree}$.

The laser energy was selected to cover a large area to increase the
comparability of the model with different laser systems. At the
conceptualization phase of this study, it was unreasonable to assume high
repetition rate experimentation with much larger systems (e.g. GSI’s PHELIX
system[51, 52]), since the currently achievable repetition rates were too low.
This might change in the future, such that higher energies are realistic, and
the model base has to be expanded under such cases.

The most critical parameter is the pulse length. As our reference value, we
selected the FWHM in the time domain of a pure Gaussian pulse. Firstly, the
approximation of a pure Gaussian pulse is not necessarily true for a
technically implemented laser. If the FWHM according to Equation 4 is not
used, then the value has to be adjusted accordingly. Due to calculation time
issues of the underlying PIC models, we reduced the selected times to the
interval from $15\text{\,}\mathrm{fs}$ to $150\text{\,}\mathrm{fs}$. We know
that this time can be significantly larger, but high repetition systems can
operate with low pulse length variables. We also acknowledge, that our lower
simulation border for the time is close to the bandwidth-limited pulse limit,
but we wanted to have some lower data points to force the interpolation into
good behavior and therefore mathematically overshot into the lower regime. The
upper pulse length boundary is also the first parameter we want to increase in
further studies since technical laser systems do need a larger pulse length to
apply this model.

The focus FWHM was then sampled according to the definition of the laser
$a_{0}$ conditions which we applied to our parameters. We made sure to stay in
the TNSA regime and prevented the laser $a_{0}$ to be smaller than one and
keep the focus still realistically small with
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ as the lower range cut-off.

The selected wavelengths are larger than those used in engineered systems and
also are somehow continuously sampled from this larger range. The reason for
this is the importance of the wavelength parameter following the similitude
relations, which are dependent on the laser wavelength in every component.

Concerning the target, we chose the thickness according to the parameters of
the physical implementation of a liquid jet, which is currently under
development. The mixture can only vary between 0 and
$100\text{\,}\mathrm{\char 37\relax}$. Again, the effective charge of the
particles plays an important role. We started with fully ionized oxygen and
found traits of the multi-species effect during our investigation. While
discussing the results, we generalized the effective charge discussion but
were not able to properly simulate different effective ionization levels. This
is due to a lack of proper ionization models (also beyond the scope of this
study) and limited numerical resources. This would also be a parameter that
could be further improved in additional studies.

##  Appendix C Transversal Lorentz Boosted 1.5D PIC Simulations

Modeling oblique laser incidence onto a target is inherently at least a 2D
problem, which requires substantially more computational power than a similar
1D geometry to simulate. Bourdier [53] thus proposed a method in which a
relativistic Lorentz boost is applied to the frame of reference in the
simulation. This method has later been employed by Gibbon et al. [54] in a PIC
code.

Here, we would like to present the implementation of this technique yet again
for a modern PIC code while also correcting some mistakes in the calculations
by Gibbon et al.. A schematic of the general principle is shown in Figure 11.
To obtain the results in the lab frame a back transformation must be applied
to the diagnostics obtained from the simulation.

Figure 11:  Schematic of the Lorentz boosted simulation frame versus the
implied lab frame. In the simulation frame, the laser appears to be at normal
incidence onto the target while the particles appear to drift in negative
$y$-direction with velocity $v_{y}$.

For finding the transformations, a simple Lorentz boost in $y$-direction by
the velocity $v_{y}=c\cdot\sin(\theta)$ is applied. In matrix form, this can
be represented by:

$\displaystyle\Lambda$
$\displaystyle=\begin{pmatrix}1/\cos(\theta)&0&-\tan(\theta)&0\\\ 0&1&0&0\\\
-\tan(\theta)&0&1/\cos(\theta)&0\\\ 0&0&0&1\end{pmatrix}\;,$ (30)

since $\gamma_{0}=\left(1-(v_{y}/c)^{2}\right)^{-1/2}=1/\cos(\theta)$. This
transformation matrix can be used to transform all the quantities of the
particles and the electromagnetic fields. Indicating quantities in the
transformed system with a prime, we find after carrying out all
transformations

$\displaystyle\begin{split}k^{\prime}_{y}&=0\\\
\omega_{\text{L}}^{\prime}&=\omega_{\text{L}}/\gamma_{0}\\\
a_{0}^{\prime}&=a_{0}\;,\end{split}$ (31)

where $k^{\prime}_{y}$ is the y-component of the wave vector in the boosted
system, showing that, indeed, the laser is now at normal incidence. Note that
the dimensionless laser amplitude $a_{0}$ is invariant under the
transformation [54]. Further, denoting PIC code units with a tilde, we find

$\displaystyle\begin{split}\tilde{x}^{\prime}&=\tilde{x}/\gamma_{0}\\\
\tilde{t}^{\prime}&=\tilde{t}/\gamma_{0}^{2}\end{split}$ (32)

giving a re-scaling of both the simulation time as well as the cell grid. The
initial particle density is also affected by

$\displaystyle\tilde{n}_{0}^{\prime}$
$\displaystyle=\tilde{n}_{0}\cdot\gamma_{0}^{3}\;.$ (33)

With these conditions, the particles can be initialized in the boosted frame.
The relative velocity $v_{y}$ is added as a permanent drift which is handled
and relativistically added to the particles by the code.
From the diagnostics in the simulation, we can obtain desired quantities via a
back transformation. For the particle kinetic energies, we find using the
energy-momentum relation:

$\displaystyle\begin{split}\tilde{E}&=\gamma_{0}(\tilde{E}^{\prime}+\tilde{v}_{y}\tilde{p}_{y}^{\prime})\\\
&=\gamma_{0}\left(\frac{1}{m_{\text{e}}c^{2}}\sqrt{{p^{\prime}}^{2}c^{2}+(m_{0}c^{2})^{2}}+\tilde{v}_{y}\tilde{p}_{y}^{\prime}\right)\\\
\Rightarrow
E_{\text{kin}}&=\gamma_{0}m_{\text{e}}c^{2}\left(\sqrt{\tilde{p}^{\prime
2}+\frac{m_{0}^{2}}{m_{\text{e}}^{2}}}+\tilde{v}_{y}\tilde{p}_{y}^{\prime}\right)-m_{0}c^{2}\;,\end{split}$
(34)

where $m_{0}$ is the particle rest mass. Noting that
$\tilde{B}_{x}^{\prime}=0$, we can find relations to recover the fields of the
laser. The non-zero fields are:

For s-polarization:
$\displaystyle\begin{split}&\tilde{E}_{z}=\tilde{E}_{z}^{\prime}\\\
&\tilde{B}_{x}=\tilde{v}_{y}\tilde{E}_{z}^{\prime}\\\
&\tilde{B}_{y}=\tilde{B}_{y}^{\prime}/\gamma_{0}\end{split}$ (35) For
p-polarization:
$\displaystyle\begin{split}&\tilde{E}_{x}=\tilde{E}_{x}^{\prime}-\tilde{v}_{y}\tilde{B}_{z}^{\prime}\\\
&\tilde{E}_{y}=\tilde{E}_{y}^{\prime}/\gamma_{0}\\\
&\tilde{B}_{z}=\tilde{B}_{z}^{\prime}-\tilde{v}_{y}\tilde{E}_{x}^{\prime}\end{split}$
(36)

Using the field transformations and assuming that reflection at the plasma
surface does not change polarization, we find for the absolute magnitude of
the Poynting vector:

$\displaystyle|\vec{S}|$
$\displaystyle=|\vec{S}^{\prime}|\cdot\gamma_{0}^{2}\;,$ (37)

with which the relative absorption of the laser into the plasma can be
calculated by dividing the incoming Poynting flux by the outgoing Poynting
flux.

It should be noted that while the Lorentz boosted frame method can replicate
incidence angle-based behavior, it cannot replace a 2D or even 3D simulation
on all accounts [54]. Firstly, in the general case, all physical quantities
depend separately on the transformed coordinates $x,y,z,t,p_{x},p_{y},p_{z}$.
Thus, the Lorentz-boosted simulation can only be used for a problem
independent of $y$ and $z$. Additionally, reducing the geometry after the
boost to 1D limits the spatial dynamics of the particles. Since only the
$x$-axis is present, all particles (while having 3D velocities) can only move
along a straight line (i.e. have only 1 spatial dimension). This disregards
the angular spread at the back of the target such that the particles can be
accelerated for longer times and thus end up with higher energies compared to
a similar 2D simulation. Distinctly 2D effects such as hole boring can also
not be modeled accurately.

To illustrate the capabilities of this method, however, the relative laser
absorption of a p-polarized laser impinging on a hydrogen plasma target was
measured for varying laser incidence angles using the above method in the
Smilei PIC code. The resulting absorption curve is shown in Figure 12. The
results agree well with 2D simulations by Cui et al. [55] using a similar
target and laser (see Figure 14 for a comparison).

Figure 12:  Simulation of laser incidence angles between
$0\text{\,}\mathrm{\SIUnitSymbolDegree}$ and
$85\text{\,}\mathrm{\SIUnitSymbolDegree}$. The plot shows the incidence angle
versus the relative absorption of a laser into a hydrogen plasma target
(dotted red line) as well as the maximum proton kinetic energy behind the
target (blue solid line). The laser impinges on the target with p-polarized
fields. Classical resonance absorption, also known as the Desinov curve [56],
is shown as a dashed line.

### Explicit Lorentz Boost for oblique Laser Incidence

In the following section we discuss the full transformation in more detail,
and explicitly calculate the relations we mentioned before. Starting from the
transformation matrix in Equation 30 the full derivation will be done for all
quantities in the system.

Firstly, the four-position $R$, the four-momentum $P$, the four-wave vector
$K$ and the four-current $J$ are given as follows:

$\displaystyle\begin{split}R&=\left(ct,x,y,z\right)^{\top}\\\ P&=\left(\gamma
m_{0}c,p_{x},p_{y},p_{z}\right)^{\top}\\\
K&=\left(\omega_{\text{L}}/c,k\cdot\cos(\theta),k\cdot\sin(\theta),0\right)^{\top}\\\
J&=\left(c\rho,j_{x},j_{y},j_{z}\right)^{\top}\end{split}$ (38)

where $k=\omega_{\text{L}}/c$ is the magnitude of the wave vector, $\rho$ is
the charge density and $\vec{j}$ is the current density. Here, the geometry of
the wave vector from Figure 11 has already been applied, reducing the wave
vector to two spatial dimensions. By left multiplication of $\Lambda$ these
quantities can be transformed into the boosted frame. This multiplication
yields

$\displaystyle\begin{split}R^{\prime}&=\left(\gamma_{0}(ct-y\beta_{0}),x,\gamma_{0}(y-ct\beta_{0}),z\right)^{\top}\\\
P^{\prime}&=\left(\gamma_{0}(\gamma
m_{0}c-p_{y}\beta_{0}),p_{x},\gamma_{0}(p_{y}-\gamma
m_{0}c\beta_{0}),p_{z}\right)^{\top}\\\
K^{\prime}&=\left(\omega_{\text{L}}/(c\gamma_{0}),k_{0}/\gamma_{0},0,0\right)^{\top}\\\
J^{\prime}&=\left(\gamma_{0}(c\rho-
j_{y}\beta_{0}),j_{x},\gamma_{0}(j_{y}-c\rho\beta_{0}),j_{z}\right)^{\top}\end{split}$
(39)

where a prime indicates quantities in the transformed system and
$\beta_{0}=v_{y}/c=\sin(\theta)$. Most importantly here we find
$k^{\prime}_{y}=0$ and
$\omega^{\prime}_{\text{L}}=\omega_{\text{L}}/\gamma_{0}$. Also, since the
particles are assumed cold at $t=0$, we find for the initial density
$\rho^{\prime}_{0}=\gamma_{0}\rho_{0}$.
The next transformation is for the electromagnetic fields. Here, we
differentiate between s- and p-polarized incidence lasers. To transform the
electric and magnetic fields of the incoming laser, the electromagnetic tensor
is used:

$\displaystyle F_{\text{s-pol}}^{\mu\nu}$
$\displaystyle=\begin{pmatrix}0&0&0&-E_{z}/c\\\ 0&0&0&B_{y}\\\ 0&0&0&-B_{x}\\\
E_{z}/c&-B_{y}&B_{x}&0\end{pmatrix}$ (40) $\displaystyle
F_{\text{p-pol}}^{\mu\nu}$
$\displaystyle=\begin{pmatrix}0&-E_{x}/c&-E_{y}/c&0\\\ E_{x}/c&0&-B_{z}&0\\\
E_{y}/c&B_{z}&0&0\\\ 0&0&0&0\end{pmatrix}\;.$ (41)

The Lorentz transformation of such a tensor is given by:

$\displaystyle F^{\mu^{\prime}\nu^{\prime}}$
$\displaystyle=\Lambda^{\mu^{\prime}}{}_{\mu}\Lambda^{\nu^{\prime}}{}_{\nu}F^{\mu\nu}\;,$
(42)

where a prime again indicates quantities in the transformed system. The
calculated fields are

$\displaystyle\left.\begin{aligned} E^{\prime}_{x}&=0\\\ E^{\prime}_{y}&=0\\\
E^{\prime}_{z}&=\gamma_{0}(E_{z}-v_{y}B_{x})\\\
B^{\prime}_{x}&=\gamma_{0}(B_{x}-E_{z}v_{y}/c^{2})\\\ B^{\prime}_{y}&=B_{y}\\\
B^{\prime}_{z}&=0\end{aligned}\right\\}\;\text{s-pol}$ (43)
$\displaystyle\left.\begin{aligned}
E^{\prime}_{x}&=\gamma_{0}(E_{x}+v_{y}B_{z})\\\ E^{\prime}_{y}&=E_{y}\\\
E^{\prime}_{z}&=0\\\ B^{\prime}_{x}&=0\\\ B^{\prime}_{y}&=0\\\
B^{\prime}_{z}&=\gamma_{0}(B_{z}+E_{x}v_{y}/c^{2})\end{aligned}\right\\}\;\text{p-pol}$
(44)

where $E^{\prime}_{x},B^{\prime}_{x}\overset{!}{=}0$ since the laser is at
normal incidence in the boosted system. For absorption measurements, it is
useful to have a look at the transformation of the Poynting Vector $\vec{S}$.
We first define in vacuum

$\displaystyle\vec{S}$
$\displaystyle=\frac{1}{\mu_{0}}\vec{E}\times\vec{B}\;,$ (45)

where $\mu_{0}$ is the vacuum permeability. As an example, we will only
present the calculation in the p-polarization case. The s-polarization
calculation is equivalent. We find

$\displaystyle\vec{S}_{\text{p-pol}}$
$\displaystyle=\frac{1}{\mu_{0}}\left(E_{y}B_{z},-E_{x}B_{z},0\right)^{\top}$
(46) $\displaystyle\Rightarrow\vec{S}^{\prime}_{\text{p-pol}}$
$\displaystyle=\frac{1}{\mu_{0}}(E^{\prime}_{y}B^{\prime}_{z},\underbrace{-E^{\prime}_{x}B^{\prime}_{z}}_{=0},0)^{\top}$
(47)
$\displaystyle=\frac{1}{\mu_{0}}\left(E^{\prime}_{y}B^{\prime}_{z},0,0\right)^{\top}\;.$
(48)

We hence find for the magnitude of the transformed Poynting Vector

$\displaystyle|\vec{S}^{\prime}_{\text{p-pol}}|$
$\displaystyle=\frac{1}{\mu_{0}}\sqrt{E_{y}^{\prime 2}B_{z}^{\prime 2}}$ (49)
$\displaystyle=\frac{1}{\mu_{0}c}E_{y}^{\prime 2}\;,$ (50)

since $|\vec{B}|=|\vec{E}|/c$. On the other hand, inserting the transformation
into $\vec{S}$, we find

$\displaystyle\vec{S}_{\text{p-pol}}$
$\displaystyle=\frac{1}{\mu_{0}}\begin{pmatrix}E_{y}^{\prime}\gamma_{0}\left(B^{\prime}_{z}-E^{\prime}_{x}\frac{v_{y}}{c^{2}}\right)\\\
-\gamma_{0}\left(E^{\prime}_{x}-v_{y}B^{\prime}_{z}\right)\gamma_{0}\left(B^{\prime}_{z}-E^{\prime}_{x}\frac{v_{y}}{c^{2}}\right)\\\
0\end{pmatrix}$ (51)
$\displaystyle=\frac{1}{\mu_{0}}\begin{pmatrix}\gamma_{0}E^{\prime}_{y}B^{\prime}_{z}\\\
\gamma_{0}^{2}v_{y}B_{z}^{\prime 2}\\\ 0\end{pmatrix}$ (52)

such that for the magnitude we have

$\displaystyle|\vec{S}_{\text{p-pol}}|$
$\displaystyle=\frac{\gamma_{0}}{\mu_{0}}\sqrt{E_{y}^{\prime 2}B_{z}^{\prime
2}+\gamma_{0}^{2}v_{y}^{2}B_{z}^{\prime 4}}$ (53)
$\displaystyle=\frac{E_{y}^{\prime
2}\gamma_{0}}{\mu_{0}c}\sqrt{1+\gamma_{0}^{2}\frac{v_{y}^{2}}{c^{2}}}\;.$ (54)

The term in the square root can be resolved elegantly once we remind ourselves
of the definition of $v_{y}$:

$\displaystyle 1+\gamma_{0}^{2}\frac{v_{y}^{2}}{c^{2}}$
$\displaystyle=1+\frac{1}{\cos^{2}(\theta)}\frac{c^{2}\sin^{2}(\theta)}{c^{2}}$
(55) $\displaystyle=\frac{1}{\cos^{2}(\theta)}$ (56)
$\displaystyle=\gamma_{0}^{2}\;,$ (57)

and with that we have

$\displaystyle|\vec{S}_{\text{p-pol}}|$ $\displaystyle=\frac{E_{y}^{\prime
2}\gamma_{0}}{\mu_{0}c}\;\gamma_{0}$ (58)
$\displaystyle=|\vec{S}^{\prime}_{\text{p-pol}}|\cdot\gamma_{0}^{2}\;.$ (59)

Next, let us consider the transformed quantities in code units, so as to
initialize the particles correctly in the PIC code. For the space coordinate,
we find

$\displaystyle\tilde{x}^{\prime}$
$\displaystyle=\frac{\omega^{\prime}x^{\prime}}{c}=\tilde{x}/\gamma_{0}\;,$
(60)

while for the time coordinate, since $\tilde{y}^{\prime}\overset{!}{=}0$:

$\displaystyle 0\overset{!}{=}\tilde{y}^{\prime}$
$\displaystyle=\frac{\omega^{\prime}}{c}\gamma_{0}(y-ct\beta_{0})$ (61)
$\displaystyle=\tilde{y}-\tilde{t}\tilde{v}_{y}$ (62)
$\displaystyle\Rightarrow\tilde{t}^{\prime}$
$\displaystyle=\omega^{\prime}\gamma_{0}\left(t-\frac{y}{c}\beta_{0}\right)$
(63) $\displaystyle=\tilde{t}-\tilde{y}\tilde{v}_{y}$ (64)
$\displaystyle=\tilde{t}-\tilde{t}\tilde{v}_{y}^{2}$ (65)
$\displaystyle=\tilde{t}/\gamma_{0}^{2}$ (66)

Finally, the critical density transforms as

$\displaystyle\frac{n^{\prime}_{\text{c}}}{n_{\text{c}}}$
$\displaystyle=\frac{\omega^{\prime
2}}{\omega^{2}}=\frac{1}{\gamma_{0}^{2}}\;,$ (67)

such that the initial particle densities in code units become

$\displaystyle\tilde{n}^{\prime}_{0}$
$\displaystyle=\frac{n^{\prime}_{0}}{n^{\prime}_{\text{c}}}=\tilde{n}_{0}\cdot\gamma_{0}^{3}\;.$
(68)

A verification plot for the Lorentz Boost method is displayed in Figure 13 for
irradiation under an oblique angle.

Figure 13: Example result for a Lorentz-boosted simulation with an angle of 40
degrees. The dashed line denotes the fit of Mora’s model [20] to the data.
Figure 14: Plot of the absorption of energy in the Lorentz boosted simulation
in comparison to the data from Cui et al. [55].

##  Appendix D Laser Conversion Efficiency

The laser conversion efficiency is an important quantity to characterize
particle acceleration and especially laser-plasma acceleration. In order to
retrieve information about the energy in the output spectrum of a TNSA
experiment, consider first a spectrum $\mathrm{d}N/\mathrm{d}E$ recorded in
multiple energy bins of width $\Delta E$. In this case, the number of
particles in bin $i$ is given by the bin’s height multiplied by its width,
i.e.

$\displaystyle N_{i}$
$\displaystyle=\left(\frac{\mathrm{d}N}{\mathrm{d}E}\right)_{i}\cdot\Delta
E\;.$ (69)

Hence, the total energy of the particles within the bin could be approximated
by multiplying the particles in the bin by the bin’s central energy $E_{i}$.
Summing over all bins yields the total energy of the particles

$\displaystyle E_{\text{tot}}$
$\displaystyle=\sum_{i}\left(\frac{\mathrm{d}N}{\mathrm{d}E}\right)_{i}\cdot\Delta
E\cdot E_{i}\;,$ (70)

which can be generalized in the continuous limit $\Delta E\rightarrow 0$,
giving

$\displaystyle E_{\text{tot}}$
$\displaystyle=\int_{0}^{\infty}\frac{\mathrm{d}N}{\mathrm{d}E}\cdot
E\,\mathrm{d}E\;.$ (71)

Concretely, adjusting for the output format of the neural network models the
total energy is given by

$\displaystyle E_{\text{tot}}$
$\displaystyle=V\cdot\int_{0}^{E_{\text{max}}}\exp\left(\ln\left(\frac{\mathrm{d}n}{\mathrm{d}E}\right)\right)\cdot
E\,\mathrm{d}E\;,$ (72)

where $\ln\left(\frac{\mathrm{d}n}{\mathrm{d}E}\right)$ and $E_{\text{max}}$
are given by the neural network models and $V$ is a unit volume. To obtain a
measure for the energy conversion efficiency then, the above integral should
be weighted by the laser pulse energy $E_{\text{L}}$, resulting in the
maximization problem shown in Eq. (6).

##  Appendix E Neural Network Training and preparation

In this section, we discuss the chosen parameter ranges for the surrogate
models based on neural networks.

Training surrogate models is a tedious and numerically expensive task. This
means that we have to be clear about the parameters and data used for the
training process. We will first focus on the data preparation task, and second
on the numerical hyperparameters chosen for our model. Both parts are
important if we want to create fast converging models.

###  E.1 Data preparation

Neural networks can only be as good as the data used for training them.
Convergence is important and data, therefore, has to be prepared properly. We
can only investigate the multi-species effect and subsequent optimizations if
we take the full spectrum into account.

The spectral data for the output spectrum is taken on a logarithmic scale
since the count rates vary over several orders of magnitude. The logarithmic
data can directly be used to train a model. We tried using the data directly,
but convergence was problematic. This is due to the noise of the data and the
mixture-depending shifts of multi-species plateaus. The signal variation in
both cases is similar and it is therefore difficult for the network to fit the
dependencies. To mitigate this we applied a Savitzky-Golay filter [57] with a
window size of 7 points and a 3rd-order polynomial. This filter decreased the
noise-based fluctuations and allowed subsequent convergence. We display a
comparison for filtered and unfiltered data in Figure 15, which shows, that
the major behavior of the curves is reproduced but the bin-to-bin fluctuations
in the mid to high energy range are minimized.

Figure 15: Savitzky-Golay filter with window size 7 and a 3rd order
polynomial. Blue indicates the raw data, projected to a wider range for
visibility and black indicates the filtered signal.

###  E.2 Numerical Parameters, Training and Topology

As none of the architectural parameters for these models were known, some
outlying hyperparameters were decided first. For a regression problem, the
Rectified Linear Unit (ReLU) activation function is widely used and was added
to every layer of the network except the output layer which used the identity
activation. Similarly, we chose the mean squared error, suited for regression
problems, as loss and it was minimized using the Adam optimizer with
$\beta_{1}=0.9$, $\beta_{2}=0.999$ and
$\epsilon=$1\text{\times}{10}^{-07}\text{\,}$$. The initial learning rate was
0.001 which was lowered to a minimum of 0.0001 during training should the
optimizer detect a plateau in the validation loss value (Keras’
ReduceLROnPlateau feature). In order for the physical parameters to be more
manageable numerically, all parameters were divided by the maximum value in
their range (see Table 1) before being given to the model.
With these outlying parameters in place, the architecture of the FCNs, i.e.
the number of layers and the number of neurons in each layer, was left
variable and was optimized for the problem using a hyperparameter tuning
method. Keras Tuner allows for extensive hyperparameter tuning using various
optimization algorithms [31].

Recalling subsection II.2, each simulation output contains information about
100 locations in the energy spectrum of the particles. Hence, for the reduced
continuous model, the available data length was $68973\times 100=6897300$ data
points. Of these, 81% were used for training, 9% were used for validation, and
10% were used for testing. Running Keras Tuner on Google Cloud Compute Engine
API from a Google Colab Notebook, Bayesian Optimization could be performed for
the hyperparameters of the continuous model of hydrogen ions. In order to find
a model architecture that most accurately describes the simulation results the
number of layers and the number of neurons for each layer was first optimized
to achieve the lowest possible training loss. Every training used a batch size
of 256 and an early stopping mechanism. After 50 trials, each running training
twice in order to lower the chance of a bad local minimum, a suitable
architecture was found. However, this optimized model was only tuned to
minimize the training loss of the model without considering the validation
data at all. To generalize the model, hyperparameter tuning was run again on
the optimized architecture, this time with L1 and L2 regularization on each
layer as the hyperparameters to be tuned and with the tuning objective set to
the mean squared error on the validation set. Each hidden layer in the network
has L1 regularization strength of $1.98\text{\times}{10}^{-06}\text{\,}$ and
L2 regularization strength of $3.07\text{\times}{10}^{-08}\text{\,}$. The
network achieved a mean squared error of 3.38 on the 620 757 randomly selected
validation data points. As a reminder, this number is equal to the mean
squared error on the $\ln\left(\frac{\mathrm{d}n}{\mathrm{d}E}(E)\right)$
prediction for input parameters $\\{E,\text{[physical parameters]}\\}$.
Equivalently, the second model predicting the maximum ion energy could be
tuned and optimized. Since the maximum energy is only predicted per simulation
and not per energy bin of the energy spectra, the second model was trained on
68 973 unique data points. This significantly smaller dataset made the model
training on a home computer feasible. The optimized model for the maximum
energy found L1 regularization strength of $2.3\text{\times}{10}^{-4}$ and L2
regularization strength of $1.1\text{\times}{10}^{-7}$.

## References

  * Katayama  _et al._ [2022] R. Katayama, E. Kako, S. Yamaguchi, S. Michizono, K. Umemori, and Y. Kondo, “Design study of compact medical accelerator using superconducting rf quadrupole for boron neutron capture therapy,” Physical Review Accelerators and Beams 25 (2022), 10.1103/physrevaccelbeams.25.021601.
  * Eickhoff  _et al._ [2003] H. Eickhoff, R. Bar, A. Dolinskii, T. Haberer, B. Schlitt, P. Spiller, and U. Weinrich, “HICAT - the german hospital-based light ion cancer therapy project,” in  _Proceedings of the 2003 Bipolar/BiCMOS Circuits and Technology Meeting (IEEE Cat. No.03CH37440)_ (IEEE, 2003).
  * Schardt, Elsässer, and Schulz-Ertner [2010] D. Schardt, T. Elsässer, and D. Schulz-Ertner, “Heavy-ion tumor therapy: Physical and radiobiological benefits,” Reviews of Modern Physics 82, 383–425 (2010).
  * Linz and Alonso [2007] U. Linz and J. Alonso, “What will it take for laser driven proton accelerators to be applied to tumor therapy?” Physical Review Special Topics - Accelerators and Beams 10 (2007), 10.1103/physrevstab.10.094801.
  * Favalli  _et al._ [2019] A. Favalli, N. Guler, D. Henzlova, S. Croft, K. Falk, D. C. Gautier, K. D. Ianakiev, M. Iliev, S. Palaniyappan, M. Roth, J. C. Fernandez, and M. T. Swinhoe, “Characterizing laser-plasma ion accelerators driving an intense neutron beam via nuclear signatures,” Scientific Reports 9 (2019), 10.1038/s41598-019-39054-z.
  * Kleinschmidt  _et al._ [2018] A. Kleinschmidt, V. Bagnoud, O. Deppert, A. Favalli, S. Frydrych, J. Hornung, D. Jahn, G. Schaumann, A. Tebartz, F. Wagner, G. Wurden, B. Zielbauer, and M. Roth, “Intense, directed neutron beams from a laser-driven neutron source at PHELIX,” Physics of Plasmas 25, 053101 (2018).
  * Roth  _et al._ [2013] M. Roth, D. Jung, K. Falk, N. Guler, O. Deppert, M. Devlin, A. Favalli, J. Fernandez, D. Gautier, M. Geissel, R. Haight, C. E. Hamilton, B. M. Hegelich, R. P. Johnson, F. Merrill, G. Schaumann, K. Schoenberg, M. Schollmeier, T. Shimada, T. Taddeucci, J. L. Tybo, F. Wagner, S. A. Wender, C. H. Wilde, and G. A. Wurden, “Bright laser-driven neutron source based on the relativistic transparency of solids,” Physical Review Letters 110 (2013), 10.1103/physrevlett.110.044802.
  * Aymar  _et al._ [2020] G. Aymar, T. Becker, S. Boogert, M. Borghesi, R. Bingham, C. Brenner, P. N. Burrows, O. C. Ettlinger, T. Dascalu, S. Gibson, T. Greenshaw, S. Gruber, D. Gujral, C. Hardiman, J. Hughes, W. G. Jones, K. Kirkby, A. Kurup, J.-B. Lagrange, K. Long, W. Luk, J. Matheson, P. McKenna, R. McLauchlan, Z. Najmudin, H. T. Lau, J. L. Parsons, J. Pasternak, J. Pozimski, K. Prise, M. Puchalska, P. Ratoff, G. Schettino, W. Shields, S. Smith, J. Thomason, S. Towe, P. Weightman, C. Whyte, and R. Xiao, “LhARA: The laser-hybrid accelerator for radiobiological applications,” Frontiers in Physics 8 (2020), 10.3389/fphy.2020.567738.
  * Zimmer [2020] M. Zimmer,  _Laser-Driven Neutron Sources - A Compact Approach to Non-Destructive Material Analysis_ , Ph.D. thesis, Technische Universität Darmstadt, Darmstadt (2020).
  * Huebl  _et al._ [2020] A. Huebl, M. Rehwald, L. Obst-Huebl, T. Ziegler, M. Garten, R. Widera, K. Zeil, T. E. Cowan, M. Bussmann, U. Schramm, and T. Kluge, “Spectral control via multi-species effects in PW-class laser-ion acceleration,” Plasma Physics and Controlled Fusion 62, 124003 (2020).
  * Berglund  _et al._ [1998] M. Berglund, L. Rymell, H. M. Hertz, and T. Wilhein, “Cryogenic liquid-jet target for debris-free laser-plasma soft x-ray generation,” Review of Scientific Instruments 69, 2361–2364 (1998).
  * Wieland  _et al._ [2001] M. Wieland, T. Wilhein, M. Faubel, C. Ellert, M. Schmidt, and O. Sublemontier, “EUV and fast ion emission from cryogenic liquid jet target laser-generated plasma,” Applied Physics B 72, 591–597 (2001).
  * George  _et al._ [2019] K. M. George, J. T. Morrison, S. Feister, G. K. Ngirmang, J. R. Smith, A. J. Klim, J. Snyder, D. Austin, W. Erbsen, K. D. Frische, J. Nees, C. Orban, E. A. Chowdhury, and W. M. Roquemore, “High-repetition-rate $(\geq$khz$)$ targets and optics from liquid microjets for high-intensity laser–plasma interactions,” High Power Laser Science and Engineering 7 (2019), 10.1017/hpl.2019.35.
  * Djordjević  _et al._ [2021] B. Djordjević, A. Kemp, J. Kim, R. Simpson, S. Wilks, T. Ma, and D. Mariscal, “Modeling laser-driven ion acceleration with deep learning,” Physics of Plasmas 28, 043105 (2021).
  * Zimmer  _et al._ [2021] M. Zimmer, S. Scheuren, T. Ebert, G. Schaumann, C. Rödel, M. Roth, B. Schmitz, V. Bagnoud, and J. Hornung, “Analysis of laser-proton acceleration experiments for development of empirical scaling laws,” Physical Review E 104, 045210 (2021).
  * Ma  _et al._ [2021] T. Ma, D. Mariscal, R. Anirudh, T. Bremer, B. Z. Djordjevic, T. Galvin, E. Grace, S. Herriot, S. Jacobs, B. Kailkhura, R. Hollinger, J. Kim, S. Liu, J. Ludwig, D. Neely, J. J. Rocca, G. G. Scott, R. A. Simpson, B. S. Spears, T. S. Spinka, K. Swanson, J. J. Thiagarajan, B. V. Essen, S. Wang, S. C. Wilks, G. J. Williams, J. Zhang, M. C. Herrmann, and C. Haefner, “Accelerating the rate of discovery: toward high-repetition-rate HED science,” Plasma Physics and Controlled Fusion 63, 104003 (2021).
  * Djordjevic  _et al._ [2021] B. Z. Djordjevic, A. Kemp, J. Kim, S. Wilks, R. Simpson, T. Ma, and D. Mariscal, “Parameter space exploration of short-pulse laser-driven ion acceleration via ensemble simulations and neural networks,” in  _High Power Lasers for Fusion Research VI_, edited by C. L. Haefner and A. A. Awwal (SPIE, 2021).
  * gsi [2020] “Hpc infrastructure virgo,” https://hpc.gsi.de/virgo (2020).
  * Derouillat  _et al._ [2018] J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech, “Smilei : A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation,” Computer Physics Communications 222, 351–373 (2018).
  * Mora [2003] P. Mora, “Plasma expansion into a vacuum,” Phys. Rev. Lett. 90, 185002 (2003).
  * Alejo  _et al._ [2014] A. Alejo, S. Kar, H. Ahmed, A. G. Krygier, D. Doria, R. Clarke, J. Fernandez, R. R. Freeman, J. Fuchs, A. Green, J. S. Green, D. Jung, A. Kleinschmidt, C. L. S. Lewis, J. T. Morrison, Z. Najmudin, H. Nakamura, G. Nersisyan, P. Norreys, M. Notley, M. Oliver, M. Roth, J. A. Ruiz, L. Vassura, M. Zepf, and M. Borghesi, “Characterisation of deuterium spectra from laser driven multi-species sources by employing differentially filtered image plate detectors in thomson spectrometers,” Review of Scientific Instruments 85, 093303 (2014).
  * Obst  _et al._ [2017] L. Obst, S. Göde, M. Rehwald, F.-E. Brack, J. Branco, S. Bock, M. Bussmann, T. E. Cowan, C. B. Curry, F. Fiuza, M. Gauthier, R. Gebhardt, U. Helbig, A. Huebl, U. Hübner, A. Irman, L. Kazak, J. B. Kim, T. Kluge, S. Kraft, M. Loeser, J. Metzkes, R. Mishra, C. Rödel, H.-P. Schlenvoigt, M. Siebold, J. Tiggesbäumker, S. Wolter, T. Ziegler, U. Schramm, S. H. Glenzer, and K. Zeil, “Efficient laser-driven proton acceleration from cylindrical and planar cryogenic hydrogen jets,” Scientific Reports 7 (2017), 10.1038/s41598-017-10589-3.
  * Cappa  _et al._ [2005] C. D. Cappa, W. S. Drisdell, J. D. Smith, R. J. Saykally, and R. C. Cohen, “Isotope fractionation of water during evaporation without condensation,” The Journal of Physical Chemistry B 109, 24391–24400 (2005).
  * smi [2023] “Smilei 4.7 - write a namelist,”  (2023).
  * Lécz [2013] Z. Lécz,  _Laser ion acceleration from a double-layer metal foil_ , Ph.D. thesis, Technische Universität, Darmstadt (2013).
  * Yoo, Jette, and Grondona [2003] A. B. Yoo, M. A. Jette, and M. Grondona, “Slurm: Simple linux utility for resource management,” in  _Workshop on job scheduling strategies for parallel processing_ (Springer, 2003) pp. 44–60.
  * Lécz, Boine-Frankenheim, and Kornilov [2013] Z. Lécz, O. Boine-Frankenheim, and V. Kornilov, “Target normal sheath acceleration for arbitrary proton layer thickness,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 727, 51–58 (2013).
  * Sinigardi, Babaei, and Turchetti [2018] S. Sinigardi, J. Babaei, and G. Turchetti, “TNSA proton maximum energy laws for 2d and 3d PIC simulations,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 909, 438–440 (2018).
  * Chollet  _et al._ [2015] F. Chollet _et al._ , “Keras,” https://keras.io (2015).
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,”  (2015), software available from tensorflow.org.
  * O’Malley  _et al._ [2019] T. O’Malley, E. Bursztein, J. Long, F. Chollet, H. Jin, L. Invernizzi, _et al._ , “Kerastuner,” https://github.com/keras-team/keras-tuner (2019).
  * Bishop [2006] C. M. Bishop,  _Pattern Recognition and Machine Learning_ (Springer-Verlag New York Inc., 2006).
  * VEGA [2022] C. P. L. S. VEGA,  _https://www.clpu.es/en/VEGA_Introduction_ (2022).
  * Schramm  _et al._ [2017] U. Schramm, M. Bussmann, A. Irman, M. Siebold, K. Zeil, D. Albach, C. Bernert, S. Bock, F. Brack, J. Branco, J. Couperus, T. Cowan, A. Debus, C. Eisenmann, M. Garten, R. Gebhardt, S. Grams, U. Helbig, A. Huebl, T. Kluge, A. Köhler, J. Krämer, S. Kraft, F. Kroll, M. Kuntzsch, U. Lehnert, M. Loeser, J. Metzkes, P. Michel, L. Obst, R. Pausch, M. Rehwald, R. Sauerbrey, H. Schlenvoigt, K. Steiniger, and O. Zarini, “First results with the novel petawatt laser acceleration facility in dresden,” Journal of Physics: Conference Series 874, 012028 (2017).
  * Spears  _et al._ [2018] B. K. Spears, J. Brase, P.-T. Bremer, B. Chen, J. Field, J. Gaffney, M. Kruse, S. Langer, K. Lewis, R. Nora, J. L. Peterson, J. J. Thiagarajan, B. V. Essen, and K. Humbird, “Deep learning: A guide for practitioners in the physical sciences,” Physics of Plasmas 25, 080901 (2018).
  * Wang and Deng [2018] M. Wang and W. Deng, “Deep visual domain adaptation: A survey,” Neurocomputing 312, 135–153 (2018).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods 17, 261–272 (2020).
  * Lalee, Nocedal, and Plantenga [1998] M. Lalee, J. Nocedal, and T. Plantenga, “On the implementation of an algorithm for large-scale equality constrained optimization,” SIAM Journal on Optimization 8, 682–706 (1998).
  * Sobol$\prime$ [2001] I. Sobol$\prime$, “Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates,” Mathematics and Computers in Simulation 55, 271–280 (2001).
  * Saltelli [2002] A. Saltelli, “Making best use of model evaluations to compute sensitivity indices,” Computer Physics Communications 145, 280–297 (2002).
  * Saltelli  _et al._ [2010] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola, “Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index,” Computer Physics Communications 181, 259–270 (2010).
  * Pianosi and Wagener [2015] F. Pianosi and T. Wagener, “A simple and efficient method for global sensitivity analysis based on cumulative distribution functions,” Environmental Modelling & Software 67, 1–11 (2015).
  * Baroni and Francke [2020] G. Baroni and T. Francke, “An effective strategy for combining variance- and distribution-based global sensitivity analysis,” Environmental Modelling & Software 134, 104851 (2020).
  * Ghanem and Spanos [1991] R. G. Ghanem and P. D. Spanos, “Stochastic finite element method: Response statistics,” in  _Stochastic Finite Elements: A Spectral Approach_ (Springer New York, 1991) pp. 101–119.
  * Schmitz, Metternich, and Boine-Frankenheim [2022] B. Schmitz, M. Metternich, and O. Boine-Frankenheim, “Automated reconstruction of the initial distribution of laser accelerated ion beams from radiochromic film (RCF) stacks,” Review of Scientific Instruments 93, 093306 (2022).
  * Leshno  _et al._ [1993] M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function,” Neural networks 6, 861–867 (1993).
  * Birdsall and Langdon [2004] C. Birdsall and A. Langdon,  _Plasma Physics via Computer Simulation_ , Series in Plasma Physics and Fluid Dynamics (Taylor & Francis, 2004).
  * Mulser and Bauer [2010] P. Mulser and D. Bauer,  _High power laser-matter interaction_ , Vol. 238 (Springer, 2010).
  * Arber  _et al._ [2015] T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, “Contemporary particle-in-cell approach to laser-plasma modelling,” Plasma Physics and Controlled Fusion 57, 113001 (2015).
  * Buckingham [1914] E. Buckingham, “On physically similar systems & illustrations of the use of dimensional equations,” Physical Review 4, 345–376 (1914).
  * Bagnoud  _et al._ [2009] V. Bagnoud, B. Aurand, A. Blazevic, S. Borneis, C. Bruske, B. Ecker, U. Eisenbarth, J. Fils, A. Frank, E. Gaul, S. Goette, C. Haefner, T. Hahn, K. Harres, H.-M. Heuck, D. Hochhaus, D. H. H. Hoffmann, D. Javorková, H.-J. Kluge, T. Kuehl, S. Kunzer, M. Kreutz, T. Merz-Mantwill, P. Neumayer, E. Onkels, D. Reemts, O. Rosmej, M. Roth, T. Stoehlker, A. Tauschwitz, B. Zielbauer, D. Zimmer, and K. Witte, “Commissioning and early experiments of the PHELIX facility,” Applied Physics B 100, 137–150 (2009).
  * Busold  _et al._ [2014] S. Busold, A. Almomani, V. Bagnoud, W. Barth, S. Bedacht, A. Blažević, O. Boine-Frankenheim, C. Brabetz, T. Burris-Mog, T. Cowan, O. Deppert, M. Droba, H. Eickhoff, U. Eisenbarth, K. Harres, G. Hoffmeister, I. Hofmann, O. Jaeckel, R. Jaeger, M. Joost, S. Kraft, F. Kroll, M. Kaluza, O. Kester, Z. Lecz, T. Merz, F. Nürnberg, H. Al-Omari, A. Orzhekhovskaya, G. Paulus, J. Polz, U. Ratzinger, M. Roth, G. Schaumann, P. Schmidt, U. Schramm, G. Schreiber, D. Schumacher, T. Stoehlker, A. Tauschwitz, W. Vinzenz, F. Wagner, S. Yaramyshev, and B. Zielbauer, “Shaping laser accelerated ions for future applications – the LIGHT collaboration,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 740, 94–98 (2014).
  * Bourdier [1983] A. Bourdier, “Oblique incidence of a strong electromagnetic wave on a cold inhomogeneous electron plasma. relativistic effects,” Physics of Fluids 26, 1804 (1983).
  * Gibbon  _et al._ [1999] P. Gibbon, A. Andreev, E. Lefebvre, G. Bonnaud, H. Ruhl, J. Delettrez, and A. R. Bell, “Calibration of one-dimensional boosted kinetic codes for modeling high-intensity laser–solid interactions,” Physics of Plasmas 6, 947–953 (1999).
  * Cui  _et al._ [2013] Y.-Q. Cui, W.-M. Wang, Z.-M. Sheng, Y.-T. Li, and J. Zhang, “Laser absorption and hot electron temperature scalings in laser–plasma interactions,” Plasma Physics and Controlled Fusion 55, 085008 (2013).
  * Gibbon [2005] P. Gibbon,  _Short Pulse Laser Interactions with Matter_ (IMPERIAL COLLEGE PRESS, 2005).
  * Savitzky and Golay [1964] A. Savitzky and M. J. E. Golay, “Smoothing and differentiation of data by simplified least squares procedures.” Analytical Chemistry 36, 1627–1639 (1964), https://doi.org/10.1021/ac60214a047 .