# Simulation of ion temperature gradient driven modes with 6D kinetic Vlasov
code

M. Raeth<EMAIL_ADDRESS>Max Planck Institute for Plasma Physics,
Boltzmannstr. 2, 85748 Garching, Germany  K. Hallatschek  Max Planck Institute
for Plasma Physics, Boltzmannstr. 2, 85748 Garching, Germany  K. Kormann
Ruhr-Universität Bochum, Universitätsstraße 150, 44780 Bochum, Germany

###### Abstract

With the increase in computational capabilities over the last years it becomes
possible to simulate more and more complex and accurate physical models.
Gyrokinetic theory has been introduced in the 1960s and 1970s in the need of
describing a plasma with more accurate models than fluid equations, but
eliminating the complexity of the fast gyration about the magnetic field
lines. Although results from current gyrokinetic computer simulations are in
fair agreement with experimental results in core physics, crucial assumptions
made in the derivation make it unreliable in regimes of higher fluctuations
and stronger gradient, such as the tokamak edge. With our novel optimized and
scalable semi-Lagrangian solver we are able to simulate ion-temperature
gradient modes with the 6D kinetic model including the turbulent saturation.
After thoroughly testing our simulation code against analytical computations
and gyrokinetic simulations (with the gyrokinetic code GYRO), it has been
possible to show first plasma properties that go beyond standard gyrokinetic
simulations. This includes the explicit description of the complete
perpendicular energy fluxes and the excitation of high frequency waves (around
the Larmor frequency) in the nonlinear saturation phase.

††preprint: AIP/123-QED

##  I Introduction

Kinetic models are capable of describing physical phenomena in a tokamak
plasma on all scales from the size of the device ($\sim$1m) down to
microscales of the order of the electron Larmor radius ($\sim 10^{-4}$m).
However, until recent years, the computational capabilities did not allow the
computation of the full 6D-kinetic equation for a magnetized plasma with the
time resolution required to resolve the fast orbital motion of the articles
around the magnetic field. To circumvent this problem, many gyrokinetic models
have been developed littlejohn_hamiltonian_1981 , which reduce the
dimensionality and the required time resolution by eliminating frequencies of
the order of the Larmor frequency $\omega_{c}$.
Gyrokinetic simulations codes such as GENE jenko_massively_2000 , GS2
kotschenreuther_comparison_1995 , GYRO candy_high-accuracy_2016  or ORB5
jolliet_global_2007  are in good agreement with experiments in the core of
fusion devices where only small perturbation amplitudes and gradients (in
density and temperature) are present. Nonetheless, in regimes of large
gradients and high fluctuation amplitudes, such as the plasma edge of a
tokamak, the gyrokinetic approximations are debatable and at least those
models based on a $\delta f$ approximation break down completely. Moreover,
the limited frequency range in gyrokinetic simulations precludes the accurate
representation of physical processes associated with high-frequency modes. To
illustrate this point, Craddock et al. demonstrated the suppression of
turbulence by ion Bernstein waves in their study craddock_theory_1994 .
With the rise of modern high-performance computing capabilities, it becomes
more and more viable to simulate the full 6D kinetic equation. For this
purpose, the massively parallel semi-Lagrangian code BSL6D for the Vlasov
equation in 6D phase space has been developed kormann_massively_2019 . The
code, BSL6D, is capable of simulating plasma turbulence across a broad range
of frequencies, including those beyond the Larmor frequency, and with
arbitrarily large fluctuation amplitudes.
Past efforts to venture into regimes beyond gyrokinetic theory include the
drift-cyclotron model by Waltz and Deng waltz_nonlinear_2013  and a non-
gyrokinetic magnetized plasma turbulence code developed by Sturdevant et al.
sturdevant_implicit_2016 ; sturdevant_low_2017 ; miecnikowski_nonlinear_2018 .
However, they either still use a reduced kinetic model, or are limited to a
single toroidal model, effectively resulting in a 5D simulation. This work
shows unprecedented 6D nonlinear kinetic simulations of the slab ITG
instability, examining the energy fluxes inherent in the model.
In Section II, we provide a concise introduction to the physical model
implemented in the code for slab Ion Temperature Gradient (ITG) simulations.
Subsequently, we conduct a verification study in the linear regime,
encompassing a comparison with the analytically derived dispersion relation
(section 2) and the computation of quasi-linear energy fluxes (section 1).
Finally, in section V, we present results from simulations capturing the
nonlinear saturation phase of the ITG instability, comparing them to results
obtained from the gyrokinetic code CGYRO. Intriguingly, the 6D kinetic
simulations reveal that the saturation phase leads to the excitation of high-
frequency ion Bernstein waves. We conducted our nonlinear studies employing
both a Boussinesq approximation for the gradients (corresponding to the
gyrokinetic description of the gradients) and a fully nonlinear treatment.
This allows us to investigate the saturation phase under different conditions
and gain a deeper understanding of the underlying physics. The discovery of
high-frequency ion Bernstein waves during the saturation phase challenges
conventional understanding and opens new avenues for research into the
nonlinear dynamics of plasma turbulence.

##  II 6D semi-Lagrangian kinetic turbulence code for magnetized plasmas

The Vlasov equation describes the motion of a plasma in presence of
electromagnetic fields. We consider ions in a constant and homogeneous
magnetic field $\bm{B}=\hat{\bm{z}}$ represented by a 6D distribution function
$f$, with an electric field $\bm{E}$ originating from the interaction with
adiabatic electrons. The target of our 6D simulations is the ion kinetic
equation in dimensionless variables ($\rho_{i}\;(\text{Ion Larmor
radius})=v_{\text{th}}\;(\text{Ion thermal velocity})=n\;(\text{Ion background
density})=T\;(\text{Ion background temperature})=1$)

$\partial_{t}f+\bm{v}\cdot\nabla
f+(-\nabla\phi+\bm{v}\times\bm{\hat{}}{\bm{z}})\cdot\nabla_{\bm{v}}f=0.$ (1)

For the simulation of ion of temperature gradient (ITG) driven modes, the
gradient is introduced in Boussinesq approximation. For this purpose, we
assume that a background distribution function $f_{0}$ exists, which fulfills

$\displaystyle\bm{v}\cdot\nabla
f_{0}+\bm{v}\times\hat{\bm{z}}\cdot\nabla_{v}f_{0}=0.$ (2)

This condition is met, when the background distribution is parameterized by an
arbitrary function $g$ with
$f_{0}(\bm{r},\bm{v})=f_{0}(\bm{R}+\bm{\rho},\bm{v})=g(\bm{R},v_{\perp},v_{z})$
(with $v_{z}=\bm{v}\cdot\hat{\bm{z}}$, $v_{\perp}=|\bm{v}\times\hat{\bm{z}}|$
and $\bm{R}=\bm{r}-\bm{\rho}=\bm{r}+\bm{v}\times\bm{\hat{}}{\bm{z}}$ is the
location of the gyrocenter of a given particle, where $\bm{r}$ is the
configurations space coordinate and $\bm{\rho}=\hat{\bm{z}}\times\bm{v}$ is
the Larmor radius vector). The background distribution is chosen to resemble a
Maxwellian with varying temperature $T(\bm{R})$, and thus reads

$\displaystyle g_{0}(\bm{R},v_{\perp},v_{z}):=\left(\frac{1}{2\pi
T(\bm{R})}\right)^{\frac{3}{2}}\exp\left(-\frac{[v_{\perp}^{2}+v_{z}^{2}]}{2T(\bm{R})}\right).$
(3)

When splitting the distribution function in the Vlasov equation into a
perturbation $\delta f=f-g_{0}$ and the background $g_{0}$, one obtains

$\displaystyle\partial_{t}\delta f+\bm{v}\cdot\nabla\delta
f+\left[-\nabla\phi+(\bm{v}\times\hat{\bm{z}})\right]\cdot\nabla_{v}\delta
f=\nabla\phi\cdot\nabla_{v}g_{0},$ (4)

where the gradient on the right-hand side is given by

$\displaystyle\nabla_{v}g_{0}$ $\displaystyle=\hat{\bm{z}}\times\nabla
g_{0}+\nabla_{v}g_{0}=:{\bm{v}^{*}}g_{0}+\nabla_{v}g_{0}$
$\displaystyle=-g_{0}\hat{\bm{z}}\times\left(\frac{\nabla
T}{T}\frac{3-(v_{\perp}^{2}+v_{z}^{2})}{2}\right)+\nabla_{v}g_{0}.$ (5)

The limit that the gradient length is infinite $T(\bm{R})\rightarrow T=1$ is
taken, such that the background distribution is independent of $\bm{R}$. Thus,
$g_{0}$ is replaced with the homogeneous Maxwellian distribution
$f_{M}=\frac{1}{(2\pi)^{\frac{3}{2}}}e^{-\frac{v^{2}}{2}}$.
For simplicity, we assume that the electrons are adiabatic
$f_{e}=e^{-\phi}f_{M}\approx f_{M}-\phi f_{M}$ and the plasma is quasineutral
$n_{e}=n$ (with the electron density $n_{e}$), the complete system of equation
is given by

$\displaystyle\partial_{t}f+\bm{v}\cdot\nabla
f+\left(-\nabla\phi+\bm{v}\times\hat{\bm{z}}\right)\cdot\nabla_{v}f=\bm{v}^{*}\cdot\nabla\phi
f_{M},$ (6) $\displaystyle\phi=\delta n,\quad\quad\quad\vspace{2cm}\delta
n=\int h\mathrm{d}^{3}v.$ (7)

