# A Parallelized Vlasov-Fokker-Planck-Solver for Desktop PCs

Patrik Schönfeldt<EMAIL_ADDRESS>Miriam Brosi  Markus Schwarz
Johannes L. Steinmann  Anke-Susanne Müller  Karlsruhe Institute of Technology,
Kaiserstraße 12, 76131 Karlsruhe, Germany

###### Abstract

The numerical solution of the Vlasov-Fokker-Planck equation is a well
established method to simulate the dynamics, including the self-interaction
with its own wake field, of an electron bunch in a storage ring. In this paper
we present _Inovesa_ , a modularly extensible program that uses OpenCL to
massively parallelize the computation. It allows a standard desktop PC to work
with appropriate accuracy and yield reliable results within minutes. We
provide numerical stability-studies over a wide parameter range and compare
our numerical findings to known results. Simulation results for the case of
coherent synchrotron radiation will be compared to measurements that probe the
effects of the micro-bunching instability occurring in the short bunch
operation at ANKA. It will be shown that the impedance model based on the
shielding effect of two parallel plates can not only describe the instability
threshold, but also the presence of multiple regimes that show differences in
the emission of coherent synchrotron radiation.

##  I Introduction

At synchrotron light sources electron bunches of a length of a few millimeters
are used to produce coherent synchrotron radiation (CSR) in the terahertz
(THz) frequency range. Due to the coherent emission, the intensity scales with
the number of emitting particles squared, instead of linearly as for
incoherent emission. In storage rings the spatial compression is achieved by
using magnet optics with a small momentum compaction factor $\alpha_{c}$. The
compression leads to the micro-bunching instability. On the one hand, this
instability limits the electron bunch charge that can be used in stable
operation; on the other hand the emerging substructures emit coherent
radiation also in wavelength smaller than the electron bunch length.

First observations of micro-bunching in a storage ring as well as of the
increase of coherent emission in the spectral range of interest were made at
the NSLS VUV ring Carr  _et al._ (1999, 2001). It then has been studied both
experimentally, to map out the parameters governing the bursting behavior
(e.g. at ANKA Müller (2010), BESSY II Abo-Bakr  _et al._ (2003), DIAMOND
Shields  _et al._ (2012), MLS Wüstefeld  _et al._ (2010), and SOLEIL Evain
_et al._ (2012a)), and theoretically to predict thresholds Venturini and
Warnock (2002); Stupakov and Heifets (2002); Bane  _et al._ (2010); Cai
(2011a) and to simulate the dynamics. To simulate the dynamics, it is possible
to solve the Vlasov-Fokker-Planck equation for the longitudinal phase space
density Venturini  _et al._ (2005); Roussel  _et al._ (2014), or, using
supercomputers, to do particle tracking with one million macro particles or
more Evain  _et al._ (2012b).

Recent advances in detector development and readout electronics facilitate
fast mapping of the micro-bunching instability Brosi  _et al._ (2016) over a
wide range of physical parameters. To cover these settings in simulation as
well, it calls for an ultra-fast simulation technique. Also, the simulation
tool should be designed such that the influences of both the simulated physics
and of numerical effects can be studied and separated. As we are interested in
simulating an instability, in particular sources of numerical instabilities
should be ruled out. In this paper we present _Inovesa_ (_Inovesa Numerical
Optimized Vlasov-Equation Solver Application_ , available at git (2016)), a
Vlasov-Fokker-Planck solver that runs on standard desktop PCs and yields
robust results within minutes.

Section II summarizes the theoretical description of the problem. The actual
implementation is described in section III, while section IV presents
numerical studies that show the robustness of the implementation and compares
results of simulation and measurements.

##  II Longitudinal Phase-Space Dynamics

###  II.1 Vlasov-Fokker-Planck equation

The phenomenon of micro-bunching happens in the longitudinal phase space,
which is spanned by the position $z$ relative to the synchronous particle and
the energy $E$. Taking the particle density $\psi(z,E,t)$ of electrons in a
storage ring to be a smooth function, its evolution with time $t$ can be
described by the Vlasov-Fokker-Planck equation (VFPE). Following the notation
of Bane  _et al._ (2010) it reads

$\frac{\partial\psi}{\partial\theta}+\frac{\partial H}{\partial
p}\frac{\partial\psi}{\partial q}-\frac{\partial H}{\partial
q}\frac{\partial\psi}{\partial p}=\beta\frac{\partial}{\partial
p}\left(p\psi+\frac{\partial\psi}{\partial p}\right),$ (1)

with the time given in multiples of synchrotron periods $\theta=f_{s}t$, the
normalized coordinates $q={z}/{\sigma_{z,0}}$, and
$p={(E-E_{0})}/{\sigma_{\delta,0}}$, the Hamiltonian $H$, the reference
particle’s energy $E_{0}$, and $\beta=1/(f_{s}\tau_{d})$, where $\tau_{d}$ is
the longitudinal damping time. The quantities $\sigma_{\delta,0}$ and
$\sigma_{z,0}$ describe, respectively, energy spread and bunch length in the
equilibrium state that exists for small bunch charges.

To solve this partial differential equation there exists a widely used
formalism by Warnock and Ellison Warnock and Ellison (2000). It uses a grid to
discretize $\psi(q,p)$ and assumes that the collective force due to self-
interaction with the bunch’s own coherent synchrotron radiation is constant
for small time steps. The perturbation due to the collective effects is
described as a perturbation to the Hamiltonian

