Comparison of Analytical IBS Growth Rates and Emittances Evolution

This example shows a comparison of the performance of BjorkenMtingwaIBS and NagaitsevIBS vs MAD-X when computing emittances evolutions from analytical IBS growth rates.

We will do the comparison by using the case of the CERN LHC, with protons at top energy (aka 6 800 GeV at the moment), including the crossing schemes, which introduce vertical dispersion through the machine.

Note

This example requires the acc-models-lhc repository next to this file, which you can get with the following command:

git clone -b 2023 https://gitlab.cern.ch/acc-models/acc-models-lhc.git --depth 1

import logging
import warnings

from dataclasses import dataclass
from pprint import pprint

import matplotlib.pyplot as plt
import numpy as np

from helpers import _get_dummy_ibs_from_madx_rates, prepare_all

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 setting up the lattice in MAD-X (we are using cpymad for this) and our IBS classes. The following helper functions calls a configuration file from this repository’s tests, sets up the lattice and beam in MAD-X and creates corresponding BjorkenMtingwaIBS and NagaitsevIBS.

params, madx, BM_IBS, NAG_IBS = prepare_all(
    configname="lhc_top_protons",  # fetches the config file with this name
    extrafile="acc-models-lhc/strengths/ATS_Nominal/2023/ats_30cm.madx",  # to have xing
    stdout=False,  # don't want the MAD-X output
)

Setting up MAD-X for lhc_top_protons case
Calling extra file at 'acc-models-lhc/strengths/ATS_Nominal/2023/ats_30cm.madx'

Here params hold the properties we want to be analytically tracking. We can check that crossing angles are present in the MAD-X lattice, and lead to vertical dispersion:

pprint(params)
print(f"Crossing angles in IP1/5: {madx.globals['on_x1']}, {madx.globals['on_x5']} urad")

table = madx.table.twiss.dframe()
plt.figure(figsize=(11, 6))
plt.plot(table.s, table.dx, lw=2, label="Vertical")
plt.plot(table.s, table.dy, lw=2, label="Horizontal")
plt.ylabel(r"$D_{x,y}$ [m]")
plt.xlabel("Longitudinal location [m]")
plt.legend()
plt.tight_layout()
plt.show()

Params(geom_epsx=2.483661433330319e-10,
       geom_epsy=2.483661433330319e-10,
       sig_delta=4.7170567118141664e-05,
       bunch_length=0.0374740568932693)
Crossing angles in IP1/5: 160.0, 160.0 urad

Comparing Analytical Evolution for a Time Period

We will analytically look at the evolution through time by looping, just as in the Bjorken-Mtingwa and Nagaitsev analytical examples, respectively. Let’s do so for 10 hours of beam time, recomputing the growth rates every 30 minutes.

# Duration to track for and frequency at which to re-compute the growth rates, both in [s]
nsecs = 10 * 3_600  # that's 10h
ibs_step = 30 * 60  # re-compute rates every 30min
seconds = np.linspace(0, nsecs, nsecs).astype(int)
beta_rel: float = BM_IBS.beam_parameters.beta_rel


# 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 dataclasses for each approach & store the initial values
madx_tbt = Records.init_zeroes(nsecs)
bm_tbt = Records.init_zeroes(nsecs)
nag_tbt = Records.init_zeroes(nsecs)

madx_tbt.update_at_turn(0, params.geom_epsx, params.geom_epsy, params.sig_delta, params.bunch_length)
bm_tbt.update_at_turn(0, params.geom_epsx, params.geom_epsy, params.sig_delta, params.bunch_length)
nag_tbt.update_at_turn(0, params.geom_epsx, params.geom_epsy, params.sig_delta, params.bunch_length)

With the settings above and this loop, we compute the emittances every second but we only update the growth rates every 30 minutes.