For all presented results, the electron temperature has been chosen to be the
same as the ion temperature $T_{e}=1$. The system is integrated in time using
a semi-Lagrangian method sonnendrucker_semi-lagrangian_1999 , which combines
the advantages of grid-based (Eulerian) and particle-based (Lagrangian)
methods. In this method, the distribution function is represented as a
discrete function on a grid. However, instead of explicitly computing partial
derivatives (as in Eulerian methods), the distribution function is advected
along the characteristics of the Vlasov equation. This is achieved by solving
the characteristic equations for particle trajectories and interpolating the
distribution function at their new positions. This approach circumvents the
strict CFL (Courant-Friedrichs-Lewy) condition imposed on Eulerian methods,
allowing for the selection of a larger time step. Additionally, the noise
levels are significantly lower compared to particle methods, making the semi-
Lagrangian method a more robust and efficient choice for simulating the
dynamics of plasmas. Implementation details can be found in the work of
Kormann et al. kormann_massively_2019 .

##  III Dispersion relation

###  1 Analytical description

The response of the non-adiabatic distribution, given by $h=f-f_{M}+\phi
f_{M}$, to a present electrostatic potential $\phi$ in the linearized system
is given by the Gordeev integral gordeev_notitle_1952 . This computation
yields an expression for the velocity distribution
($\bm{v}=(v_{x},v_{y},v_{z})^{T}$) of a specific Fourier mode of the
distribution function $h_{k}$, influenced by an electrostatic wave
$\phi=\phi_{k}e^{i(\bm{k}\cdot\bm{r}-\omega t)}$ with frequency $\omega$ and
wave vector $\bm{k}$ (where $k_{z}=\bm{k}\cdot\hat{\bm{z}}$,
$\bm{k}_{\perp}=\bm{k}-k_{z}\hat{\bm{z}}$, and $\alpha$ is the azimuth angle
in velocity space with $\hat{\bm{z}}\times\bm{k}$ corresponding to
$\alpha=\frac{\pi}{2}$) brambilla_kinetic_1998

$\displaystyle
h_{k}=\phi_{k}(\omega-\bm{k}\cdot\bm{v^{*}})f_{M}\sum_{m,p\in\mathbb{Z}}\frac{J_{m}(k_{\perp}v_{\perp})J_{p}(k_{\perp}v_{\perp})e^{i(p-m)\alpha}}{\omega-
k_{z}v_{z}-p},$ (8)

where $J_{p}(x)$ represents the Bessel function of the first kind, and
$\bm{v}^{*}$ is defined in equation (5). The density response $n_{k}=\int
h_{k}\mathrm{d}^{3}v-\phi_{k}$ is given by

$\displaystyle n_{k}=$
$\displaystyle\int_{-\infty}^{\infty}\int_{0}^{\infty}\int_{0}^{2\pi}i(\omega-\bm{k}\cdot\bm{v}^{*})\phi_{k}f_{M}$
$\displaystyle\sum_{m,p\in\mathds{Z}}\frac{J_{m}(k_{\perp}v_{\perp})J_{p}(k_{\perp}v_{\perp})e^{i(p-m)\alpha}}{k_{z}v_{z}-\omega-p}v_{\perp}\mathrm{d}\alpha\mathrm{d}v_{\perp}\mathrm{d}v_{z}-\phi_{k}.$
(9)

To simplify the integration, the source term $\bm{k}\cdot\bm{v}^{*}f_{M}$,
defined in equation (5) can be written in terms of a derivative with respect
to the helping variable $\xi$ (assuming $\frac{\nabla
T}{T}\propto\hat{\bm{y}}$)

$\displaystyle\bm{k}\cdot\bm{v}^{*}f_{M}(v_{\perp},v_{z})=k_{y}\frac{\nabla
T}{T}\partial_{\xi}f_{M}(\sqrt{\xi}v_{\perp},\sqrt{\xi}v_{z}),$ (10)

which is set $\xi=1$ after the derivative is calculated. The parallel velocity
integral is resolved using the plasma dispersion function

$\displaystyle
Z(x)=\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}\frac{e^{-t^{2}}}{(t-x)^{2}}\mathrm{d}t.$
(11)

After solving the perpendicular integral, we introduce
$\Gamma_{n}(x^{2})=e^{-x^{2}}I_{n}(x^{2})$, resulting in the expression
brambilla_kinetic_1998

$\displaystyle\frac{n_{k}}{\phi_{k}}=\left[\omega-k_{y}\frac{\nabla
T}{T}\partial_{\xi}\right]\frac{\sqrt{\xi}}{|k_{z}|\sqrt{{2}}}\sum_{p\in\mathds{Z}}Z\left(\frac{\omega+p}{|k_{z}|}\sqrt{\frac{\xi}{2}}\right)\Gamma_{n}\left(\frac{k_{\perp}^{2}}{\xi}\right)-1.$
(12)

The density response is introduced in the quasi-neutrality equation which
gives the electrostatic potential form equation (7). The resulting dispersion
relation, given by solution $\omega(\bm{k})$ to the equation

$\displaystyle 0=\left[\omega-k_{y}\frac{\nabla
T}{T}\partial_{\xi}\right]\frac{\sqrt{\xi}}{|k_{z}|\sqrt{{2}}}\sum_{p\in\mathds{Z}}Z\left(\frac{\omega+p}{|k_{z}|}\sqrt{\frac{\xi}{2}}\right)\Gamma_{n}\left(\frac{k_{\perp}^{2}}{\xi}\right)-2,$
(13)

describes all modes for the 6D kinetic system. In the gyrokinetic limit
$\omega\ll 1$ and $k_{z}\ll 1$ all terms of the sum vanish except the $p=0$
contribution, simplifying the dispersion relation to

$\displaystyle 0=\left[\omega-k_{y}\frac{\nabla
T}{T}\partial_{\xi}\right]\frac{\sqrt{\xi}}{|k_{z}|\sqrt{{2}}}Z\left(\frac{\omega}{|k_{z}|}\sqrt{\frac{\xi}{2}}\right)\Gamma_{0}\left(\frac{k_{\perp}^{2}}{\xi}\right)-2.$
(14)

###  2 Verification of dispersion relation

For the numerical verification of the code against the dispersion relation, a
linear simulation has been conducted. The removal of the non-linear term from
the Vlasov equation (1) results in indefinite linear growth, simplifying the
determination of growth rates. In space, a box with dimensions
$N_{x}=128\times 128\times 8$ has been selected. The box size ensures the
smallest perpendicular wave number is $k_{x,0}=k_{y,0}=0.3$
($L_{x}=L_{y}=\frac{20}{3}\pi$). The parallel length is determined by the wave
number of the fastest-growing mode in the system, dependent on the temperature
gradient. A temperature gradient of $\frac{\nabla T}{T}=0.05$ has been
employed, leading to an ideal parallel wave number of $k_{z}=\frac{1}{240}$
and, consequently, $L_{z}=480\pi$. The velocity space is symmetric in all
directions, with $v_{\mathrm{max}}=4$ and $N_{v}=32\times 32\times 32$. For
interpolation, a 7th-order Lagrange interpolation has been applied to the
velocity directions, and an 8th order for the spatial directions.

The simulation commences with a small white noise density perturbation, and
growth rates are determined by fitting the time evolution of various Fourier
modes of the density.

Figure 1: Comparison of the frequency (left) and the growth rate (right) for
various modes with $k_{z}=\frac{1}{240}$ determined from linear simulations
($\bullet$) and computed by numerical root finding from the analytical
dispersion relation (13) in solid lines

Figure 1 presents a comparison between the numerical results and the
analytical dispersion relation. The code demonstrates excellent agreement for
both the growth rate and frequency. However, discrepancies are more pronounced
for very small and larger wave numbers. For small $k_{y}$, the growth rates
are minimal, or even negative, making the numerical computation challenging.
As for larger perpendicular wave numbers, the numerical damping, attributed to
the inherent diffusion of the Lagrange interpolation, is comparable to the ITG
growth rates. Consequently, the code tends to underestimate the growth rate.
This discrepancy can be mitigated by either enhancing the resolution or
adopting a higher interpolation order rath_beyond_2023 .

##  IV Contributions to the energy flux in the 6D system

###  1 Analytical description

The total energy density $E$ of the system results from the sum of the kinetic
energy of ions and the energy stored in the electrons

$\displaystyle E$
$\displaystyle=\int\frac{v^{2}}{2}f\mathrm{d}^{3}v+\frac{n^{2}}{2}.$ (15)

The time derivative of the total energy is computed using the Vlasov equation
(1)

$\displaystyle\partial_{t}E=$
$\displaystyle\int\frac{v^{2}}{2}\partial_{t}f\mathrm{d}^{3}v+n\partial_{t}n$
$\displaystyle=$ $\displaystyle-\int\frac{v^{2}}{2}\bm{v}\cdot\nabla
f\mathrm{d}^{3}v$
$\displaystyle-\int\frac{v^{2}}{2}(-\nabla\phi+\bm{v}\times\hat{\bm{z}})\cdot\nabla_{\bm{v}}f\mathrm{d}^{3}v$
$\displaystyle+n\partial_{t}n.$ (16)

Applying integration by parts to the second integral yields

$\displaystyle\partial_{t}E$
$\displaystyle=-\int\frac{v^{2}}{2}\bm{v}\cdot\nabla
f\mathrm{d}^{3}v-\int\bm{v}\cdot(\nabla\phi)f\mathrm{d}^{3}v+n\partial_{t}n$
(17)

which farther can be modified to be expressed as a divergence of the entire
integral

$\displaystyle\partial_{t}E$
$\displaystyle=-\nabla\cdot\left(\int\frac{v^{2}}{2}\bm{v}f\mathrm{d}^{3}v\right)-\nabla\cdot\left(\phi\int\bm{v}f\mathrm{d}^{3}v\right)$
$\displaystyle+\phi\int\bm{v}\cdot\nabla f\mathrm{d}^{3}v+n\partial_{t}n.$
(18)

By introducing
$\partial_{t}n=-\nabla\cdot\Gamma:=-\nabla\cdot\int\bm{v}f\mathrm{d}^{3}v$ and
utilizing the field equation (7), the expression can be succinctly summarized
as

$\displaystyle\partial_{t}E$
$\displaystyle=-\nabla\cdot\bm{Q}-\nabla\cdot(\phi\bm{\Gamma}).$ (19)