$\displaystyle H(q,p,t)$
$\displaystyle=\underbrace{H_{e}(q,p,t)}_{\mathrm{external\;fields}}+\underbrace{H_{c}(q,t)}_{\mathrm{collective\;effects}}$
$\displaystyle=\frac{1}{2}\left(q^{2}+p^{2}\right)+Q_{c}\times
V_{c}(Z_{c},q,t),$ (2)

where $Q_{c}$ is the charge involved in the perturbation, and $V_{c}$ is the
potential due to the collective effect, which can be expressed in terms of an
impedance $Z_{c}$.

It is then possible to use the homogeneous solution, which in the unperturbed
case is represented by a rotation in phase space, and add the influence of
diffusion and damping as a particular solution. To model the perturbation, the
influence of $V_{c}$ is implemented as a ‘kick’ along the energy axis.

###  II.2 Micro-Bunching Instability

To calculate the effect of the perturbation term introduced in Eq. 2, one
needs the electric field $E(q,s)$ at the longitudinal position $s$. It is
convenient to express it via a wake potential

$\displaystyle V(q)$ $\displaystyle=\int E(q,s)\,\mathrm{d}s,$ (3)

which directly gives the energy difference for the electrons (in
$\mathrm{eV}$). The wake potential can be obtained from the wake function
$W(q)$, which describes the field produced by one single particle. The wake
potential $V(q)$ then is obtained by convolving it with the charge density
$\varrho$ Warnock and Ellison (2000)

$V(q)=\int_{-\infty}^{\infty}W(q-q^{\prime})\varrho(q^{\prime})\,\mathrm{d}q^{\prime}.$
(4)

As in frequency space closed and smooth expressions exist for many commonly
used impedances, we decided to work in frequency space. Then the wake
potential can be deduced directly from the impedance $Z(k)$ in every time step
using

$V(q)=\int_{-\infty}^{\infty}Z(k)\tilde{\varrho}(k)e^{ikq}\,\mathrm{d}k,$ (5)

where $\tilde{\varrho}(k)$ is the Fourier transform of the bunch profile.

This method allows to implement different impedance models for _Inovesa_ in
just a few lines of code. As we are mostly interested in CSR-driven dynamics,
we currently implemented two cases. The first, the free space CSR impedance,
describes the effect of coherent synchrotron radiation of particles traveling
on a curved path in vacuum. A good approximation for the CSR impedance in a
perfect circle is Murphy  _et al._ (1997)

$Z(n)\approx
Z_{0}\frac{\Gamma(2/3)}{3^{1/3}}\left(\frac{\sqrt{3}}{2}+\frac{i}{2}\right)n^{1/3},$
(6)

where $Z_{0}$ is the vacuum impedance, and $n=f/f_{rev}$ is frequency
expressed in multiples of the revolution frequency $f_{rev}$.

The second implemented impedance approximates the shielding effect of the beam
pipe by two parallel plates with distance $g$ Murphy  _et al._ (1997). It can
be approximated with the Airy functions Ai and Bi Cai (2011b):

$\displaystyle Z(n,R/g)$
$\displaystyle\approx\frac{4\pi^{2}\,2^{1/3}}{\epsilon_{0}c}\left(\frac{R}{g}\right)n^{-1/3}$
$\displaystyle\times\sum\limits_{p}\mathrm{Ai}^{\prime}(u_{p})\mathrm{Ci}^{\prime}(u_{p})+u_{p}\;\mathrm{Ai}(u_{p})\mathrm{Ci}(u_{p}),$
(7)

where the prime marks the first derivative with respect to the argument, $R$
is the beam path’s radius, $\mathrm{Ci}:=\mathrm{Ai}-j\,\mathrm{Bi}$, and

$u_{p}:=\frac{\pi^{2}\left(2p+1\right)^{2}}{2^{2/3}}\left(\frac{R}{g}\right)^{2}n^{-4/3}.$
(8)

Figure 1:  Unshielded (free space) CSR impedance (FS) and CSR impedance
shielded by parallel plates (PP), calculated for the case of an accelerator
with $f_{rev}=8.582\,\mathrm{GHz}$, and (for the shielded case)
$g=32\,\mathrm{mm}$. For high frequencies both impedances converge, as short
wavelength are not effected by the shielding. For $f\rightarrow
0\,\mathrm{Hz}$ both impedances approach $Z=0\,\mathrm{\Omega}$.

This more complex impedance has already proven to describe the micro-bunching
instability threshold within the uncertainty of the measurements Brosi  _et
al._ (2016). For the limit $g\rightarrow\infty$ it converges to the free space
impedance. Figure 1 shows these two impedances, for the case of an accelerator
with $f_{rev}=8.582\,\mathrm{GHz}$, and (for the shielded case)
$g=32\,\mathrm{mm}$.

##  III Implementation

###  III.1 Discretization

Following the approach of Warnock and Ellison Warnock and Ellison (2000), the
VFPE is discretized on a grid to be solved numerically. For _Inovesa_ , we
define a grid starting from a minimum value that can be expressed for each
dimension. For $q>q_{\mathrm{min}}$ and $p>p_{\mathrm{min}}$, the generalized
coordinates $q,p\in\mathbb{R}$ become the grid coordinates
$x_{r},y_{r}\in\mathbb{R}_{\geq 0}$. With $\Delta p$, the granularity of the
grid in energy direction, the transformation between the coordinate systems
can be expressed as

