Nagaitsev Formalism - Analytical Growth Rates and Emittance Evolution

This example shows how to use the NagaitsevIBS class to calculate IBS growth rates and emittances evolutions analytically.

We will demonstrate using an xtrack.Line of the CLIC damping ring, for a positron beam.

import logging
import warnings

from dataclasses import dataclass

import matplotlib.pyplot as plt
import numpy as np
import xpart as xp
import xtrack as xt

from xibs.analytical import NagaitsevIBS
from xibs.formulary import _bunch_length, _geom_epsx, _geom_epsy, _sigma_delta
from xibs.inputs import BeamParameters, OpticsParameters

warnings.simplefilter("ignore")  # for this tutorial's clarity
logging.basicConfig(
    level=logging.WARNING,
    format="[%(asctime)s] [%(levelname)s] - %(module)s.%(funcName)s:%(lineno)d - %(message)s",
    datefmt="%H:%M:%S",
)
plt.rcParams.update(
    {
        "font.family": "serif",
        "font.size": 20,
        "axes.titlesize": 20,
        "axes.labelsize": 20,
        "xtick.labelsize": 20,
        "ytick.labelsize": 20,
        "legend.fontsize": 15,
        "figure.titlesize": 20,
    }
)

Let’s start by defining the line and particle information, as well as some parameters for later use:

line_file = "lines/chrom-corr_DR.newlattice_2GHz.json"
bunch_intensity = 4.4e9
sigma_z = 1.58e-3
nemitt_x = 5.6644e-07
nemitt_y = 3.7033e-09
n_part = int(5e3)  # for this example let's not use too many particles

Setting up line and particles

Let’s start by loading the xtrack.Line, activating RF cavities and generating a matched Gaussian bunch:

line = xt.Line.from_json(line_file)
line.build_tracker(extra_headers=["#define XTRACK_MULTIPOLE_NO_SYNRAD"])

# ----- Power accelerating cavities ----- #
cavities = [element for element in line.elements if isinstance(element, xt.Cavity)]
for cavity in cavities:
    cavity.lag = 180

p0 = xp.Particles(mass0=xp.ELECTRON_MASS_EV, q0=1, p0c=2.86e9)
line.particle_ref = p0
twiss = line.twiss()

particles = xp.generate_matched_gaussian_bunch(
    num_particles=n_part,
    total_intensity_particles=bunch_intensity,
    nemitt_x=nemitt_x,
    nemitt_y=nemitt_y,
    sigma_z=sigma_z,
    particle_ref=p0,
    line=line,
)

Loading line from dict:   0%|          | 0/71796 [00:00<?, ?it/s]
Loading line from dict:   4%|▍         | 3033/71796 [00:00<00:02, 30325.98it/s]
Loading line from dict:   9%|▊         | 6142/71796 [00:00<00:02, 30770.65it/s]
Loading line from dict:  13%|█▎        | 9266/71796 [00:00<00:02, 30979.79it/s]
Loading line from dict:  17%|█▋        | 12364/71796 [00:00<00:01, 30967.85it/s]
Loading line from dict:  22%|██▏       | 15518/71796 [00:00<00:01, 31163.69it/s]
Loading line from dict:  26%|██▌       | 18635/71796 [00:00<00:02, 24490.33it/s]
Loading line from dict:  30%|███       | 21746/71796 [00:00<00:01, 26310.30it/s]
Loading line from dict:  35%|███▍      | 24855/71796 [00:00<00:01, 27661.10it/s]
Loading line from dict:  39%|███▊      | 27751/71796 [00:00<00:01, 27996.62it/s]
Loading line from dict:  43%|████▎     | 30643/71796 [00:01<00:01, 27463.73it/s]
Loading line from dict:  47%|████▋     | 33454/71796 [00:01<00:01, 27195.63it/s]
Loading line from dict:  50%|█████     | 36218/71796 [00:01<00:01, 27027.26it/s]
Loading line from dict:  55%|█████▍    | 39331/71796 [00:01<00:01, 28213.32it/s]
Loading line from dict:  59%|█████▉    | 42464/71796 [00:01<00:01, 29119.27it/s]
Loading line from dict:  64%|██████▎   | 45651/71796 [00:01<00:00, 29929.40it/s]
Loading line from dict:  68%|██████▊   | 48878/71796 [00:01<00:00, 30621.57it/s]
Loading line from dict:  73%|███████▎  | 52082/71796 [00:01<00:00, 31042.16it/s]
Loading line from dict:  77%|███████▋  | 55368/71796 [00:01<00:00, 31578.41it/s]
Loading line from dict:  82%|████████▏ | 58533/71796 [00:02<00:00, 24344.20it/s]
Loading line from dict:  86%|████████▌ | 61728/71796 [00:02<00:00, 26221.51it/s]
Loading line from dict:  90%|████████▉ | 64570/71796 [00:02<00:00, 26433.23it/s]
Loading line from dict:  94%|█████████▍| 67369/71796 [00:02<00:00, 26381.33it/s]
Loading line from dict:  98%|█████████▊| 70116/71796 [00:02<00:00, 26369.81it/s]
Loading line from dict: 100%|██████████| 71796/71796 [00:02<00:00, 27821.33it/s]
Done loading line from dict.
Compiling ContextCpu kernels...
Done compiling ContextCpu kernels.