Apart from $\bm{Q}=\frac{1}{2}\int\bm{v}v^{2}f\mathrm{d}^{3}v$, an additional
energy flux $\bm{S}=\phi\bm{\Gamma}$ appears which, can be identified as a
Poynting like energy flux. The energy flux in the gyrokinetic system is
completely described by the $\bm{E}\times\bm{B}$ heat flux. Thus, in the
gyrokinetic limit, the energy flux in the 6D system has to be equal the
$\bm{E}\times\bm{B}$ heat flux. The various contributions to the energy flux
$\bm{Q}=\int\frac{v^{2}}{2}f\mathrm{d}^{3}v$ are obtained by computing the
time derivative with the Vlasov equation (1). The Lorentz force term of the
Vlasov equation is altered by partial integration

$\displaystyle\int$
$\displaystyle\frac{v^{2}}{2}\bm{v}\left[(-\nabla\phi+\bm{v}\times\hat{\bm{z}})\cdot\nabla_{\bm{v}}f\right]\mathrm{d}^{3}v$
$\displaystyle=$
$\displaystyle\int\left(\nabla\phi\cdot\bm{v}\right)\bm{v}f\mathrm{d}^{3}v+\int\frac{v^{2}}{2}\nabla\phi
f\mathrm{d}^{3}v-\int(\bm{v}\times\hat{\bm{z}})\frac{v^{2}}{2}f\mathrm{d}^{3}v$
(20)

After introducing the expression for kinetic energy density
$\epsilon=\int\frac{v^{2}}{2}f\mathrm{d}^{3}v$, the stress tensor
$\bm{\mathrm{\Pi}}=\int\bm{v}\bm{v}f\mathrm{d}^{3}v$,
$\bm{\mathrm{\Pi}}^{*}=\int\frac{v^{2}}{2}\bm{v}\bm{v}f\mathrm{d}^{3}v$ and
computing the cross product with the magnetic field
$\partial_{t}\bm{Q}\times\hat{\bm{z}}$, the Grassmann identify can be applied
to obtain an expression for the perpendicular energy flux

$\displaystyle\partial_{t}{\bm{Q}}\times\hat{\bm{z}}=-$
$\displaystyle\bm{Q}_{\perp}-(\nabla\cdot\bm{\mathrm{\Pi}}^{*})\times\hat{\bm{z}}$
$\displaystyle-(\bm{\mathrm{\Pi}}\cdot\nabla\phi)\times\hat{\bm{z}}-\nabla\phi\times\hat{\bm{z}}\epsilon.$
(21)

For low frequency waves $\omega\ll 1$, the time derivative can be neglected
$\partial_{t}\bm{Q}\approx 0$ and without boundary contributions, the flux
surface ($y-z-$plane,
$\langle\cdot\rangle:=\frac{1}{L_{y}L_{z}}\int\cdot\mathrm{d}y\mathrm{d}z$)
average results in

$\displaystyle\langle\bm{Q}_{\perp}\rangle={-\langle\nabla\phi\times\hat{\bm{z}}\epsilon\rangle}-\langle(\bm{\mathrm{\Pi}}\cdot\nabla\phi)\times\hat{\bm{z}}\rangle.$
(22)

The overall change in energy density can be written as

$\displaystyle\langle\partial_{t}E\rangle=\nabla\cdot\left(\langle\nabla\phi\times\hat{\bm{z}}\epsilon\rangle+\langle(\bm{\mathrm{\Pi}}\cdot\nabla\phi)\times\hat{\bm{z}}\rangle-\langle\bm{S}\rangle\right).$
(23)

In conclusion, the mean perpendicular energy flux consists of three
contributions

  1. 1.

Poynting flux $\bm{S}=\phi\bm{\Gamma}$

  2. 2.

$\bm{E}\times\bm{B}$ \- energy flux
$\bm{Q}^{\bm{E}\times\bm{B}}=-(\nabla\phi\times\hat{\bm{z}})\epsilon$

  3. 3.

Stress tensor induced energy flux
$\bm{Q}^{\bm{\mathrm{\Pi}}}=(\bm{\mathrm{\Pi}}\cdot\nabla\phi)\times\hat{\bm{z}}$

This depiction of the energy flux offers two significant advantages. Firstly,
it provides a direct mean of calculating the individual contributions to the
energy flux from the code. In the gyrokinetic limit, the Poynting flux
$\bm{S}$ and the stress-induced heat flux $\bm{Q}^{\bm{\mathrm{\Pi}}}$ nullify
each other, leaving the $\bm{E}\times\bm{B}$-heat flux as the sole remaining
form of energy transport. The equivalency between these two energy fluxes
serves as a useful test for verifying the accuracy of the energy flux in the
gyrokinetic regime. Secondly, it facilitates the direct computation of
contributions to the energy flux, rendering it less susceptible to gyro-
oscillations and, consequently, easier to calculate.
In the following, we analytically compute the contributions to the energy flux
for a given solution of the dispersion relation.

#### 1\. The Poynting flux

$\bm{S}=\phi\bm{\Gamma}$ is given by the particle flux multiplied by the
electrostatic potential. The derivation of the particle flux, denoted by
$\bm{\Gamma}=\int\bm{v}f\mathrm{d}^{3}v$, follows the same approach as for the
density response in equation (12) using the ansatz for the distribution
function from equation (8).
For the perpendicular flux, computations are carried out independently in the
two directions $(v_{1},v_{2})=v_{\perp}(-\sin\alpha,\cos\alpha)$, where
$\alpha$ is the azimuth angle in velocity space, ensuring that the angle of
$\hat{\bm{z}}\times\bm{k}$ corresponds to $\alpha=\frac{\pi}{2}$. For a more
general calculation, the integral is evaluated for an arbitrary complex
Fourier mode $P_{l}=\int v_{\perp}e^{ila}f\mathrm{d}^{3}v$. We recall the
definition of the non-adiabatic perturbation of the distribution function
$h=f-f_{M}+\phi f_{M}$, which allows us to compute the integral $P_{l}$ from
equation (8)

$\displaystyle P_{l}$ $\displaystyle=\int
v_{\perp}h_{k}e^{il\alpha}\mathrm{d}^{3}v$
$\displaystyle=\int\left[\omega-\bm{k}\cdot\bm{v}^{*}\right]\phi_{k}f_{M}$
$\displaystyle\hskip
7.11317pt\times\sum_{m,p\in\mathds{Z}}\frac{J_{m}(k_{\perp}v_{\perp})J_{p}(k_{\perp}v_{\perp})e^{i(p-m)\alpha}}{\omega-
k_{z}v_{z}-p}v_{\perp}e^{il\alpha}\mathrm{d}^{3}v,$ (24)

and subsequently assemble the particle flux $\bm{\Gamma}$. The double sum
simplifies when the angle integral is evaluated with
$\int_{0}^{2\pi}exp({i(p-m+l)\alpha})\mathrm{d}\alpha=2\pi\delta_{m,p+l}$
which leads to the expression

$\displaystyle P_{l}=2\pi\phi_{k}$
$\displaystyle\int_{-\infty}^{\infty}\int_{0}^{\infty}\left[\omega-\bm{k}\cdot\bm{v}^{*}\right]f_{M}$
$\displaystyle\sum_{p\in\mathds{Z}}\left[\frac{J_{p+l}(v_{\perp}k_{\perp})J_{p}(v_{\perp}k_{\perp})}{(\omega-
k_{z}v_{z}-p)}\right]v_{\perp}^{2}\mathrm{d}v_{\perp}\mathrm{d}v_{z}.$ (25)

The velocity integrals are solved numerically for a given set of parameters
and the respective solution ($\omega,\bm{k}$) of the dispersion relation (13).
When computing the flux in a spatial direction, the definition of the velocity
angle $\alpha$ has to be kept in mind. The angle is defined such that
$\bm{k}\cdot\hat{\bm{z}}$ corresponds to $\alpha=\frac{\pi}{2}$. Thus, the
fluxes

$\displaystyle\Gamma_{1}$
$\displaystyle=\frac{i}{2}\left(P_{1}-P_{-1}\right),$ (26)
$\displaystyle\Gamma_{2}$
$\displaystyle=\frac{1}{2}\left(P_{-1}+P_{1}\right),$ (27)

point in the direction perpendicular (for $\Gamma_{1}$) and parallel (for
$\Gamma_{2}$) of the wave vector. In the case of $\bm{k}=k\hat{\bm{y}}$,
$\Gamma_{1}$ and $\Gamma_{2}$ corresponds to the $x$ and $y$ direction
respectively. For a general wave vector, the fluxes need to be combined to
compute the correct fluxes in $x$ and $y$ direction. The resulting velocity
integrals

$\displaystyle\frac{\Gamma_{1}}{\phi_{k}}=\Xi i\pi\int_{-\infty}^{\infty}$
$\displaystyle\int_{0}^{\infty}\mathrm{d}v_{\perp}\mathrm{d}v_{z}v_{\perp}^{2}f_{M}(\sqrt{\xi}v_{z},\sqrt{\xi}v_{\perp})$
$\displaystyle\times\left[\sum_{p\in\mathds{Z}}\left(\frac{J_{p+1}(v_{\perp}k_{\perp})J_{p}(v_{\perp}k_{\perp})}{\omega-
k_{z}v_{z}-p}\right.\right.$ $\displaystyle\hskip
28.45274pt\left.\left.-\frac{J_{p-1}(v_{\perp}k_{\perp})J_{p}(v_{\perp}k_{\perp})}{\omega-
k_{z}v_{z}-p}\right)\right],$
$\displaystyle\frac{\Gamma_{2}}{\phi_{k}}=\Xi\pi\int_{-\infty}^{\infty}$
$\displaystyle\int_{0}^{\infty}\mathrm{d}v_{\perp}\mathrm{d}v_{z}v_{\perp}^{2}f_{M}(\sqrt{\xi}v_{z},\sqrt{\xi}v_{\perp})$
$\displaystyle\times\left[\sum_{p\in\mathds{Z}}\left(\frac{J_{p-1}(v_{\perp}k_{\perp})J_{p}(v_{\perp}k_{\perp})}{\omega-
k_{z}v_{z}-p}\right.\right.$ $\displaystyle\hskip
28.45274pt\left.\left.+\frac{J_{p+1}(v_{\perp}k_{\perp})J_{p}(v_{\perp}k_{\perp})}{\omega-
k_{z}v_{z}-p}\right)\right],$ (28)