$y_{r}(p)=(p-p_{\mathrm{min}})/\Delta p,$ (9)

and accordingly for $x_{y}(q)$. Function values are only stored at integer
coordinates $m$, where $m$ refers to either $x$ or $y$. When a function value
at an arbitrary non-integer coordinate $m_{r}$ is needed, Inovesa approximates
it by interpolation using

$f(m_{r})\approx\vec{f}_{N}(m)\cdot\vec{p}_{N}(\\{m\\})=:P_{N}(m_{r})$ (10)

were $\\{m\\}=m_{r}-m$ denotes the fractional part of $m_{r}$. The
interpolation multiplies the vector of the function values at the $N$
surrounding mesh points

$\displaystyle\vec{f}_{N}(m)=($ $\displaystyle
f(m-\lfloor(N-1)/2\rfloor),\dots,f(m),$
$\displaystyle\dots,f(m+\lceil(N-1)/2\rceil))^{T}$ (11)

with a vector containing $N$ interpolation coefficients

$\vec{p}_{N}(\\{m\\})=(l_{0,N}(\\{m\\}),\dots,l_{N-1,N}(\\{m\\}))^{T},$ (12)

where $l_{\nu,N}(\\{m\\})$ are the Lagrange basis polynomials Press  _et al._
(2007)

$l_{\nu,N}(\\{m\\}):=\prod\limits_{k=0\neq\nu}^{(N-1)}\frac{\\{m\\}-m_{k}}{m_{\nu}-m_{k}}.$
(13)

We have implemented this for up to $N=4$, which results in a cubic
interpolation scheme. Interpolation using these polynomials sometimes
overshoots the values of the neighboring grid cells. Section IV.2 will show
that these numerical artifacts can even self-amplify and become a very
dominant feature in the simulation. To avoid overshooting, one has to use
clamped or saturated interpolation functions. A simple, clamped version of Eq.
10 reads

$\displaystyle\mathrm{clamp}\left(f(m_{r})\right)=$
$\displaystyle\max\\{\min[\min(f_{N}(m),f_{N}(m+1)),$
$\displaystyle\vec{f}_{N}(m)\cdot\vec{p}_{N}(\\{m\\})],$
$\displaystyle\max(f_{N}(m),f_{N}(m+1))\\}.$ (14)

Keeping in mind that not only $q,p\in\mathbb{R}$ is discretized to
$x,y\in\mathbb{N}$ but also $\psi(x,y)$ is discretized in the computer’s
memory, we analyzed the effect of this second discretization. To be able to
change the accuracy in small steps, we used a fixed point representation fpm
(2009) where the number of bits for the fractional part could be chosen
freely. The result of this test will be shown in the convergence studies
(section IV.1).

###  III.2 Simulation steps

Each time step $f:\psi_{t}(q,p)\rightarrow\psi_{t+\Delta t}(q,p)$ is composed
of a number of simulation steps that model rotation, kick, damping and
diffusion $f=f_{rot}+f_{kick}+f_{damp,diff}$. _Inovesa_ splits the direct
implementation of each of these simulation steps into two sub-steps. The
information on the actual coordinates ($q,p$) is used by the first half
simulation step. We call its result a ‘source map’ (SM). Then, in the second
half step only the grid coordinates $x,y$ are used. Further we define
$z=N_{y}\times x+y$ with the number of grid cells in energy direction $N_{y}$.
Doing so, the function $f:\psi_{t}(\mathbb{R}^{2})\rightarrow\psi_{t+\Delta
t}(\mathbb{R}^{2})$ becomes a one dimensional function
$f_{SM}:\psi_{t}(\mathbb{R}^{1})\rightarrow\psi_{t+\Delta t}(\mathbb{R}^{1})$,
which depends only on $z$.

This method – by construction – produces the same results as the single-step
implementation. Practically speaking, the source map expresses the information
which data of the current simulation step contributes to a grid point for the
next simulation step directly in terms of position in the computer’s memory.
For many functions – such as rotation – the SM will look the same for the
whole runtime of the program. For that reason, it only has to be computed once
during a simulation run and can be kept for multiple usages.

The source map formalism does not only allow to keep intermediate results, it
also gives a handy interface to implement arbitrary functions on the phase
space. Furthermore, the reduction of the problem’s dimension also leads to a
speedup of the calculations. Results of a benchmark of the computational
performance for the particular case of rotation are shown in Fig. 2.

Figure 2: Computational time needed by an Intel Core i5 4258U for 1000
cubically interpolated rotation steps using different implementations. In this
test case, the grid has 512 points per axis, no optimizations (besides SM)
were used. The first bar represents the computational time a direct
implementation of the rotation takes. The second and third bar show the time
the implementation using the source map formalism take – running sequentially
or parallel on the integrated graphics processor.

In this case source mapping halves the runtime. Parallelization using OpenCL
ope (2016) allows further speedup. One advantage of OpenCL is that it can
utilize not only multi-core CPUs but also graphic processors. In total, a non-
optimized program takes days for a typical simulation run. Using the method
described above, _Inovesa_ can reduce this to 15 minutes when running on a
customer grade graphics card.