We can have a look at the generated particles (see the xsuite user guide for more)

fig, (axx, axy, axz) = plt.subplots(3, 1, figsize=(9, 11))

axx.plot(1e6 * particles.x, 1e5 * particles.px, ".", ms=3)
axy.plot(1e6 * particles.y, 1e6 * particles.py, ".", ms=3)
axz.plot(1e3 * particles.zeta, 1e3 * particles.delta, ".", ms=3)

axx.set_xlabel(r"$x$ [$\mu$m]")
axx.set_ylabel(r"$p_x$ [$10^{-5}$]")
axy.set_xlabel(r"$y$ [$\mu$m]")
axy.set_ylabel(r"$p_y$ [$10^{-3}$]")
axz.set_xlabel(r"$z$ [$10^{-3}$]")
axz.set_ylabel(r"$\delta$ [$10^{-3}$]")
fig.align_ylabels([axx, axy, axz])
plt.tight_layout()
plt.show()

We can compute initial (geometrical) emittances as well as the bunch length from the xtrack.Particles object:

geom_epsx = _geom_epsx(particles, twiss.betx[0], twiss.dx[0])
geom_epsy = _geom_epsy(particles, twiss.bety[0], twiss.dy[0])
bunch_l = _bunch_length(particles)
sig_delta = _sigma_delta(particles)

Computing Elliptic Integrals and IBS Growth Rates

Let us instantiate the NagaitsevIBS class. It is fairly simple and is done from both the beam parameters of the particles and the optics parameters of the line. For each of these a specific dataclass is provided in the inputs module.

beam_params = BeamParameters(particles)
optics = OpticsParameters(twiss)
IBS = NagaitsevIBS(beam_params, optics)

As a first step, all calculations in rely on the computing of Nagaitsev integrals [Nag05], which is done with a dedicated function. These are returned, but also stored internally in the IBS object and will be updated internally each time they are computed.

integrals = IBS.integrals(geom_epsx, geom_epsy, sig_delta)
print(integrals)
print(IBS.elliptic_integrals)

NagaitsevIntegrals(Ix=8.086835088033967e+19, Iy=1.5092530205187414e+17, Iz=6.0325356795434126e+23)
NagaitsevIntegrals(Ix=8.086835088033967e+19, Iy=1.5092530205187414e+17, Iz=6.0325356795434126e+23)

From these the IBS growth rates can be computed, which is again done by calling a dedicated function. If the integrals mentioned above have not been computed yet, an message will be logged for the user and they will be computed first. Once again, these are returned but also stored internally and updated internally each time they are computed.

growth_rates = IBS.growth_rates(geom_epsx, geom_epsy, sig_delta, bunch_l)
print(growth_rates)
print(IBS.ibs_growth_rates)

IBSGrowthRates(Tx=370.6264156506234, Ty=106.49624313657692, Tz=87.01055361087455)
IBSGrowthRates(Tx=370.6264156506234, Ty=106.49624313657692, Tz=87.01055361087455)

Please note that the IBS.growth_rates method by default re-computes the integrals before computing the growth rates. This is convenient as usually, one needs to update the integrals when they want to update the growth rates. It can be disabled by setting the compute_integrals argument to False.

Computing New Emittances from Growth Rates

From these one can compute the emittances at the next time step. For the emittances at the next turn, once should use \(1 / f_{rev}\) as the time step, which is the default value used if it is not provided.