for sec in range(1, nsecs):
    # ----- Potentially re-compute the IBS growth rates ----- #
    if (sec % ibs_step == 0) or (sec == 1):
        print(f"At {sec}s: re-computing growth rates")
        # For MAD-X - call the 'twiss' then 'ibs' commands (computes from beam attributes)
        MAD_IBS = _get_dummy_ibs_from_madx_rates(madx)  # calls twiss, ibs, gets rates and puts them in returned 'dummy IBS'
        MAD_IBS.beam_parameters = BM_IBS.beam_parameters  # we will want to access these
        MAD_IBS.optics = BM_IBS.optics  # we will want to access these
        # For Bjorken-Mtingwa - compute from values at the previous sec
        BM_IBS.growth_rates(
            bm_tbt.epsilon_x[sec - 1],
            bm_tbt.epsilon_y[sec - 1],
            bm_tbt.sig_delta[sec - 1],
            bm_tbt.bunch_length[sec - 1],
        )
        # For Nagaitsev - compute from values at the previous sec
        NAG_IBS.growth_rates(
            nag_tbt.epsilon_x[sec - 1],
            nag_tbt.epsilon_y[sec - 1],
            nag_tbt.sig_delta[sec - 1],
            nag_tbt.bunch_length[sec - 1],
        )

    # ----- Compute the new parameters using current growth rates ----- #
    # For MAD-X - compute from values at previous second and specify dt=1s
    madx_epsx, madx_epsy, madx_sigd, madx_bl = MAD_IBS.emittance_evolution(
        madx_tbt.epsilon_x[sec - 1],
        madx_tbt.epsilon_y[sec - 1],
        madx_tbt.sig_delta[sec - 1],
        madx_tbt.bunch_length[sec - 1],
        dt=1,
    )
    # For Bjorken-Mtingwa - compute from values at previous second and specify dt=1s
    bm_epsx, bm_epsy, bm_sigd, bm_bl = BM_IBS.emittance_evolution(
        bm_tbt.epsilon_x[sec - 1],
        bm_tbt.epsilon_y[sec - 1],
        bm_tbt.sig_delta[sec - 1],
        bm_tbt.bunch_length[sec - 1],
        dt=1,
    )
    # For Nagaitsev - compute from values at previous second and specify dt=1s
    nag_epsx, nag_epsy, nag_sigd, nag_bl = NAG_IBS.emittance_evolution(
        nag_tbt.epsilon_x[sec - 1],
        nag_tbt.epsilon_y[sec - 1],
        nag_tbt.sig_delta[sec - 1],
        nag_tbt.bunch_length[sec - 1],
        dt=1,
    )

    # ----- Update the records with the new values ----- #
    madx_tbt.update_at_turn(sec, madx_epsx, madx_epsy, madx_sigd, madx_bl)
    bm_tbt.update_at_turn(sec, bm_epsx, bm_epsy, bm_sigd, bm_bl)
    nag_tbt.update_at_turn(sec, nag_epsx, nag_epsy, nag_sigd, nag_bl)

    # ----- Update the beam attributes for MAD-X ----- #
    madx.sequence.lhcb1.beam.ex = madx_epsx
    madx.sequence.lhcb1.beam.ey = madx_epsy
    madx.sequence.lhcb1.beam.sige = madx_sigd * beta_rel**2
    madx.sequence.lhcb1.beam.sigt = madx_bl

At 1s: re-computing growth rates
At 1800s: re-computing growth rates
At 3600s: re-computing growth rates
At 5400s: re-computing growth rates
At 7200s: re-computing growth rates
At 9000s: re-computing growth rates
At 10800s: re-computing growth rates
At 12600s: re-computing growth rates
At 14400s: re-computing growth rates
At 16200s: re-computing growth rates
At 18000s: re-computing growth rates
At 19800s: re-computing growth rates
At 21600s: re-computing growth rates
At 23400s: re-computing growth rates
At 25200s: re-computing growth rates
At 27000s: re-computing growth rates
At 28800s: re-computing growth rates
At 30600s: re-computing growth rates
At 32400s: re-computing growth rates
At 34200s: re-computing 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))

# Plotting horizontal emittances
axs["epsx"].plot(seconds / 3600, 1e10 * madx_tbt.epsilon_x, lw=2.5, label="MAD-X")
axs["epsx"].plot(seconds / 3600, 1e10 * bm_tbt.epsilon_x, lw=1.5, label="BjorkenMtingwaIBS")
axs["epsx"].plot(seconds / 3600, 1e10 * nag_tbt.epsilon_x, lw=1, label="NagaitsevIBS")

# Plotting vertical emittances
axs["epsy"].plot(seconds / 3600, 1e10 * madx_tbt.epsilon_y, lw=2.5, label="MAD-X")
axs["epsy"].plot(seconds / 3600, 1e10 * bm_tbt.epsilon_y, lw=1.5, label="BjorkenMtingwaIBS")
axs["epsy"].plot(seconds / 3600, 1e10 * nag_tbt.epsilon_y, lw=2, label="NagaitsevIBS")

# Plotting momentum spread
axs["sigd"].plot(seconds / 3600, 1e4 * madx_tbt.sig_delta, lw=2.5, label="MAD-X")
axs["sigd"].plot(seconds / 3600, 1e4 * bm_tbt.sig_delta, lw=1.5, label="BjorkenMtingwaIBS")
axs["sigd"].plot(seconds / 3600, 1e4 * nag_tbt.sig_delta, lw=1, label="NagaitsevIBS")

# Plotting bunch length
axs["bl"].plot(seconds / 3600, 1e2 * madx_tbt.bunch_length, lw=2.5, label="MAD-X")
axs["bl"].plot(seconds / 3600, 1e2 * bm_tbt.bunch_length, lw=1.5, label="BjorkenMtingwaIBS")
axs["bl"].plot(seconds / 3600, 1e2 * nag_tbt.bunch_length, lw=1, label="NagaitsevIBS")

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

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("Duration [h]")

for axis in axs.values():
    axis.xaxis.set_major_locator(plt.MaxNLocator(8))

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

# Figure parameters
fig.suptitle("LHC Top Protons w/ Xing")
plt.tight_layout()
plt.show()

Notice how we observe a great agreement between different implementations, except for the vertical emittances. This is expected, as the lattice setup includes vertical dispersion which is not taken into consideration by the NagaitsevIBS formalism.

References

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

• analytical: BjorkenMtingwaIBS, NagaitsevIBS, growth_rates, emittance_evolution