Using the SM formalism, we implemented different versions of the simulation
steps necessary to solve Eq. 1. For the rotation, we provide a direct
implementation (as in Warnock and Ellison (2000), ‘standard rotation’).
Additionally, we use a symplectic integrator Wolski (2014) to implement the
rotation. This method offers two advantages: First, since the method is
symplectic, it is automatically area preserving – also for truncated power
series. Using the map from an infinitesimal rotation given by a direct
implementation, symplecticity is easily lost in the numeric treatment. Second,
the symplectic method provides additional numerical stability also for the
interpolation. The standard rotation algorithm requires a two-dimensional
interpolation, which involves $N\times N$ data points, leading to
interpolation coefficients
$l_{\nu,N}^{2}\propto\\{x\\}^{N-1}\times\\{y\\}^{N-1}\ll 1$. The symplectic
map, in contrast, splits the rotation into a energy-dependent drift followed
by a location dependent RF kick. Each requires a one-dimensional
interpolation, involving $N$ data points. The resulting coefficients are
$l_{\nu,N}^{1}\propto\\{m\\}^{N-1}\gg l_{\nu,N}^{2}$, minimizing the
vulnerability to numerical absorption. Since a point is rotated by moving on
straight lines in perpendicular directions, in the following the symplectic
approach will be referred to as ‘Manhattan rotation’.

The damping and diffusion terms (right hand side of Eq. 1) need numerical
differentiation. We found that the same type of artifacts that we found for
the interpolation (see section IV.2) can also occur because of the
differentiation. This can be explained by the fact that the algorithm usually
used for numerical differentiation Press  _et al._ (2007) is equivalent to
differentiating the quadratic interpolation polynomial $P_{3}$ (see Eq. 10)

$\frac{\partial f(x)}{\partial
x}\biggr{\rvert}_{x_{0}}\approx\frac{f(x_{0}+\Delta x)-f(x_{0}-1)}{2\Delta
x}=\frac{\partial P_{3}(x)}{\partial x}\biggr{\rvert}_{x_{0}},$

where the distance between the sampling points in our case is $\Delta x=1$. As
a consequence, we target this by using the cubic interpolation polynomial
$P_{4}$, and obtain

$\frac{\partial f(x)}{\partial
x}\biggr{\rvert}_{x_{0}}\approx\frac{-2f(x_{0}-\Delta x)-3f(x_{0})}{6\Delta
x}\\\ +\frac{6f(x_{0}+\Delta x)-f(x_{0}+2\Delta x)}{6\Delta x}=\frac{\partial
P_{4}(x)}{\partial x}\biggr{\rvert}_{x_{0}}.$ (15)

Aside from this improvement, we proceed analogously to Warnock and Ellison
(2000).

To implement the perturbation via a kick, we just had to translate the wake
potential (see Eq. 5) to the grid coordinate system using Eq. 9. Furthermore,
our implementation uses the fact that both $Z(k)$ and $\tilde{\varrho}(k)$ are
Hermitian. This means that optimized algorithms like the ones from Ref. fft
(2016); clf (2016) only need explicit function values for $k\geq 0$ to perform
the inverse Fourier transform. Using this symmetry also brings an improvement
in both speed and memory usage by roughly a factor of two fft (2016).

##  IV Results

###  IV.1 Convergence studies

In this section we compare the effect of different numerical settings, which
ideally should not affect the physical result. We also compare different
implementations of the rotation as described in section III.2. For the sake of
simplicity, we go to the unperturbed case (meaning $H_{c}=0$ in Eq. 2). So any
starting distribution should exponentially converge to a Gaussian distribution
with $\sigma_{q}=\sigma_{p}=1$.

Figure 3:  Exponential damping of the bunch length $\sigma_{z}$ as a function of time for different data types. There is an initial jitter due to quantization noise in the test pattern that has been read in from a 16 bit _PNG_ file and due to the fact that the initial distribution is not fully covered by the grid. Even with this problematic starting conditions, after $T=60\,\mathrm{T_{s}}$, all data types have converged to a constant value. Note that all simulation curves are almost perfectly overlapping, there is no noticeable difference coming from the used data type. The values the simulations converge to are listed in Tab. 1.  Table 1: Deviation of the converged simulation results for the different used data types shown in Fig. 3 from the analytic result ($1\,\mathrm{\sigma_{z,0}}$). Data Type | Deviation / $\sigma_{z,0}$
---|---
bin64 | $0.991\times 10^{-4}$
bin32 | $0.998\times 10^{-4}$
fix59 | $0.991\times 10^{-4}$
fix54 | $0.991\times 10^{-4}$
fix44 | $0.999\times 10^{-4}$
fix34 | $8.877\times 10^{-4}$

At first, we investigate the effect on the results of using different data
types. To do so, we observe the evolution of a Gaussian distributed charge
density with $\sigma_{q}=\sigma_{p}=2.8$, when the damping time is set to five
synchrotron periods $\tau=5\,T_{s}$. We start by reading the distribution from
a 16-bit grayscale PNG. This brings initial quantization noise, but also
provides well defined starting conditions for every data type: Initial
rounding errors will be the same in the different runs. Figure 3 shows the
different simulation results. We find that the results obtained using fixed
point representations with a fractional part of at least 44 bits, and the
tested floating point representations (binary32 and binary64 IEE (2008), often
referred to as ‘single precision’ and ‘double precision’) show a common
difference from the analytic result of about $10^{-4}$. The relative
differences between the data types are less then $10^{-6}$ and therefore can
be neglected. As most libraries are developed focusing on floating point data
types, we implemented the calculation of the wake potential only for those
types. In the following, we default to binary32 because it is much faster than
binary64.