new_geom_epsx, new_geom_epsy, new_sig_delta, new_bunch_length = IBS.emittance_evolution(
    geom_epsx, geom_epsy, sig_delta, bunch_l
)
print(f"Next time step geometrical epsilon_x = {new_geom_epsx} m")
print(f"Next time step geometrical epsilon_y = {new_geom_epsy} m")
print(f"Next time step sigma_delta = {new_sig_delta}")
print(f"Next time step bunch length = {new_bunch_length} m")

Next time step geometrical epsilon_x = 9.863604502224669e-11 m
Next time step geometrical epsilon_y = 6.404092907988679e-13 m
Next time step sigma_delta = 0.0017699984202813742
Next time step bunch length = 0.0016073564228652843 m

Analytical Evolution for Many Turns

One can then analytically look at the evolution through many turns by looping over this calculation. Let’s do this for 10 000 turns, re-computing the elliptic integrals and the IBS growth rates every 200 turns. The more frequent this update the more physically accurate the results will be, but the longer the simulation as this computation is the most compute-intensive process.

nturns = 10_000  # number of turns to loop for
ibs_step = 200  # frequency at which to re-compute the growth rates in [turns]
turns = np.arange(nturns, dtype=int)  # array of turns


# Set up a dataclass to store the results
@dataclass
class Records:
    """Dataclass to store (and update) important values through tracking."""

    epsilon_x: np.ndarray  # geometric horizontal emittance in [m]
    epsilon_y: np.ndarray  # geometric vertical emittance in [m]
    sig_delta: np.ndarray  # momentum spread
    bunch_length: np.ndarray  # bunch length in [m]

    @classmethod
    def init_zeroes(cls, n_turns: int):
        return cls(
            epsilon_x=np.zeros(n_turns, dtype=float),
            epsilon_y=np.zeros(n_turns, dtype=float),
            sig_delta=np.zeros(n_turns, dtype=float),
            bunch_length=np.zeros(n_turns, dtype=float),
        )

    def update_at_turn(self, turn: int, epsx: float, epsy: float, sigd: float, bl: float):
        """Works for turns / seconds, just needs the correct index to store in."""
        self.epsilon_x[turn] = epsx
        self.epsilon_y[turn] = epsy
        self.sig_delta[turn] = sigd
        self.bunch_length[turn] = bl


# Initialize the dataclass & store the initial values
turn_by_turn = Records.init_zeroes(nturns)
turn_by_turn.update_at_turn(0, geom_epsx, geom_epsy, sig_delta, bunch_l)


# ----- We loop here now ----- #

for turn in range(1, nturns):
    # ----- Potentially re-compute the elliptic integrals and IBS growth rates ----- #
    if (turn % ibs_step == 0) or (turn == 1):
        print(f"Turn {turn}: re-computing the Nagaitsev integrals and growth rates")
        # We compute from values at the previous turn
        IBS.growth_rates(  # this recomputes the integrals by default
            turn_by_turn.epsilon_x[turn - 1],
            turn_by_turn.epsilon_y[turn - 1],
            turn_by_turn.sig_delta[turn - 1],
            turn_by_turn.bunch_length[turn - 1],
        )

    # ----- Compute the new emittances ----- #
    new_emit_x, new_emit_y, new_sig_delta, new_bunch_length = IBS.emittance_evolution(
        epsx=turn_by_turn.epsilon_x[turn - 1],
        epsy=turn_by_turn.epsilon_y[turn - 1],
        sigma_delta=turn_by_turn.sig_delta[turn - 1],
        bunch_length=turn_by_turn.bunch_length[turn - 1],
        # dt = 1.0 / IBS.optics.revolution_frequency,  # default value
    )

    # ----- Update the records with the new values ----- #
    turn_by_turn.update_at_turn(turn, new_emit_x, new_emit_y, new_sig_delta, new_bunch_length)