where $\Xi=\left[\omega-k_{y}\left(\frac{\nabla
T}{T}\partial_{\xi}\right)\right]$, are then computed numerically.

#### 2\. The $\bm{E}\times\bm{B}$ heat flux

$\bm{Q}_{\perp}^{E\times B}={-\nabla\phi\times\hat{\bm{z}}\epsilon}$ can be
computed analogously to the Poynting flux. In Fourier space, the flux can be
written as
${{\bm{Q}}_{\perp}^{\bm{E}\times\bm{B}}=i\phi\bm{k}\times\hat{\bm{z}}\epsilon_{k}}$,
where the energy density response
$\epsilon_{k}=\int\frac{v^{2}}{2}h_{k}\mathrm{d}^{3}v$ is computed from the
distribution function (8) resulting in

$\displaystyle\epsilon_{k}=\phi_{k}\left[\omega\partial_{\xi}\Psi-\hat{\bm{z}}\times\frac{\nabla
T}{T}\partial_{\xi}\left(\xi^{\frac{3}{2}}\partial_{\xi}\Psi\right)\right],\text{
with }$ (29) $\displaystyle\hskip
42.67912pt\Psi=\frac{1}{|k_{z}|}\sqrt{\frac{\xi}{2}}e^{-\frac{k_{\perp}^{2}}{\xi}}\sum_{p\in\mathds{Z}}Z\left(\frac{\omega+p}{|k_{z}|}\sqrt{\frac{\xi}{2}}\right)I_{p}\left(\frac{k_{\perp}^{2}}{\xi}\right).$

#### 3\. The stress induced energy flux

is computed from the stress tensor response
$\bm{\mathrm{\Pi}}_{k}=\int\bm{v}\bm{v}h_{k}\mathrm{d}^{3}v$. When the
velocity components
$(v_{1},v_{2},v_{z})=(-v_{\perp}\sin\alpha,v_{\perp}\cos\alpha,v_{z})$ are
introduced, the expression for the energy flux component reads

$\displaystyle{\bm{Q}}_{\perp}^{\Pi}$
$\displaystyle=i\phi_{k}\left[(\bm{\mathrm{\Pi}}_{k}\cdot\bm{k})\times\hat{\bm{z}}\right]=i\phi\int\frac{1}{2}\left(\begin{array}[]{c}{k_{y}}v_{\perp}^{2}-k_{y}v_{\perp}^{2}\cos(2\alpha)+2{k_{z}}v_{\perp}v_{z}\sin\alpha+k_{x}v_{\perp}^{2}\sin(2\alpha)\\\
-k_{x}v_{\perp}^{2}-k_{x}v_{\perp}^{2}\cos(2\alpha)-{k_{y}}v_{\perp}^{2}\sin(2\alpha)-2{k_{z}}v_{z}v_{\perp}\sin\alpha\\\
0\end{array}\right)f\mathrm{d}^{3}v.$ (33)

The energy flux is assembled using the response derived in equation (25). The
velocity integrals have been computed numerically for the given parameters.

###  2 Verification of energy fluxes

In the numerical testing, the contribution to the energy flux derived in
section 1 are computed from the numerically determined particle flux
$\bm{\Gamma}=\int\bm{v}f\mathrm{d}^{3}v$, the energy density
$\epsilon=\int\frac{v^{2}}{2}f\mathrm{d}^{3}v$, and the components of the
stress tensor $\mathbf{\Pi}=\int\bm{v}\bm{v}f\mathrm{d}^{3}v$. Subsequently,
these numerical results are compared to the analytical descriptions for the
individual contributions to the total energy flux as derived in equations
(23). The analytical calculations involve the introduction of the particle
flux (28), energy density (1), and the stress tensor response, as described in
equation (23).

Figure 2: Comparison (left column) of energy fluxes (Poynting flux $\bm{S}$,
$\bm{E}\times\bm{B}$ heat flux $\bm{Q}^{\bm{E}\times\bm{B}}$, and stress
induced heat flux $\bm{Q}^{\Pi}$) opposite to the temperature gradient for
various wave numbers between (1) results from numerical simulations
($\bullet$) and (2) analytical computations (-). Respective deviation from the
analytical results are shown in the second column

In the simulations detailed in section 2, where the growth rate and frequency
of various modes have been determined, a linearized model is used. In the
linearized limit, the velocity advection reduces to
$\nabla\phi\cdot\nabla_{v}f\approx-\nabla\phi\cdot\bm{v}f_{M}$ to remove
nonlinear effects. However, while this simplified model accurately reproduces
the particle flux, higher moments of the distribution function are not
faithfully represented. For these simulations, the full advection of the
background distribution by the electric field must be considered. To prevent
saturation, a very small initial perturbation amplitude (here, $\delta n\sim
10^{-8}$) must be chosen. All other parameters remain consistent with the
simulation of the dispersion relation in Section 2 ($N=128\times 128\times
8\times 32\times 32\times 32$, $L_{x}=L_{y}=\frac{20}{3}\pi$, $L_{z}=480\pi$,
$\Delta t=0.03$), and the distribution function is initialized with a white
noise density perturbation.
The results for the contributions of the energy flux are shown in figure 2 in
comparison to the analytically obtained results together with the absolute
error in the left column. The $\bm{E}\times\bm{B}$ heat flux is an order of
magnitude larger than the other two fluxes. For large wave numbers, the
$\bm{E}\times\bm{B}$ heat flux is proportional to the $y$-component of the
wave vector $k_{y}$. The values of the Poynting flux and the stress induced
heat flux are nearly identical, as anticipated, given that the gradient has
been selected close to the gyrokinetic limit. In gyrokinetic theory, the
entire energy flux is determined by the $\bm{E}\times\bm{B}$ heat flux.
The code results demonstrate strong agreement with analytical predictions.
Absolute errors are less than 10% for both the pointing flux and the stress-
induced heat flux while the relative error of the $\bm{E}\times\bm{B}$ heat
flux is approximately 1%.

##  V Non-linear ITG simulations

###  1 Comparison with gyrokinetic simulation

Having demonstrated the effective performance of the simulations in the linear
regime, we have delved into the investigation of the nonlinear saturation of
the ITG instability. For this purpose, a simulation similar to section 2 has
been conducted. However, the grid size has been adjusted from $N=128\times
128\times 8$ to $N=64\times 64\times 64$, specifically targeting the necessity
for a higher resolution in $z$-direction due to the emergence of high-
frequency waves during the nonlinear phase. This adjustment ensures that the
various $x-y-$slices remained connected.
Figure 3 displays snapshots of the density perturbation in the $x-y$ plane
during three distinct phases of the simulation. In the first plot, a distinct
mode structure is evident, characteristic of a slab ITG instability where the
system is dominated by modes with $\bm{k}\perp\nabla T$ (in this case, $\nabla
T/|\nabla T|=\hat{\bm{x}}$). The growth rates decrease rapidly as the wave
number parallel to the temperature gradient increases. The second snapshot
depicts the density during the transition from the linear to the nonlinear
phase. The third image represents the perturbation in a fully saturated
turbulent state.
The nonlinear saturation of the simulation is compared to the gyrokinetic code
CGYRO candy_high-accuracy_2016 . To achieve this, a simulation with the same
spatial parameters (length and resolution) is conducted. In gyrokinetic codes,
the velocity space is not fully resolved as in 6D kinetic simulations but is
parameterized by two coordinates: the kinetic energy,
$\epsilon=[0,\epsilon_{\text{max}}]$, and the cosine of the pitch angle,
$\xi\in[-1,1]$. The pitch angle characterizes the ratio between parallel and
perpendicular velocity. The chosen parameters for the simulation are
$\epsilon_{\text{max}}=8$, $N_{\epsilon}=12$, and $N_{\xi}=32$. This
simulation is initiated with a random noise density perturbation having the
same amplitude as the kinetic simulation.

Figure 3:  Snapshots of the particle density in the $x$-$y$-plane (with fixed
$z$) for various points in time

For this comparison, emphasis will be placed on the heat flux induced during
both the linear and nonlinear phases. Section 1 has established that the sole
heat flux present in gyrokinetic theory is the $\bm{E}\times\bm{B}$ heat flux.
Therefore, this serves as the metric for comparison. In Figure 4 (top), the
heat flux parallel to the gradient is depicted for the two simulations.
Following a linear phase (up to $t=3000$), during which the flux experiences
exponential growth, the energy flux reaches saturation. The saturation levels
are similar for both simulations; however, in both cases, the levels are
susceptible to significant fluctuations. This is partially attributed to the
small domain size. The smallest wave number in the system is $k_{0}=0.3$,
whereas the dominant growing mode has a wave number $k_{y}=0.9$. Consequently,
only three wavelengths of the dominant mode fit into the system. Increasing
the domain size would lead to a stronger averaging and thus, a more stable
saturation level. In addition to expanding the domain, a longer runtime would
enhance the comparison by providing more data and enabling a statistical
analysis. However, due to the substantial computational cost associated with
the 6D kinetic simulation, this option has been foregone, and the obtained
results are considered satisfactory.

For further verification, we conducted a comparison of the ratio between the
heat flux and the $L_{2}$-norm of the electrostatic potential, as depicted in
Figure 4 (bottom). For the dominant mode, with
$\bm{k}=(\frac{9}{10},0,\frac{1}{240})$, a ratio of approximately 0.9 is
anticipated, as indicated in Figure 2 (2nd row, left), which is in close
agreement with the simulated results. Following the nonlinear phase, the ratio
decreases to approximately 0.25, a consistent result for both models. In
summary, a very good agreement is observed between the two simulations.