Figure 4:  Exponential damping of the RMS energy spread for different numbers
of computation steps per synchrotron period. (Here the grid size is set to
256.) To emphasize the differences, the analytical, exponential damping
($\sigma_{\delta,t}$) has been subtracted from the simulated values. Every
color represents a different number of simulation steps per synchrotron
period. Until $T=100\,\mathrm{T_{s}}$, all simulations reproduce the
exponential behavior. Manhattan rotation (dashed lines) reproduces the set
values (solid gray line) better than standard rotation (other solid lines)
independently of the number of steps. However, when the number of steps per
synchrotron period becomes too large, the tiny step size leads to systematic
drifts also for this more robust method.

To investigate the influence of the size of the time steps, we observe the
evolution of a Gaussian distributed charge density with
$\sigma_{q}=\sigma_{p}=1.45$. For this we set the damping time to
$\tau=45\,T_{s}$. The Number of simulation steps ($\Delta T$) per $T_{s}$ is
varied between 500 and 4000. All settings reproduce the exponential damping.
However when $\sigma_{q}$ approaches $1$ some of the runs start to diverge. As
depicted in Fig. 4 for Manhattan rotation the error on the reconstructed
values is significantly lower (usually $\ll 1\%$). Furthermore the Manhattan
rotation is more robust against changing the step sizes. Note that the error
increases when the number of steps per synchrotron period becomes too large.
For standard rotation on the other hand there is no obvious optimum for the
size of the time steps.

Figure 5:  Evolution of the RMS energy spread for different grid sizes and
$\Delta t=1/1000\,\mathrm{T_{s}}$. Manhattan rotation (dashed lines) shows
some initial wiggles but reproduces the expected value
($\sigma_{\delta,t}=\sigma_{\delta,0}$, solid grey line) very well. For
standard rotation (other solid lines) on the other hand the divergence is even
worse than shown in Fig. 4.

In another example to test the influence of the grid size we study the
evolution of a Gaussian distributed charge density with
$\sigma_{q}=\sigma_{p}=1$. From the physics point of view it is expected that
the distribution stays constant with time. However, as illustrated in Fig. 5,
it is observed that $\sigma_{p}$ converges to different values or even
diverges depending on the numerical parameters. For the standard rotation with
the specific time step of $\Delta t=1/1000\,\mathrm{T_{s}}$, there is no mesh
size where the simulation converges. Manhattan rotation on the other hand
shows initial jitter with a relative amplitude smaller than $1\%$ and reliably
converges afterwards.

Table 2: Typical relative errors observed in the convergence studies using Manhattan rotation. Parameter | Error
---|---
Switching data types | $<0.1\%$
Changing number of time steps | $<1\%$
Changing number of grid points | $<1\%$

It can be seen that different combinations of the relative time step $\theta$,
the damping factor $\beta$, and the grid size may converge to slightly
different values – or even diverge. We found that Manhattan rotation is much
more robust than the standard approach. Provided that $\theta$ and $\beta$ do
not become too small, the relative error can be reliably kept below $1\%$ (see
Tab. 2).

###  IV.2 Numerical artifacts

Besides the numerical inaccuracies discussed above there are also numerical
artifacts. For those we identify two main sources: interpolation and numerical
differentiation.

As an example, we do two simulation runs with the same current distribution,
one run using cubic interpolation, the second using quadratic interpolation.
For the simulation run that uses cubic interpolation, the distribution stays
in a relatively calm state with just little oscillation. Note that this is not
equilibrium: As expected for the unshielded CSR case and $\xi>0.5$, we have
$\sigma_{\delta}>\sigma_{\delta,0}$ Stupakov and Heifets (2002) and an
oscillation with $f\approx 2f_{s}$. The complete set of simulation parameters
is listed in Table 3.

Table 3: Parameters for an example run to check for numerical artifacts. For the given set no artifacts were observed. Changing to quadratic interpolation however, triggered the occurrence of (non-physical) structures with a period length of two grid cells (see Fig. 7).  Parameter | Value
---|---
Grid points per axis | 256
Steps per $T_{s}$ | 4000
Interpolation method | cubic
Impedance model | free space
Scaled current ($\xi$) | 0.516
Damping time | $200\;T_{s}$

If there was no influence of the different interpolation schemes, one would
expect no difference between the two runs. However, as Fig. 6 depicts, this is
not guaranteed.

Figure 6:  Evolution of the energy spread over time. For the case where cubic
interpolation is used (solid black line), the energy spread stays at
$\sigma_{\delta}\approx 1.01\sigma_{\delta,0}$. When using quadratic
interpolation (dashed blue line), an increase in energy spread can be observed
at $T=1\;T_{s}$. Looking at the bunch profiles at the points in time marked by
the disks (see Fig. 7) reveals that this increase is a numerical artifact.

The energy spread simulated using quadratic interpolation rapidly increases to
a higher value after one synchrotron period ($T=1\;T_{s}$). On a longer time
scale it will damp down again, and after that a new numerical instability
might rise. If the aim is to find out whether the simulated conditions are
above the micro-bunching threshold these artifacts might not matter – below
the threshold the initial, numerical modulation should not be amplified.
However, we want to track the evolution of the charge distribution, and this
numerical artifact might be interpreted as an unphysical slow bursting
frequency. So we have to avoid numerical artifacts as they occur in the run
using quadratic interpolation (depicted in Fig. 7): The period length of the
ripples is exactly two grid cells, and they continue to exist even to a
position where they create negative charge densities. In this particular case,
the instability is clearly triggered by numerical artifacts (overshoots) of
the interpolation.

