Include quantum excitation

[jquentino]

Introduction

This PR contains a set of modifications in the TrackCPP structure to implement the quantum excitation (or quantum diffusion) in our tracking algorithms. Including this effect will make it possible to reach equilibrium conditions by particle tracking. The related PR in Pyaccel is: lnls-fac/pyaccel#113

In trackcpp, the radiation_on property was changed to accept the integers 0, 1, or 2, with each number meaning a different radiation state: no radiation, damping and damping+quantum excitation, respectively.

In Pyaccel, it is possible to set the property with an integer, string, or a boolean (due to backward compatibility), with the following options:

• 0, False, "off" = No radiative effects.
• 1, True, "damping" = Turns on radiation damping, without quantum excitation.
• 2, "full" = Turns on radiation damping with quantum excitation

The random quantum kicks can be generated from a normal distribution or from a uniform distribution. In Pyaccel, this can be controlled by the function pyaccel.utils.set_distribution. It is also possible to set a seed in the pseudo-random number generator with pyaccel.utils.set_random_seed, allowing the repeatability of the simulations. By default, the distribution is 'normal' and the seed is randomly chosen.

In the following sections, I will show some performance tests and the results of some simulations.

Performance Tests

In addition to the inclusion of radiation effects, some pass methods were refactored to improve the tracking performance and reduce the simulation time. Special thanks to @fernandohds564, responsible for the great part of this refactorization, and to @xresende who helped to build a fast random number generator.

The below table shows the meantime to simulate a single particle for 500 turns in the SIRIUS lattice for different configurations of radiative effects. The mean was computed over 10 simulations and the uncertainty corresponds to the standard deviation.

Radiative effects	Time (New pass methods)	Time (Old pass methods)
No radiative effect	660 +/- 50 ms	740 +/- 40 ms
Radiation damping	960 +/- 50 ms	1010 +/- 30 ms
Rad. Damping + Quant. Diffusion (uniform)	1100 +/- 80 ms	---
Rad. Damping + Quant. Diffusion (normal)	1380 +/- 80 ms	---

CPU details:

• RAM: 15.5 GiB
• Processor: Intel® Core™ i7-8700 CPU @ 3.20GHz × 12
• Graphics: Intel® UHD Graphics 630 (CFL GT2)

Convergence to equilibrium parameters

The following Figure shows the beam's second momenta evolution in the SIRIUS model with fitted coupling and vertical dispersion (lnls-fac/pymodels#85). The simulation was carried out with 2000 particles for 25000 turns, with all particles starting from the fixed point $\approx$ [0,0,0,0,0,0]. The equilibrium values (the straight black lines in the figure) were taken from the equilibrium beam envelope matrix.

Normal distributed quantum kicks:

Uniform distributed quantum kicks:

The non-diagonal terms of the beam envelope matrix were observed to oscillate around the equilibrium values. To express the oscillation amplitude in terms of more intuitive quantities, they were normalized by the correspondent diagonal terms (normalized covariance):

$$ R_{ij} = \frac{\sigma^2_{i, j}}{\sqrt{\sigma^2_{i, i} \sigma^2_{j, j}}}, $$

which ranges from -1 to 1. $i$ and $j$ are any 6D phase space coordinates pair. The results are shown below.

Normal distributed quantum kicks:

Uniform distributed quantum kicks:

Excitation rate and damping times fitting.

The above simulation was repeated in a model without coupling and for fewer turns (5000). In the absence of coupling, the beam size growth and the energy deviation evolution are simply described by the following equations:

$$ \frac{d(\sigma_x^2)}{dt} = (\frac{Q_x\beta_x}{2} - \frac{2\sigma_x^2}{\tau_x}) $$

$$ \frac{d(\sigma_{\delta}^2)}{dt} = (\frac{Q_{\epsilon}}{2 E_0^2} - \frac{2\sigma_{\delta}^2}{\tau_\epsilon}) $$

I integrated them and fitted $Q_x$, $Q_{\epsilon}$ and the damping times. The Figures below show the results:

Beam horizontal size	Energy deviation

Comparing the obtained values with the nominal values:

Quantity	Fitting	Nominal
$Q_x$ [ pm/s]	0.532	0.598
$\tau_x$ [ ms]	19.4	16.8
$Q_{\epsilon}$ [MeV^2/s]	2.058	2.029
$\tau_{\epsilon}$ [ms]	16.2	12.8

The low precision in fitting the damping times can be justified by the fact that quantum excitation dominates over radiation damping for the applied initial conditions.

[murilobalves]

good job @jquentino, congrats!

[fernandohds564]

very nice @jquentino ! Congrats! I just have a few questions and suggestions in the code.

[jquentino]

To keep the group members that were on vacation updated about this PR when they return:

In some talks with @fernandohds564, we decide to change the radiation_on property from a boolean to an enumerate, with the following options:

• 0 = no radiation,
• 1 = radiation damping,
• 2 = damping+quantum_excitation.

The idea is that the user can set one of the three options by passing the correspondent integer number or a string.
In this way, we can use radiation_on to control all available radiative effects without creating a new property.

The feature is being developed in the branch modify-rad-properties and will be merged into this PR when it is working.

[xresende]

should we not redefine pnorm = (1 + pos.de) here too?
I am a bit confused: the two code sections seem inconsistent to me.

[fernandohds564]

@xresende, the idea of defining the second appearance of pnorm was to avoid the unnecessary floating point division operation and substitute the other divisions by multiplications.

[xresende]

the optimization idea I understand. but in one case the variable pnorm is set to (1 + pos.de) and in another to 1 / (1 + pos.de). There seems to be a name inconsistency. maybe we should use pnorm_invin one of the cases.

[fernandohds564]

I didn't want to create a new variable. The inconsistency is pointed out in the comment in front of the definition, for clarity.