A notable distinction between the two simulations is evident. When examining
the system in the turbulent state, one observes the presence of high-frequency
oscillations with a frequency closely aligned with the Larmor frequency. The
existence of such high-frequency oscillations in the saturation of nonlinear
ion temperature gradient simulations has been previously acknowledged in the
literature miecnikowski_nonlinear_2018 . However, a comprehensive description,
let alone an investigation into their excitation, is still lacking. Figure 5
displays a spectrogram of the electrostatic potential. The ITG intensity is
visible at the bottom of the figures with frequencies close to zero. Modes
with frequencies close to the harmonics of the Larmor frequency are excited
during the nonlinear saturation phase (at $t\sim 3500$). The modes are
clustered in frequency bands close to the harmonics with a mean frequency
slightly larger than the harmonic. The distribution of the high frequency
waves suggest the presence of ion Bernstein waves (IBWs) candy_high-
accuracy_2016 . More precisely neutralized IBWs are electrostatic waves in the
ions with the presence of adiabatic electrons. Various mechanisms have been
proposed to explain the instability of Ion Bernstein Waves (IBWs)
raeth_high_2023 ; yoon_bernstein_2014 ; noreen_ion_2019 . Our analysis
indicates that the excitation is triggered by local negative velocity space
gradients of the distribution function, induced by the temperature gradient
source term rath_beyond_2023 . These findings underscore the necessity of a
nonlinear treatment of the temperature gradient in the full-f 6D model.

Figure 4: Comparison between BSL6D and gyrokinetic simulation (CGYRO) for the
heat flux (a) and the ratio from heat flux and electrostatic potential  Figure
5: Spectrogram of electrostatic potential amplitude $\phi(\omega,t)$

###  2 Nonlinear treatment of gradients

This section transitions from a treatment of the gradients in a local limit,
using the Boussinesq approximation to a nonlinear approach for handling
gradients in our simulation code. The chosen initial condition in the
simulation is designed to allow the distribution function to exhibit a density
and temperature profile

$\displaystyle f_{0}(\bm{r},\bm{v})=\frac{n(\bm{R})}{(2\pi
T(\bm{R}))^{\frac{3}{2}}}e^{-\frac{v^{2}}{2T(\bm{R})}}.$ (34)

The background profiles $n(\bm{r}-\bm{\rho})$ and $T(\bm{r}-\bm{\rho})$ are
defined in gyrocenter coordinates $\bm{R}:=\bm{r}-\bm{\rho}$ (where
$\bm{\rho}:=\hat{\bm{z}}\times\bm{v}$ is the Larmor radius vector) to
establish a background that does not oscillate with the Larmor frequency. To
simplify the treatment of boundary conditions, the profiles are periodically
set up using a sine-profile in the $x$-direction

$\displaystyle n(\bm{R})$
$\displaystyle=1+\kappa_{n}\sin(k_{0}\left(x-v_{y}\right)),$ (35)
$\displaystyle T(\bm{R})$
$\displaystyle=1+\kappa_{T}\sin(k_{0}\left(x-v_{y}\right)).$ (36)

To prevent the background density gradient from generating an electric field,
all modes with wave numbers parallel to the gradient are removed from the
electrostatic potential by subtracting the flux surface average over the
$y-z$-plane

$\displaystyle\phi=n-\langle n\rangle_{y,z}=n-\frac{1}{L_{y}L_{z}}\int
n\mathrm{d}y\mathrm{d}z.$ (37)

For the simulation, the parameters $\kappa_{n}=0$, $\kappa_{T}=0.5$, and
$k_{0}=0.2$ have been chosen, resulting in a temperature gradient
$\max_{x\in[0,L_{x}]}\frac{\partial_{x}T(x)}{T(x)}=0.115$, where $L_{x}$
represents the box length in the $x$-direction. The simulation has been
conducted on a box with a length of $L=10\pi\times\frac{5}{2}\pi\times
240\pi$, $N=128\times 32\times 16\times 32\times 32\times 16$ and $\delta
t=0.02$).

Figure 6: Snapshots of electrostatic potential in simulation of an unstable
ITG modes with nonlinear treatment of temperature gradient at two different
times during the linear (up) and nonlinear (down) phase

A snapshot of the perturbation in the linear phase is displayed in figure 6.
The perturbation exhibits two peaks situated at the maxima of the normalized
gradient $\frac{\partial xT(x)}{T(x)}$. To determine the wavenumber, a sine
wave is fitted to the envelope of the unstable within the full width half
maximum (FWHM) of the perturbation (in x-direction), resulting in a wavenumber
of approximately $k_{x}\approx 0.44$. To facilitate a comparison between the
nonlinear gradient and the analytical calculation, the mean of the gradient is
computed across the mode profile. Upon determining the wavenumber and
effective gradient, the solution of the dispersion relation is obtained from
the analytical expression in equation (13). The anticipated complex frequency
for the fastest-growing mode is $\omega=0.01850+0.00420i$.

Figure 7: Root mean square of electrostatic potential compared to expected
growth rate

Figure 7 displays the root-mean-square (RMS) of the electrostatic potential on
a logarithmic scale (shown in blue), alongside the expected linear growth rate
represented by a dashed gray line. The results indicate that the analytically
predicted growth rate is slightly higher than the expected value. The growth
rate and frequency for the simulation are derived from the Fourier-transformed
electrostatic potential $\phi(\bm{k},t)$, resulting in

$\displaystyle\omega_{\text{BSL}}=0.01809+0.00409i.$ (38)

The growth rate and frequency exhibit an error of approximately 2% when
compared to the analytically calculated values.

After the linear phase, the simulation gradually reaches saturation around
$t\approx 6000$ and settles into a fully saturated turbulent state (compare
figure 6(below)). This simulation marks the first of its kind in modeling
developed ITG turbulence within a 6D kinetic model with a nonlinear treatment
of the gradients.

##  VI Summary

Six-dimensional kinetic simulations, which enable the representation of modes
with frequencies around the Larmor frequency (such as ion Bernstein waves),
serve as a valuable tool for investigating physics beyond conventional
gyrokinetic models. The development of our semi-Lagrangian solver for the
Vlasov system marks the initial stride in constructing a comprehensive kinetic
simulation code for plasma simulations at all tokamak relevant parameters.

The verification tests conducted demonstrated excellent agreement in the
growth rate and frequency for the gyrokinetic modes induced by the Ion
Temperature Gradient (ITG) instability, with the values computed from the
analytical dispersion relation (13). This consistency extends across a broad
range of wave numbers. Additionally, we established that the code accurately
reproduces the quasi-linear fluxes throughout the domain. In the course of
this, a novel formulation of the energy flux in the 6D kinetic system has been
developed, representing a sum of its distinct components (namely, Poynting
flux, $\bm{E}\times\bm{B}$ energy flux, and stress-induced energy flux). The
explicit expression for the energy flux establishes a connection to
gyrokinetic theory, as the $\bm{E}\times\bm{B}$ heat flux can be recognized as
one of its components. Additionally, the individual contributions can be
analytically computed, facilitating a comparison with the simulation results.
Beyond the linear verification of our model, we have conducted simulations
extending well into the nonlinear phase. Comparisons with the nonlinear
gyrokinetic code CGYRO indicate that the code accurately depicts the nonlinear
saturation of a slab ITG instability, yielding consistent saturation
amplitudes. Examination of the saturation phase and the ensuing turbulent
state has led to the identification of the excitation of high-frequency modes.
The gyrokinetic model, by definition, does not account for these excited ion
Bernstein waves. The potential impact of these modes on the saturation process
and energy transport remains unexplored. Some studies from the 1990s suggested
that these waves could suppress turbulence levels in the plasma edge
craddock_theory_1994 ; crawford_absolute_1965 . Additionally, alongside
simulations employing a Boussinesq limit of the gradients, we have
demonstrated the reproducibility of ITG simulations with a nonlinear treatment
of the gradients. The possession a tool that enables the simulation of such
modes presents a valuable opportunity to uncover new physics in regimes beyond
commonly utilized models.

###### Acknowledgements.

This work has been carried out partly within the framework of the EUROfusion
Consortium, funded by the European Union via the Euratom Research and Training
Programme (Grant Agreement No 101052200 – EUROfusion). Support has also been
received by the EUROfusion High Performance Computer (Marconi-Fusion). Views
and opinions expressed are however those of the author(s) only and do not
necessarily reflect those of the European Union or the European Commission.
Neither the European Union nor the European Commission can be held responsible
for them. Numerical simulations were performed at the MARCONI-Fusion
supercomputer at CINECA, Italy and at the HPC system at xthe Max Planck
Computing and Data Facility (MPCDF), Germany.