Figure 7:  Bunch profiles of the simulation run using quadratic interpolation
at selected points in time (cf. discs in Fig. 6). The bunch profiles computed
during the initial increase of the energy spread show large ripples with a
period length of two grid cells. The earlier and later profiles do not show
such structures. This implies that the increase of energy spread is driven by
a numerical instability.

Note that also higher order polynomials show overshoots and that also
numerical differentiation can be expressed in terms of interpolation
polynomials. In our tests, however, we did not find any case where a more
complex differentiation method than the one described in Eq. 15 was needed to
avoid artifacts. In contrast to the differentiation, there were rare cases
where we observed interpolation artifacts even when using cubic interpolation
polynomials. Those we had to suppress by clamping (Eq. 14) – which restricts
interpolation to the range of the neighboring values.

###  IV.3 Comparison with measurements

For our comparison with measurements results, we scale the quantities that are
a function of the revolution frequency by a factor of $2\pi R/C$ where $R$ is
the bending radius and $C$ is the circumference of ANKA. This way measurements
are comparable to the simulations – which assume an isomagnetic ring with the
same bending radius. The parameters used here are listed in Table 4. Note that
this is just one set of possible parameters because the magnet optics (and
thus $f_{s}$) as well as the RF voltage $V_{RF}$ can be gradually changed at
ANKA.

Table 4: Parameters of an isomagnetic accelerator comparable to ANKA ($t_{d},h,f_{s}$, and $f_{rev}$ are scaled by $2\pi R/C=0.316$). Parameter | Symbol | Value
---|---|---
Beam energy | $E_{0}$ | $1.285\,\mathrm{GeV}$
Energy spread | $\sigma_{\delta,0}$ | $0.47\times 10^{-3}$
Damping time | $t_{d}$ | $3.353\,\mathrm{ms}$
Harmonic number | $h$ | $58.21$
RF voltage | $V_{RF}$ | $1048\,\mathrm{kV}$
Revolution frequency | $f_{rev}$ | $8.582\,\mathrm{MHz}$
Synchrotron frequency | $f_{s}$ | $28.13\,\mathrm{kHz}$
Vacuum chamber height | $g$ | $32\,\mathrm{mm}$

Figure 8:  Example for a simulated (left) and a measured (right) bursting
spectrogram. For the simulated spectrogram the axis are scaled with a factor
of $2\pi R/C$ to correct the mismatch due to the isomagnetic approximation.
There are small differences, e.g. in the threshold current and in the
frequencies. However, keeping in mind the very simple model, the general
structure of the spectrograms matches quite well: There is an isolated finger
pointing down (at $f\approx 35\,\mathrm{kHz}$); the fingertip ($I\approx
0.21\,\mathrm{mA}$) marks the instability threshold. For slightly higher
currents, fluctuations in the lower frequency range ($f<10\,\mathrm{kHz}$)
start and the finger broadens. For the highest currents displayed here, there
is a regime showing parallel frequency lines that stay approximately constant
with changing current.

There are some effects we neglect for the simulation such as the frequency
response of the used detector, and the coherent tune shift observed in the
measurement. Also, contributions from other impedances than shielded CSR (e.g.
geometric impedance) are not studied.

We do separate simulation runs for about 150 different currents between
$I=1.3\,\mathrm{mA}$ and $I=0.5\,\mathrm{mA}$. For the first one, we start
with the highest current and a Gaussian charge distribution that is
significantly broader than the expected distribution. To compensate for this,
we allow some extra time for convergence. For the following simulation runs,
we take the final charge distribution of the run before that has the
(slightly) higher current as starting parameters. (Different approaches to
create starting distributions are discussed in the Appendix A.) For each of
these runs, the simulation time on a AMD Radeon R9 290 graphics card is a bit
more then ten minutes, which makes a total simulation time of about 19 hours.

The spectrum of the emitted CSR is calculated using

$P(t,k)\propto\Re({Z_{k}})\times|\tilde{\rho}(k,t)|^{2},$ (16)

where $\Re({Z_{k}})$ is the real part of the impedance and
$|\tilde{\rho}(k,t)|^{2}$ the form factor of the bunch profile. It is then
integrated to obtain the power a detector would measure

$P(t)\propto\int P(t,k)\mathrm{d}k.$ (17)

In analogy to what is done for the measured data, the resulting signal over
time is Fourier transformed to obtain a spectrogram of the ‘bursting’
frequencies. Figure 8 shows this spectrogram of $P(t)$. The general structure
of the simulation and the measurement results agree very well: The instability
threshold is marked by the occurrence of an isolated finger pointing down down
to $I\approx 0.21\,\mathrm{mA}$). For slightly higher currents, the finger
broadens and fluctuations in the lower frequency range ($f<10\,\mathrm{kHz}$)
start. A third regime is observed at the highest currents. It shows parallel
frequency lines that stay approximately constant with changing current. There
are slight mismatches, e.g. in the threshold currents and in the frequencies,
but most features are well reproduced by the simulation.