Turn 1: re-computing the Nagaitsev integrals and growth rates
Turn 200: re-computing the Nagaitsev integrals and growth rates
Turn 400: re-computing the Nagaitsev integrals and growth rates
Turn 600: re-computing the Nagaitsev integrals and growth rates
Turn 800: re-computing the Nagaitsev integrals and growth rates
Turn 1000: re-computing the Nagaitsev integrals and growth rates
Turn 1200: re-computing the Nagaitsev integrals and growth rates
Turn 1400: re-computing the Nagaitsev integrals and growth rates
Turn 1600: re-computing the Nagaitsev integrals and growth rates
Turn 1800: re-computing the Nagaitsev integrals and growth rates
Turn 2000: re-computing the Nagaitsev integrals and growth rates
Turn 2200: re-computing the Nagaitsev integrals and growth rates
Turn 2400: re-computing the Nagaitsev integrals and growth rates
Turn 2600: re-computing the Nagaitsev integrals and growth rates
Turn 2800: re-computing the Nagaitsev integrals and growth rates
Turn 3000: re-computing the Nagaitsev integrals and growth rates
Turn 3200: re-computing the Nagaitsev integrals and growth rates
Turn 3400: re-computing the Nagaitsev integrals and growth rates
Turn 3600: re-computing the Nagaitsev integrals and growth rates
Turn 3800: re-computing the Nagaitsev integrals and growth rates
Turn 4000: re-computing the Nagaitsev integrals and growth rates
Turn 4200: re-computing the Nagaitsev integrals and growth rates
Turn 4400: re-computing the Nagaitsev integrals and growth rates
Turn 4600: re-computing the Nagaitsev integrals and growth rates
Turn 4800: re-computing the Nagaitsev integrals and growth rates
Turn 5000: re-computing the Nagaitsev integrals and growth rates
Turn 5200: re-computing the Nagaitsev integrals and growth rates
Turn 5400: re-computing the Nagaitsev integrals and growth rates
Turn 5600: re-computing the Nagaitsev integrals and growth rates
Turn 5800: re-computing the Nagaitsev integrals and growth rates
Turn 6000: re-computing the Nagaitsev integrals and growth rates
Turn 6200: re-computing the Nagaitsev integrals and growth rates
Turn 6400: re-computing the Nagaitsev integrals and growth rates
Turn 6600: re-computing the Nagaitsev integrals and growth rates
Turn 6800: re-computing the Nagaitsev integrals and growth rates
Turn 7000: re-computing the Nagaitsev integrals and growth rates
Turn 7200: re-computing the Nagaitsev integrals and growth rates
Turn 7400: re-computing the Nagaitsev integrals and growth rates
Turn 7600: re-computing the Nagaitsev integrals and growth rates
Turn 7800: re-computing the Nagaitsev integrals and growth rates
Turn 8000: re-computing the Nagaitsev integrals and growth rates
Turn 8200: re-computing the Nagaitsev integrals and growth rates
Turn 8400: re-computing the Nagaitsev integrals and growth rates
Turn 8600: re-computing the Nagaitsev integrals and growth rates
Turn 8800: re-computing the Nagaitsev integrals and growth rates
Turn 9000: re-computing the Nagaitsev integrals and growth rates
Turn 9200: re-computing the Nagaitsev integrals and growth rates
Turn 9400: re-computing the Nagaitsev integrals and growth rates
Turn 9600: re-computing the Nagaitsev integrals and growth rates
Turn 9800: re-computing the Nagaitsev integrals and growth rates

Feel free to run this simulation for more turns, or with a different frequency of the IBS growth rates re-computation. After this is done running, we can plot the evolutions across the turns:

fig, axs = plt.subplot_mosaic([["epsx", "epsy"], ["sigd", "bl"]], sharex=True, figsize=(13, 7))

axs["epsx"].plot(turns, 1e10 * turn_by_turn.epsilon_x, lw=2)
axs["epsy"].plot(turns, 1e13 * turn_by_turn.epsilon_y, lw=2)
axs["sigd"].plot(turns, 1e3 * turn_by_turn.sig_delta, lw=2)
axs["bl"].plot(turns, 1e3 * turn_by_turn.bunch_length, lw=2)

# Axes parameters
axs["epsx"].set_ylabel(r"$\varepsilon_x$ [$10^{-10}$m]")
axs["epsy"].set_ylabel(r"$\varepsilon_y$ [$10^{-13}$m]")
axs["sigd"].set_ylabel(r"$\sigma_{\delta}$ [$10^{-3}$]")
axs["bl"].set_ylabel(r"Bunch length [mm]")

for axis in (axs["epsy"], axs["bl"]):
    axis.yaxis.set_label_position("right")
    axis.yaxis.tick_right()

for axis in (axs["sigd"], axs["bl"]):
    axis.set_xlabel("Turn Number")

fig.align_ylabels((axs["epsx"], axs["sigd"]))
fig.align_ylabels((axs["epsy"], axs["bl"]))

plt.tight_layout()
plt.show()

References

The use of the following functions, methods, classes and modules is shown in this example:

• analytical: NagaitsevIBS, growth_rates, integrals, emittance_evolution

• inputs: BeamParameters, OpticsParameters