## References

  * [1] M. Brambilla.  Kinetic theory of plasma waves: homogeneous plasmas.  Number 96 in The International series of monographs on physics. Clarendon Press, 1998.
  * [2] J. Candy, E. Belli, and R. Bravenec.  A high-accuracy Eulerian gyrokinetic solver for collisional plasmas.  Journal of Computational Physics, 324:73–93, 2016.
  * [3] G. G. Craddock, P. H. Diamond, M. Ono, and H. Biglari.  Theory of ion bernstein wave induced shear suppression of turbulence.  Physics of Plasmas, 1(6):1944–1952, 1994.
  * [4] F. W. Crawford and J. A. Tataronis.  Absolute instabilities of perpendicularly propagating cyclotron harmonic plasma waves.  Journal of Applied Physics, 36(9):2930–2934, 1965.
  * [5] G. V. Gordeev.  Soviet physics, JETP, 6(660):660, 1952.
  * [6] F. Jenko.  Massively parallel Vlasov simulation of electromagnetic drift-wave turbulence.  Computer Physics Communications, 125(1):196–209, 2000.
  * [7] S. Jolliet, A. Bottino, P. Angelino, R. Hatzky, T. Tran, B. Mcmillan, O. Sauter, K. Appert, Y. Idomura, and L. Villard.  A global collisionless PIC code in magnetic coordinates.  Computer Physics Communications, 177(5):409–425, 2007.
  * [8] K. Kormann, K. Reuter, and M. Rampp.  A massively parallel semi-Lagrangian solver for the six-dimensional Vlasov–Poisson equation.  The International Journal of High Performance Computing Applications, 33(5):924–947, 2019.  _eprint: https://doi.org/10.1177/1094342019834644.
  * [9] M. Kotschenreuther, G. Rewoldt, and W. Tang.  Comparison of initial value and eigenvalue codes for kinetic toroidal plasma instabilities.  Computer Physics Communications, 88(2):128–140, 1995.
  * [10] R. G. Littlejohn.  Hamiltonian formulation of guiding center motion.  The Physics of Fluids, 24(9):1730–1749, 1981.
  * [11] M. T. Miecnikowski, B. J. Sturdevant, Y. Chen, and S. E. Parker.  Nonlinear saturation of the slab ITG instability and zonal flow generation with fully kinetic ions.  Physics of Plasmas, 25(5):055901, 2018.
  * [12] N. Noreen, F. Riaz, S. Malik, and S. Zaheer.  Ion bernstein mode instability with ring velocity distribution function.  Progress of Theoretical and Experimental Physics, 2019(5):053E01, 2019.
  * [13] M. Raeth and K. Hallatschek.  High frequency non-gyrokinetic turbulence at tokamak edge parameters.  Physical Review Letters, [submitted], 2023.
  * [14] M. Räth.  Beyond gyrokinetic theory: Excitation of high-frequency turbulence in 6d Vlasov simulations of magnetized plasmas with steep temperature and density gradients, 2023.
  * [15] E. Sonnendrücker, J. Roche, P. Bertrand, and A. Ghizzo.  The semi-lagrangian method for the numerical resolution of the Vlasov equation.  Journal of Computational Physics, 149(2):201–220, 1999.
  * [16] B. J. Sturdevant, Y. Chen, and S. E. Parker.  Low frequency fully kinetic simulation of the toroidal ion temperature gradient instability.  Physics of Plasmas, 24(8):081207, 2017.
  * [17] B. J. Sturdevant, S. E. Parker, Y. Chen, and B. B. Hause.  An implicit delta f particle-in-cell method with sub-cycling and orbit averaging for Lorentz ions.  Journal of Computational Physics, 316:519–533, 2016.
  * [18] R. E. Waltz and Z. Deng.  Nonlinear theory of drift-cyclotron kinetics and the possible breakdown of gyro-kinetics.  Physics of Plasmas, 20(1):012507, 2013.
  * [19] P. H. Yoon, F. Hadi, and A. Qamar.  Bernstein instability driven by thermal ring distribution.  Physics of Plasmas, 21(7):074502, 2014.


# Surrogate Modeling of the CLIC Final-Focus System using Artificial Neural
Networks

J. Ögren  European Organization for Nuclear Research (CERN), Geneva,
Switzerland  C. Gohil  European Organization for Nuclear Research (CERN),
Geneva, Switzerland  John Adams Institute, University of Oxford, Oxford,
United Kingdom  D. Schulte  European Organization for Nuclear Research (CERN),
Geneva, Switzerland

###### Abstract

Artificial neural networks can be used for creating surrogate models that can
replace computationally expensive simulations. In this paper, a surrogate
model was created for a subset of the Compact Linear Collider (CLIC) final-
focus system. By training on simulation data, we created a model that maps
sextupole offsets to luminosity and beam sizes, thus replacing computationally
intensive tracking and beam-beam simulations. This model was then used for
optimizing the parameters of a random walk procedure for sextupole alignment.

##  1 Introduction

The field of machine learning has grown immensely in the past decade and is
now finding its way into a wide range of fields. Particle accelerators are
controllable physical systems with rich complexities and offer interesting
applications that can profit from different machine learning methods. Examples
of applications in particle accelerators are anomaly detection in
superconducting magnets [1], automatic control of collimators [2], electron
bunch profile prediction [3, 4], electron beam size prediction in circular
accelerators[5] and ion source prediction [6]. Different advanced techniques
have also been used for optimzation and automatic control in free-electron
lasers [7, 8], longitudinal phase space control [9], dynamic aperture
optimization [10] and optics correction [11, 12]. A more comprehensive list of
machine learning applications in particle accelerators can be found in [13].

The field of particle accelerators relies to a large extent, especially in
design and optimization, on high-fidelity simulations and there is an
abundance of simulation tools. In many cases simulations are computationally
expensive, which makes optimization time consuming and online modeling
problematic. A surrogate model is an approximative model that can mimic the
system and is particularly useful for applications that require extensive
simulations. The surrogate model, which is fast to evaluate, can be used for
optimization or as an online model. In this paper we make use of artificial
neural networks and supervised deep learning to create a surrogate model for a
part of the Compact Linear Collider (CLIC) with computationally expensive
simulations: the final-focus system.

##  2 The CLIC Final-Focus System

CLIC is a proposed linear electron–positron collider [14, 15, 16] at CERN. To
produce high luminosities, linear colliders rely on ultra-small beam sizes at
the interaction point. This requires emittance preservation along the whole
machine and puts tight tolerances on alignment and other imperfections. Being
able to reach nominal performance under imperfect conditions is crucial for
the reliability of such a machine. Nominal parameters for the 380 GeV energy
stage are summarized in Table 1.

The final-focus system constitutes the final 780 m of each beamline, see
Figure 1. Two strong quadrupole magnets at the end of the system de-magnify
the beam in both transverse directions. Since particles of different energies
are bent differently in a magnetic field, a beam with an energy spread will
not be focused to a single point. This effect is called chromaticity and
without compensation the beam size at the interaction point will be greatly
inflated with diminished luminosity as a result. The chromaticity correction
involves nonlinear magnets and a delicate compensation of unwanted
aberrations. Complex systems such as this require fine-tuning and are
challenging to design.

The current version of the CLIC final-focus system has the final quadrupole
doublet mounted outside the detector [17, 18] with $L^{*}=6$ m—distance
between the final quadrupole and point of collision. The system consists of 20
quadrupole magnets, 6 sextupole magnets and 2 octupole magnets. Bending
magnets generate dispersion and in combination with the sextupole magnets, an
energy-dependent focusing is generated that compensates the chromaticity
generated in the final doublet. The final-focus system follows a local
chromaticity correction scheme [19] where two sextupoles are placed at the
final doublet and then additional sextupoles are placed upstream with optics
designed in such a way that unwanted geometric aberrations cancel. Figure 1
shows the beta functions ($\beta_{x,y}$) and the dispersion profile
($\eta_{x}$) of the final-focus system.

Table 1: 380 GeV CLIC beam parameters. Parameter | unit | value
---|---|---
Norm. emittance (end of linac) $\gamma\epsilon_{x}/\gamma\epsilon_{y}$ | [nm] | 900 / 20
Norm. emittance (IP) $\gamma\epsilon_{x}/\gamma\epsilon_{y}$ | [nm] | 950 / 30
Beta function (IP) $\beta_{x}^{*}/\beta_{y}^{*}$ | [mm] | 8.2 / 0.1
Target IP beam size $\sigma_{x}^{*}/\sigma_{y}^{*}$ | [nm] | 149 / 2.9
Bunch length $\sigma_{z}$ | [$\mu$m] | 70
rms energy spread $\delta_{p}$ | [%] | 0.35
Bunch population $N_{e}$ | [$10^{9}$] | 5.2
Number of bunches $n_{\mathrm{b}}$ |  | 352
Repetition rate $f_{\mathrm{rep}}$ | [Hz] | 50
Luminosity $L$ | [$10^{34}\mathrm{cm}^{-2}\mathrm{s}^{-1}$] | 1.5
Peak luminosity $L_{p}$ | [$10^{34}\mathrm{cm}^{-2}\mathrm{s}^{-1}$] | 0.9
Figure 1: Optical functions for the 380 GeV CLIC final-focus system with
$L^{*}$ = 6 m. Six sextupoles are placed in dispersive sections in a local
chromaticity correction scheme.

###  2.1 Tuning

The final-focus system is sensitive to imperfections and the challenge is to
achieve small beam sizes, and high luminosity, for a system with
imperfections. Successful operation of a linear collider requires effective
and fast tuning. Tuning studies have been done extensively for all parts of
CLIC and for the final-focus system in particular. A recent report [20]
considered single-beam tuning of the 380 GeV CLIC final-focus system including
static imperfections such as magnet offsets, roll errors and strength errors.
Single beam means that only one side, i.e. one of the two beamlines, of the
final-focus system is tracked and the beam at the interaction point is
mirrored for the beam–beam simulation.

In that study, the tuning procedure consisted of beam-based alignment that
corrected the linear system, followed by sextupole alignment and tuning with
sextupole and octupole knobs. Due to the nonlinearities of the system together
with synchrotron radiation, small imperfections can have dramatic impact on
beam size at the interaction point and consequently a substantial decrease in
luminosity. The purpose of the tuning simulations is to find robust algorithms
for achieving luminosity under imperfect conditions and thus ensuring the
reliability of the collider. In all the simulations we use PLACET [21] for
tracking the beam through the final-focus system lattice and GUINEA-PIG [22]
for a full beam–beam simulation.

###  2.2 Sextupole misalignment

Let us consider a single imperfection: transverse misalignment of sextupoles.
Experience has shown that this imperfection has a large impact on luminosity.
To quantify the impact we run simulations on a perfect machine and add
transverse sextupole offsets assuming a normal distribution. For different rms
values, we track 20 cases and run a beam–beam simulation to compute the
luminosity. Similar to Ref. [20] we consider the single-beam case. Figure 2
shows the average luminosity decreasing as the rms value of the transverse
sextupole offset increases. Since transverse sextupole misalignments have a
large impact, it is essential to have a quick and robust method for tuning the
sextupoles. In the next section we build a surrogate model for the impact of
sextupole misalignments on luminosity and beam size at the interaction point.