This means that for ANKA not only the thresholds (see also Brosi  _et al._
(2016)) but also the dynamics of the micro-bunching instability are governed
by an impedance that can be approximated by the parallel plates CSR impedance.
To explain the details, a more complex model will be needed. Two possibilities
are to take into account higher orders of the momentum compaction factor
$\alpha_{c}$ or additional impedance contributions.

##  V Summary

We introduced _Inovesa_ , a Vlasov-Fokker-Planck-solver that uses a runtime-
optimized implementation of the computation steps. Utilizing OpenCL for
parallelized computation, it can simulate the dynamics of the longitudinal
phase space more than a 150 times faster than a non-optimized implementation –
using a dedicated (consumer-grade) graphics card. Furthermore, we eliminated
sources of numerical artifacts and have done numerical stability studies to
show that relative errors can usually be kept clearly below 1%.

Using _Inovesa_ we were able to simulate the dynamics of the longitudinal
phase space of ANKA in the regime of the micro-bunching instability. To do so,
we used an impedance model that assumes CSR of electrons moving on a circular
path shielded by parallel plates. Considering the simplicity of the model, the
numerical results show an excellent agreement to the measurements.

###### Acknowledgements.

C. Evain is acknowledged for useful discussions, G. Stupakov for his
suggestion to investigate numerical effects of different rotation schemes. We
also thank M. Klein for giving us her code as a reference, and T. Boltz for
beta testing. M. Brosi, P. Schönfeldt, and J. L. Steinmann want to acknowledge
the support by the Helmholtz International Research School for Teratronics
(HIRST).

##  Appendix A Starting Distribution

In equilibrium the energy is distributed according to a Gaussian function and
the bunch profile is described by the Haïssinski distribution Haissinski
(1973). However, here we are interested in the dynamics of the micro-bunching
instability above the threshold current, which means in non-equilibrium. In
this physical state, there is no simple one-dimensional function for any
possible charge distribution.

One possibility is to start the simulation with a Gaussian distribution that
is broader than the expected charge distribution. It will damp down until the
physical state of the instability is reached. Although for currents well above
the instability threshold any possible starting distribution will reach the
same state, we chose these initial conditions because it shows good
convergence. When alternatively using narrower distributions unphysical
structures form and might persist for a long time.

Figure 9:  Evolution of the RMS energy spread over a time of 500 synchrotron
periods for $I=1.5\,\mathrm{mA}$. The width of the line is caused by high
frequency oscillations. For the starting distributions the Haïssinski
distribution (equilibrium at $I=675\,\mathrm{\mu A}$), a Gaussian distribution
with $\sigma_{p}=\sigma_{q}=2.25$, and the final distribution of the previous
run (“adiabatic”, $I=1.54\,\mathrm{mA}$) were used. Note that the Haïssinski
distribution is immediately blown up and becomes larger than the Gaussian
distribution, and that the high frequency oscillation is systematically higher
until $T\approx 300\,\mathrm{T_{s}}$.

In Fig. 9 the evolution of the RMS energy spread over 500 synchrotron periods
is shown for three different starting distributions at $I=1.5\,\mathrm{mA}$.
When using the final distribution of a previous run with slightly higher
current (here $I=1.54\,\mathrm{mA}$), the simulation converges quasi
instantaneously. As shown, the Haïssinski distribution is immediately blown up
and becomes larger than the Gaussian distribution that has been set to be
larger than the expected distribution. Also in the beginning
($t<300\,\mathrm{T_{s}}$) the oscillation is systematically enlarged.

##  Appendix B Running Inovesa

_Inovesa_ is implemented as a non-interactive command-line tool. To
automatically simulate all data for a spectrogram as shown in Fig. 8, a script
may be used. Using the “adiabatic” method discussed in Appendix A, an example
bash script reads:

⬇

1config=”inovesa-run123.cfg”

2lasti=”500”

3./inovesa -I ${last}e-6 –config $config -T 1500 -o $lasti.h5

4for curri in {500..100..5}

5do

6 ./inovesa -I ${curri}e-6 -i $lasti.h5 -c $config -T 500 -o $curri.h5

7 lasti=$curri

8done

The parameters used here are:

\--BunchCurrent

(or -I) for the ring current due to a single bunch given (in Ampere)

\--config

(or -c) the file name of a configuration file

\--rotations

(or -T) the total simulation time (in synchrotron period lengths $T_{s}$))

\--InitialDistFile

(or -i) the file name of an _Inovesa_ result file to use for the initial
particle distribution

\--output

(or -o) the file name to save results to

The configuration file used for the script contains all relevant parameters in
a key-value-representation. Comments may be added using “#”. Here is an
example configuration file:

⬇

BeamEnergy=1.3e+09 # in eV

BeamEnergySpread=0.00047 # relative

BendingRadius=5.559 # in m

BunchCurrent=1.2e-05 # in mA

DampingTime=0.01 # in ms

GridSize=256 # grid points per axis

HarmonicNumber=184 # f_RF/f_rev

PhaseSpaceSize=12 # sigma_z/E per axis

AcceleratingVoltage=1.4e+06 # in V

RevolutionFrequency=2.7e+06 # in Hz

alpha0=2e-4 # momentum compaction factor

VacuumGap=0.032 # full distance in m

gui=true # show live preview of results?

outstep=50 # write output every N steps

padding=8 # factor for zero-pading (FFT)