Figure 2: The average luminosity for the perfect machine with sextupole
offsets of varying rms values. As the sextupole misalignments grow, luminosity
quickly deteriorates. The shaded area show the maximum and minimum values.

##  3 Surrogate Model

We only consider transverse sextupole offsets and aim to develop a model that
outputs: luminosity, peak luminosity, horizontal beam size and vertical beam
size. The peak luminosity is the luminosity of collisions within 1% of the
nominal energy. CLIC specifications require a minimum of 60% of collisions
occur at the nominal energy, see Table 1.

An analytical approach is not feasible due to the level of complexity of the
system. The nonlinear fields of the sextupole magnets together with effects of
synchrotron radiation make it difficult to compute the resulting beam
distribution at the interaction point. Furthermore, the disruption parameter
in the beam–beam interaction is high, which means that the beams have
considerable effect on each other. In an electron–positron collider the
transverse electromagnetic fields of one beam focuses the opposing beam
(‘pinch effect’). In such circumstances analytical estimates of luminosity are
not reliable.

Artificial neural networks are large connections of small units (artificial
neurons) and by training on data the models can learn to perform different
tasks such as classification, clustering or regression. In our case we make
use of regression and train a model to map input to known outputs. In the
context of machine learning, this is known as supervised learning. We use the
machine learning system TensorFlow [23] interfaced via the Python library
Keras [24] for training the models. The goal of a model is to make accurate
predictions based on unseen data. Such models are said to generalize well.
Large data sets are often needed to train artificial neural networks
successfully.

###  3.1 Data generation

Assuming normal distributions with rms values of 5, 10 and 20 $\mu$m, we added
random offsets to the sextupoles. For each case we tracked $10^{5}$
macroparticles through the final-focus system, ran a beam-beam simulation and
saved sextupoles positions, beam sizes at the interaction point and the output
from the beam-beam simulation: luminosity and peak luminosity. The data was
generated in parallel using the CERN batch system and in a few days we had
accumulated a data set of 450,000 cases. Each case took approximately 5
minutes to run on an ordinary desktop computer (8 cores, CPU computations
only). Figure 3 shows histograms of the accumulated data. We use 100 bins in
the histogram because we found that this was sufficient to resolve the shape
of the distribution.

Since the numeric values of the inputs and outputs have very different ranges
we must transform data before training our models. We standardize the data
set, i.e. apply transformations $X_{standard}=(X-\mu)/\sigma$, so that each
column in our data set has zero mean and unit variance. The mean $\mu$ and
standard deviation $\sigma$ are computed using the training data set. To the
columns containing output data (luminosity, peak luminosity, horizontal and
vertical beam sizes) we also apply a power transformation to make
distributions closer to Gaussian, see Fig. 4. Standardization and power
transformation are common practices to make the training of artificial neural
networks numerically stable. Both these transformations are available in the
Python library Scikit-learn [25]. When testing the model, the inverse
transformations are applied to the output predictions.

Figure 3: Histograms of the data set showing the four outputs: luminosity (top
left), peak luminosity (top right), horizontal beam size (bottom left) and
vertical beam size (bottom right). The vertical line in the luminosity plot
shows a cut that was used before training the models. Figure 4: Histograms of
the data set after power transformation and standardization: luminosity (top
left), peak luminosity (top right), horizontal beam size (bottom left) and
vertical beam size (bottom right). The data is also standardized to have zero
mean and unit variance. The vertical line in the luminosity plot shows a cut
that was used before training the models.

###  3.2 Model setup

There are no clear guidelines of how to set up an artificial neural network
and how to select a model architecture since it depends on the problem. In
general, if the model has a complexity that is too low it will not be able to
mimic nonlinear behavior and if the model has too much complexity it can
suffer from overfitting, where the model does not generalize well.

First we consider a simple case were our model consists of two hidden layers
of twelve neurons each. We use a model with fully connected layers in a
feedforward configuration (every neuron in one layer connects to every neuron
in the subsequent layer). Figure 5 shows a schematic of the model. The
12-dimensional input vector can be written as

$X=(x_{1},\dots,x_{12})^{T}=(S1_{x},S1_{y},\dots,S6_{x},S6_{y})^{T}$ (1)

where $S1_{x}$ denotes the horizontal position of the first sextupole and so
on. The 4-dimensional output can be written

$Y=(y_{1},\dots,y_{4})^{T}=(L,L_{p},\sigma_{x},\sigma_{y})^{T}$ (2)

where $L$ denotes the luminosity, $L_{p}$ the peak luminosity and $\sigma_{x}$
($\sigma_{y}$) the horizontal (vertical) beam size at IP. The input to a
neuron is the sum of weighted inputs and a bias, and during the training of
the model, weights and biases of all neurons are adjusted to minimize a loss
function—in our case the mean squared error.

Figure 5: A schematic of an artificial neural network model that maps the
horizontal and vertical positions of the six sextupole magnets in the CLIC
final-focus system to the resulting luminosity, peak luminosity and horizontal
and vertical beam sizes. This model consists of two fully connected hidden
layers in a feedforward configuration.

###  3.3 Model Training

Following common practice, we split the complete data set and reserve 20% for
testing. The remaining 80% is split into a training set (90%) and a validation
set (10%). During training the model trains on the training set and evaluates
on the validation set. The purpose of validation is to monitor overfitting.
After training is done, i.e. when there is either no improvements or
overfitting starts, the model performance is evaluated on the testing
data—data never previously seen by the model.

To assess the model performance we compare the predictions to the known
outputs in the test data set. As a figure of merit we use the mean absolute
percentage error (MAPE). From Figure 3 it is clear there are a few outliers
with low luminosity. This results in large errors which can inflate the MAPE.
Therefore, we make a cut in our testing data and only consider cases with
luminosities above $5\times 10^{31}$ cm-2s-1.

Figure 6 shows histogram of the performance of the three different models. We
tested a different number of hidden layers (2-9) and the plots show three
cases: two, five and eight layers. It is clear that a model with two layers
does not have sufficient complexity and gives rather large relative errors.
There is a substantial improvement when we increase the number of layers from
two to five and eight. Table 2 shows the mean of the absolute percentage
error. Eight layers gave the lowest MAPE and there was no improvement from
increasing the number of layers beyond eight. We also tried other network
architectures, for example: we increased the number of nodes per layer; tried
reducing the number of nodes in subsequent layers and tried architectures with
a large number of nodes in the first layer followed by smaller layers. In the
end it did not seem possible to lower the MAPE further.

Figure 6: Histograms of the absolute percentage error for the four different outputs: luminosity (top left), peak luminosity (top right), horizontal beam size (bottom left) and vertical beam size (bottom right). Three different models with two, five and eight layers are compared. Table 2: Model evaluation on testing data. Mean absolute percentage error for the four model outputs for three different model architectures. | 2 layers | 5 layers | 8 layers
---|---|---|---
Luminosity [%] | 8.87 | 4.01 | 2.47
Peak luminosity [%] | 8.42 | 3.77 | 2.29
Horizontal beam size [%] | 2.06 | 0.913 | 1.02
Vertical beam size [%] | 4.74 | 1.68 | 1.04

###  3.4 How much data is needed?

The full data set used in the previous section consisted of 450,000 samples
containing misalignments drawn from a random distribution. 360,000 samples
were randomly selected for the training dataset and the remaining 90,000
samples were used for testing. To investigate the influence of number of data
points, we train the model with eight layers on training sets of different
sizes. Subsets of [1000, 5000, 50000, 100000, …, 300000] are drawn randomly
from the full training set. Each case is trained ten times and tested on the
same full testing set used in previous section. Figure 7 shows the mean and
standard deviation of the MAPE for the models trained on different number of
data points. It is clear that more data is better but with diminishing
returns. Already at 50,000 data points the model performance approaches the
final result. For a model valid for luminosities above $5\times 10^{31}$
cm-2s-1, 50,000 data points seems sufficient. On the other hand, in order to
increase the range of the model more data and more data over a wider range is
likely required.

Figure 7: Mean absolute percentage error (MAPE) for models trained on
different number of training samples. Each case is trained ten times and the
markers show the mean values and the shaded areas the standard deviation. It
is clear that more data is better but with diminishing returns.

###  3.5 Final Model

The final model consisted of eight fully connected layers in a feedforward
configuration, c.f. Figure 5. Each layer had 12 neurons that used the
rectified linear unit as an activation function. The model was trained on the
full training set of 360,000 data points using the Adam optimizer [23]—an
adaptive learning method implemented in TensorFlow—and mean squared error as
loss function. We used a batch size (number of training samples used in every
iteration) of 50 and trained for 400 epochs (the number of times the complete
training data set is passed through). The final model took a total of 8-10
hours to train on an ordinary desktop computer (8 cores, CPU computation). The
results from testing are presented in Table 2.

Surely there is potential room for improvement. For instance, the
hyperparameter optimization was not done exhaustively and other model
architectures might be better. However, finding the optimum model is not
within the scope of this paper. Instead, we want to show that the model was
good enough to be useful. To assess the final model’s usability we compared it
to full simulations.

##  4 Model Verification

To verify the model performance we compare it to the full simulation. A nice
feature with the tracking code PLACET is it can easily interface Python. This
made it possible to seamlessly load and interface the trained machine learning
model in our simulations. At each iteration we ran the full beam–beam
simulation and at the same time checked the prediction of the machine learning
model for comparison.