All parameters are optional: If a parameter is not set _Inovesa_ will fall
back to default values. You might also overwrite settings from a configuration
file by passing the same parameter as a command line argument. For short tests
it is a good idea to enable the live preview (-g true) and not to save the
results (by seting -o /dev/null).

## References

  * Carr  _et al._ (1999) G. Carr, S. Kramer, J. Murphy, J. LaVeign, R. Lobo, D. Reitze, and D. Tanner, in  _Particle Accelerator Conference, 1999. Proceedings of the 1999_ , Vol. 1 (IEEE, 1999) pp. 134–136.
  * Carr  _et al._ (2001) G. Carr, S. Kramer, J. Murphy, R. Lobo, and D. Tanner, Nuclear Instruments and Methods in Physics Research A 463, 387–392 (2001).
  * Müller (2010) A.-S. Müller, Reviews of Accelerator Science and Technology 03, 165 (2010).
  * Abo-Bakr  _et al._ (2003) M. Abo-Bakr, J. Feikes, K. Holldack, P. Kuske, and G. Wustefeld, in  _Particle Accelerator Conference, 2003. PAC 2003. Proceedings of the_ , Vol. 5 (IEEE, 2003) pp. 3023–3025.
  * Shields  _et al._ (2012) W. Shields, R. Bartolini, G. Boorman, P. Karataev, A. Lyapin, J. Puntree, and G. Rehm, Journal of Physics: Conference Series 357, 012037 (2012).
  * Wüstefeld  _et al._ (2010) G. Wüstefeld, J. Feikes, M. Hartrott, M. Ries, A. Hoehl, R. Klein, R. Müller, A. Serdyukov, and G. Ulm, in  _Proceedings of the 1st International Particle Accelerator Conference_ (2010) p. 2508.
  * Evain  _et al._ (2012a) C. Evain, J. Barros, A. Loulergue, M. A. Tordeux, R. Nagaoka, M. Labat, L. Cassinari, G. Creff, L. Manceron, J. B. Brubach, P. Roy, and M. E. Couprie, EPL (Europhysics Letters) 98, 40006 (2012a).
  * Venturini and Warnock (2002) M. Venturini and R. Warnock, Physical Review Letters 89, 224802 (2002).
  * Stupakov and Heifets (2002) G. Stupakov and S. Heifets, Phys. Rev. ST Accel. Beams 5, 054402 (2002).
  * Bane  _et al._ (2010) K. L. F. Bane, Y. Cai, and G. Stupakov, Phys. Rev. ST Accel. Beams 13, 104402 (2010).
  * Cai (2011a) Y. Cai, in  _Proceedings of the 2nd International Particle Accelerator Conference_ (2011) p. 3774.
  * Venturini  _et al._ (2005) M. Venturini, R. Warnock, R. Ruth, and J. A. Ellison, Phys. Rev. ST Accel. Beams 8, 014202 (2005).
  * Roussel  _et al._ (2014) E. Roussel, C. Evain, C. Szwaj, and S. Bielawski, Physical Review Special Topics – Accelerators and Beams 17 (2014), 10.1103/PhysRevSTAB.17.010701.
  * Evain  _et al._ (2012b) C. Evain, J. Barros, A. Loulergue, M. A. Tordeux, R. Nagaoka, M. Labat, L. Cassinari, G. Creff, L. Manceron, J. B. Brubach, P. Roy, and M. E. Couprie, EPL (Europhysics Letters) 98, 40006 (2012b).
  * Brosi  _et al._ (2016) M. Brosi, J. L. Steinmann, E. Blomley, E. Bründermann, M. Caselle, N. Hiller, B. Kehrer, M. J. Nasse, L. Rota, M. Schedler, P. Schönfeldt, M. Schuh, M. Schwarz, M. Weber, and A.-S. Müller, ArXiv e-prints, accepted for PRAB (2016), arXiv:1605.00536 [physics.acc-ph] .
  * git (2016) “Inovesa source code,” https://github.com/Inovesa/Inovesa (2016).
  * Warnock and Ellison (2000) R. L. Warnock and J. A. Ellison,  _A General Method for Propagation of the Phase Space Distribution, with Application to the Sawtooth Instability_ , Tech. Rep. SLAC-PUB-8404 (SLAC, 2000).
  * Murphy  _et al._ (1997) J. Murphy, S. Krinsky, and R. Gluckstern, Particle Accelerators 57, 9 (1997).
  * Cai (2011b) Y. Cai, in  _Proceedings of IPAC 2011_ (2011) pp. 3774–3778.
  * Press  _et al._ (2007) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes 3rd Edition: The Art of Scientific Computing_ , 3rd ed. (Cambridge University Press, New York, NY, USA, 2007).
  * fpm (2009) “Fixed point class documentation,” http://www.codeproject.com/Articles/37636/Fixed-Point-Class (2009).
  * ope (2016) “OpenCL website,” https://www.khronos.org/opencl/ (2016).
  * Wolski (2014) A. Wolski,  _Beam Dynamics in High Energy Particle Accelerators_ (Imperial College Press, London, United Kingdom, 2014).
  * fft (2016) “FFTW manual,” http://www.fftw.org/fftw2_doc/fftw_2.html (2016).
  * clf (2016) “clFTW manual,” http://clmathlibraries.github.io/clFFT/ (2016).
  * IEE (2008) IEEE Std. 754-2008 (2008).
  * Haissinski (1973) J. Haissinski, Il Nuovo Cimento 18 B, 72 (1973).