Figure 8 shows a random walk tuning of the transverse sextupole position.
Random sextupole offsets were added to a perfect machine and using a random
walk the luminosity was optimized within 350 iterations. At each iteration the
beam distribution was tracked and the luminosity was computed with a full
beam–beam simulation. At each step we also checked the prediction of the
machine learning model (plotted in solid and dotted lines). We plot luminosity
normalized to nominal luminosity, $L_{0}$, which is the luminosity of the
unperturbed lattice. Note that sometimes the luminosity in the random walk
tuning is higher than the nominal value since sextupole offsets can introduce
waist shifts that enhances luminosity via the pinch effect and the nominal
luminosity does not include this effect. Over the range of almost two orders
of magnitude there was excellent agreement between the machine learning model
and the full simulation. The MAPE in this case was 2.6% for luminosity, 3.0%
for peak luminosity, 1.8% for horizontal beam size and 2.6% for vertical beam
size. The MAPE was higher in this simulation compared to the testing dataset,
which showed a MAPE of approximately 2% for luminosity and 1% for beam size.
This is due to the small sample size in this simulation which leads to a
higher variance in the MAPE calculated. The random walk method always used the
luminosity from the full simulation and the machine learning model was only
evaluated parasitically.

The converse was also tested, where the output of the machine learning model
guided the random walk optimization and the full simulation was done
parasitically for comparison. Both cases showed similar results.

Figure 8: Example of a random walk optimization on the full simulation.

##  5 Random Walk Optimization

Running a single case on the full simulation, i.e. tracking $10^{5}$
macroparticles and running the full beam–beam simulation, on a normal desktop
computer takes a few minutes. Making a single prediction with the machine
learning model takes a few milliseconds. Such a speed-up in performance allows
for systematic studies of optimization otherwise not feasible. As an example,
we use our machine learning model to optimize a simple random walk procedure.

For the sextupoles, there are twelve degrees of freedom: six sextupoles and
two transverse directions. We devise a random walk procedure in the following
way. Out of the twelve degrees of freedom a subset of the degrees of freedom
is selected randomly. In this subset a direction is selected at random, i.e.
we specify a random linear combination of the degrees of freedom in the
subset. Four points along the selected direction, $g[-1,-0.5,0.5,1]$ with gain
$g$, are checked. If the highest luminosity of the four points exceeds the
luminosity from the previous iteration, that points is selected. Otherwise the
procedure moves to the next iteration and selects a new subset and direction.
The whole procedure is iterated until a certain luminosity is reached. There
are two parameters to be optimized: gain and size of the subset.

We performed a grid search of gains $[0.5,1.0,1.5,2.0,2.5,3.0]$ and subsets
1-12 on the eight layer artificial neural network model. For each setting we
tested 100 different initial conditions with random sextupole offsets. For
each initial condition we ran the random walk procedure 100 times and recorded
the required number of iterations to reach the target luminosity. Figure 9
shows the average number of iterations needed for the different parameter
settings. For a smaller subset a higher gain decreases the number of
iterations whereas for a larger subset a smaller gain is preferable. The
optimum seems to be a small subset of 2 or 3 and a higher gain.

Figure 9: Grid search of gain and subset used in a random walk procedure. The
points show the average values and shaded area the standard deviation. For
each setting a 100 random cases were tested 100 times.

To verify the optimum parameters from the grid search, we tested the random
walk procedure using a full simulation. Figure 10 shows three different cases:
subsets 2, 6 and 12. For each case, three different gains were tested. The
results are in good agreement with the grid search using the machine learning
model. Using a higher gain is better when the subset is small whereas a high
gain and a large subset lead to large fluctuations and requires more
iterations to reach target.

Figure 10: Random walk procedure tested on the full simulation. Three
different gains were tested using subset = 2 (top), subset = 6 (middled) and
subset = 12 (bottom). A higher gain is better when the subset is small and a
smaller gain when subset is large.

This random walk procedure was rather simple and could surely be optimized
further. In a real machine there are other constraints to consider as well and
in most cases a slow smoothly varying optimization might be preferred over a
quicker procedure with large fluctuations. The point of this exercise was to
show the usability of the surrogate model and that results from optimizations
on the surrogate model translated to the full simulation.

##  6 Conclusions

In this paper we presented a surrogate model that maps sextupole offsets to
resulting luminosities and beam sizes in the CLIC final-focus system. We
trained an artificial neural network on a large set of simulation data. The
model was verified by comparing its predictions to the full simulation with a
mean absolute percentage error of a few percent. We used the model to optimize
two parameters of a random walk procedure and compared the results to the full
simulation, with good agreement.

We considered only a small subset of the system (transverse sextupole offsets)
and the single-beam case. A more complete surrogate model could be realized by
adding more imperfections to the system. Naturally, as the degrees of freedom
increase, training the model becomes more challenging and may require a larger
data set but the principle remains the same. Nonetheless, as we showed with
our example of optimizing the random walk procedure: even a model with limited
scope can already be very useful. In our case it helped improving the
sextupole alignment procedure and reduced the overall tuning time.

## References

  * [1] M. Wielgosz, A. Skoczen, and M. Mertik, “Using LSTM recurrent neural networks formonitoring the LHC superconducting magnets,” Nucl. Instrum. Meth.A867, 40–50 (2017).
  * [2] G. Azzopardi, B. Salvachua, G. Valentino, S. Redaelli and A. Muscat, “Operational results on the fully automatic LHC collimator alignment,” Phys. Rev. Accel. Beams 22, 093001 (2019).
  * [3] C. Emma, A. Edelen, M. J. Hogan, B. O’Shea, G. White, and V. Yakimenko, “Machine learning-based longitudinal phase space prediction of particle accelerators,” Phys. Rev. Accel. Beams 21, 112802 (2018).
  * [4] A. Scheinker and S. Gessner. ”Adaptive method for electron bunch profile prediction.” Physical Review Special Topics-Accelerators and Beams 18.10, 102801 (2015).
  * [5] S. C. Leemann, S. Liu, A. Hexemer, M. A. Marcus, C. N. Melton, H. Nishimura and C. sun. “Demonstration of Machine Learning-Based Model-Independent Stabilization of Source Properties in Synchrotron Light Sources.” Physical review letters, 123, 194801 (2019).
  * [6] Y. B. Kong, M. G. Hur, E. J. Lee, J. H. Park, Y. D. Park, and S. D. Yang. “Predictive ion source control using artificial neural network for RFT-30 cyclotron.” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 806 (2016): 55-60.
  * [7] N. Bruchon, G. F enu, G. Gaio, M. Lonza, F. A. Pellegrino, and L. Saule. “Free-electron laser spectrum evaluation and automatic optimization.” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 871 20-29, (2017).
  * [8] A. Scheinker, D. Bohler, S. Tomin, R. Kammering, I. Zagorodnov, H. Schlarb, M. Scholz, B. Beutner and W. Decking, “Model-independent tuning for maximizing free electron laser pulse energy,” Phys. Rev. Accel. Beams 22, 082802 (2019).
  * [9] A. Scheinker, A. Edelen, D. Bohler, C. Emma and A. Lutman. “Demonstration of model-independent control of the longitudinal phase space of electron beams in the Linac-coherent light source with Femtosecond resolution.” Physical review letters, 121(4), 044801 (2018).
  * [10] Y. Li, W. Cheng, L. H. Yu, and R. Rainer. “Genetic algorithm enhanced by machine learning in dynamic aperture optimization.” Physical Review Accelerators and Beams, 21(5), 054601 (2018).
  * [11] Y. Li, R. Rainer, and W. Cheng. “Bayesian approach for linear optics correction.” Physical Review Accelerators and Beams 22.1, 012804 (2019).
  * [12] E. Fol, J. M. C. Portugal, G. Franchetti, and R. Tomás. “Optics corrections using Machine Learning in the LHC.” In Proceedings of the 2019 International Particle Accelerator Conference, THPRB077, (2019). doi:10.18429/JACoW-IPAC2019-THPRB077
  * [13] A. Edelen et al., “Opportunities in Machine Learning for Particle Accelerators,” arXiv [physics.acc-ph], arXiv:1811.03172, 2018.
  * [14] “A Multi-TeV linear collider based on CLIC technology: CLIC Conceptual Design Report,” edited by M. Aicheler, P. Burrows, M. Draper, T. Garvey, P. Lebrun, K. Peach, N. Phinney, H. Schmickler, D. Schulte and N. Toge, CERN-2012-007.
  * [15] “Updated baseline for a staged Compact Linear Collider,” edited by P.N. Burrows, P. Lebrun, L. Linssen, D. Schulte, E. Sicking, S. Stapnes, M.A. Thomson, CERN–2016–004 (CERN, Geneva, 2016).
  * [16] “The Compact Linear Collider (CLIC) – Project Implementation Plan,” edited by M. Aicheler, P.N. Burrows, N. Catalan, R. Corsini, M. Draper, J. Osborne, D. Schulte, S. Stapnes, M.J. Stuart, CERN-2018-010-M.
  * [17] CLICdp collaboration, “The post-CDR CLIC detector model”, CLICdp-Note-2017-001 (2017).
  * [18] F. Plassard, A. Latina, E. Marin, and R. Tomás and P. Bambade, “Quadrupole-free detector optics design for the Compact Linear Collider final focus system at 3 TeV”, Phys. Rev. Accel. Beams, vol. 21, p. 011002 (2018).
  * [19] P. Raimondi and A. Seryi, “Novel Final Focus Design for Future Linear Colliders”, Phys. Rev. Lett. 86, 3779 (2001).
  * [20] J. Ögren, A. Latina, R. Tomás and D. Schulte, “Tuning of the CLIC 380 GeV Final-Focus System with Static Imperfections”, CERN, Geneva, Switzerland, Tech. Rep. CERN-ACC- 2018-0055, Dec. 2018.
  * [21] The tracking code placet, https://clicsw.web.cern.ch/clicsw/.
  * [22] D. Schulte, “Study of Electromagnetic and Hadronic Background in the Interaction Region of the TESLA Collider”, Ph.D. thesis, DESY-TESLA-97-08, TESLA-97-08, Germany, 1996.
  * [23] Martín Abadi et al., ”TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015\. Software available from tensorflow.org.
  * [24] F. Chollet and others, Keras, 2015. https://keras.io
  * [25] F. Pedregosa et al., ”Scikit-learn: Machine Learning in Python,” JMLR 12, pp. 2825-2830, 2011.