Intra-Beam Scattering
Analytical Growth Rates
The following example illustrates how to obtain Intra-Beam Scattering growth rates in Xsuite.
The functionality is exposed directly through xtrack.TwissTable and can make use of two different formalism: Nagaitsev
and Bjorken-Mtingwa
.
The former provides a computationally efficient approach but does not account for vertical dispersion, while the latter correctly accounts for it but is slower.
See also: xtrack.twiss.TwissTable.get_ibs_growth_rates()
import json
import xtrack as xt
##########################
# Load xt.Line from file #
##########################
fname_line_particles = "../../../xtrack/test_data/lhc_no_bb/" \
"line_and_particle.json"
with open(fname_line_particles, "r") as fid:
input_data = json.load(fid)
line = xt.Line.from_json(fname_line_particles)
line.particle_ref = xt.Particles.from_dict(input_data["particle"])
tw = line.twiss(method="4d")
#####################
# Define parameters #
#####################
# Line is for LHC protons at top energy
bunch_intensity: int = int(1.8e11)
nemitt_x: float = 1.8e-6
nemitt_y: float = 1.8e-6
sigma_delta: float = 4.71e-5
bunch_length: float = 3.75e-2
###################################
# Get growth rates with Nagaitsev #
###################################
nag_growth_rates = tw.get_ibs_growth_rates(
formalism="nagaitsev",
total_beam_intensity=bunch_intensity,
nemitt_x=nemitt_x,
nemitt_y=nemitt_y,
sigma_delta=sigma_delta,
bunch_length=bunch_length,
bunched=True,
)
#########################################
# Get growth rates with Bjorken-Mtingwa #
#########################################
bm_growth_rates = tw.get_ibs_growth_rates(
formalism="bjorken-mtingwa", # also accepts "b&m"
total_beam_intensity=bunch_intensity,
nemitt_x=nemitt_x,
nemitt_y=nemitt_y,
sigma_delta=sigma_delta,
bunch_length=bunch_length,
bunched=True,
)
##########################################################
# Compare: we expect Nagaitsev to be wrong in horizontal #
##########################################################
print()
print("Computed from normalized emittances:")
print("------------------------------------")
print(f"Nagaitsev: {nag_growth_rates}")
print(f"Bjorken-Mtingwa: {bm_growth_rates}")
# Computed from normalized emittances:
# ------------------------------------
# Nagaitsev: IBSAmplitudeGrowthRates(Kx=3.12e-05, Ky=-1.14e-09, Kz=0.000155)
# Bjorken-Mtingwa: IBSAmplitudeGrowthRates(Kx=3.11e-05, Ky=5.52e-07, Kz=0.000155)
# Complete source: xfields/examples/005_ibs/001_growth_rates_with_vdisp.py
Amplitude and Emittance Conventions
For consistency with the computed Synchrotron Radiation damping times, in Xsuite the IBS amplitude growth rates are computed. The following short example shows how to switch between the amplitude and emittance growth rates should one need to.
import xtrack as xt
##########################
# Load xt.Line from file #
##########################
fname_line_particles = "../../../xtrack/test_data/sps_ions/line_and_particle.json"
line = xt.Line.from_json(fname_line_particles)
tw = line.twiss(method="4d")
#####################
# Define parameters #
#####################
# Line is for SPS ions at injection
bunch_intensity: int = int(3.5e8)
nemitt_x: float = 1.2612e-6
nemitt_y: float = 0.9081e-6
sigma_delta: float = 3.59e-4
bunch_length: float = 19.51e-2
####################
# Get growth rates #
####################
# There is no vertical dispersion so Nagaitsev
# will be correct in vertical
amp_growth_rates = tw.get_ibs_growth_rates(
formalism="nagaitsev",
total_beam_intensity=bunch_intensity,
nemitt_x=nemitt_x,
nemitt_y=nemitt_y,
sigma_delta=sigma_delta,
bunch_length=bunch_length,
bunched=True,
)
##########################################################
# Converting between amplitude and emittance conventions #
##########################################################
# Notice how, when printing the returned object, it states
# the growth rates are given in amplitude convention
print(amp_growth_rates)
# IBSAmplitudeGrowthRates(Kx=0.000518, Ky=0.00552, Kz=0.00402)
# Methods are implemented to convert to the emittance convention
emit_growth_rates = amp_growth_rates.to_emittance_growth_rates()
print(emit_growth_rates)
# IBSEmittanceGrowthRates(Kx=0.00104, Ky=0.011, Kz=0.00803)
# It is also possible to convert back to the amplitude convention
print(f"Initial: {amp_growth_rates}")
print(f"Converted twice: {emit_growth_rates.to_amplitude_growth_rates()}")
# Initial: IBSAmplitudeGrowthRates(Kx=0.000518, Ky=0.00552, Kz=0.00402)
# Converted twice: IBSAmplitudeGrowthRates(Kx=0.000518, Ky=0.00552, Kz=0.00402)
####################################################
# Converting between growth rates and growth times #
####################################################
# Should one want the growth times, a method is available in both
# conventions to perform this conversion, although it returns a tuple
print(f"Amp times from amp rates: {amp_growth_rates.to_amplitude_growth_times()}")
print(f"Amp times from emit rates: {emit_growth_rates.to_amplitude_growth_times()}")
# Amp times from amp rates: (1930.7146824847905, 181.11747760500302, 248.968512633387)
# Amp times from emit rates: (1930.7146824847905, 181.11747760500302, 248.968512633387)
# And it is of course possible to get the emittance
# growth times from any of the two conventions
print(f"Emit times from amp rates: {amp_growth_rates.to_emittance_growth_times()}")
print(f"Emit times from emit rates: {emit_growth_rates.to_emittance_growth_times()}")
# Emit times from amp rates: (965.3573412423953, 90.55873880250151, 124.4842563166935)
# Emit times from emit rates: (965.3573412423953, 90.55873880250151, 124.4842563166935)
# Complete source: xfields/examples/005_ibs/002_growth_rates_conventions.py
Steady State Emittances in the Presence of Synchrotron Radiation, Quantum Excitation and Intra-Beam Scattering
The steady-state emittances in the presence of Synchrotron Radiation (SR), Quantum Excitation (QE), and Intra-Beam Scattering (IBS) emerge from a dynamic equilibrium, where the combined effect of these three phenomena balances each other out. These emittances can be calculated in Xsuite by numerically solving a system of ordinary differential equations while enforcing constraints on the transverse emittances. The ODE solved by the function are detailed in the Physics guide.
See also: xtrack.twiss.TwissTable.get_ibs_and_synrad_emittance_evolution()
In the following example, steady state emittances are calculated in the presence of a coupling constraint between transverse planes. Notice how the coupling constraint of round beams (coupling factor = 1) is respected through the evolution of the emittances to the steady-state.
import matplotlib.pyplot as plt
import xtrack as xt
from scipy.constants import e
##########################
# Load xt.Line from file #
##########################
fname_line_particles = "../../../xtrack/test_data/bessy3/bessy3.json"
line = xt.Line.from_json(fname_line_particles) # has particle_ref
line.build_tracker()
########################
# Twiss with Radiation #
########################
# We need to Twiss with Synchrotron Radiation enabled to obtain
# the SR equilibrium emittances and damping constants
line.matrix_stability_tol = 1e-2
line.configure_radiation(model="mean")
line.compensate_radiation_energy_loss()
tw = line.twiss(eneloss_and_damping=True)
######################################
# Steady-State Emittance Calculation #
######################################
bunch_intensity = 1e-9 / e # 1C bunch charge
emittance_coupling_factor = 1 # round beam
# If not providing starting emittances, the function will
# default to the SR equilibrium emittances from the TwissTable
result = tw.get_ibs_and_synrad_emittance_evolution(
formalism="nagaitsev", # can also be "bjorken-mtingwa"
total_beam_intensity=bunch_intensity,
# gemitt_x=..., # defaults to tw.eq_gemitt_x
# gemitt_y=..., # defaults to tw.eq_gemitt_x
# gemitt_zeta=..., # defaults to tw.eq_gemitt_zeta
emittance_coupling_factor=emittance_coupling_factor,
emittance_constraint="coupling",
)
# The returned object is an xtrack Table
print(result)
# Table: 1089 rows, 9 cols
# time gemitt_x nemitt_x gemitt_y nemitt_y gemitt_zeta nemitt_zeta sigma_zeta sigma_delta Kx Ky Kz
# 1.2104374139405318e-06 6.8355e-11 3.34418e-07 6.8355e-11 3.34418e-07 3.38167e-06 0.0581296 0.00344698 0.000981052 58.4838 -0.0165806 17.0318
# 9.306647341196751e-05 6.87219e-11 3.36214e-07 6.87219e-11 3.36214e-07 3.39225e-06 0.0583114 0.00345237 0.000982586 58.4736 -0.0165775 17.0292
# 0.0001849225094099945 6.90808e-11 3.37969e-07 6.90808e-11 3.37969e-07 3.40267e-06 0.0584907 0.00345767 0.000984095 57.7058 -0.0163509 16.8361
# ...
# 0.0997568655312697 8.4492e-11 4.13367e-07 8.4492e-11 4.13367e-07 4.51774e-06 0.0776582 0.00398413 0.00113393 29.799 -0.00861232 8.14787
# 0.09984872156726772 8.4492e-11 4.13367e-07 8.4492e-11 4.13367e-07 4.51774e-06 0.0776583 0.00398413 0.00113393 29.799 -0.00861232 8.14786
# 0.09994057760326575 8.44919e-11 4.13367e-07 8.44919e-11 4.13367e-07 4.51775e-06 0.0776584 0.00398414 0.00113393 29.799 -0.00861232 8.14784
######################################
# Comparison with analytical results #
######################################
# These are analytical estimate (from the last step's IBS growth rates)
# The factor below is to be respected with the coupling constraint
factor = 1 + emittance_coupling_factor * (tw.partition_numbers[1] / tw.partition_numbers[0])
analytical_x = result.gemitt_x[0] / (1 - result.Kx[-1] / (tw.damping_constants_s[0] * factor))
analytical_y = result.gemitt_y[0] / (1 - result.Kx[-1] / (tw.damping_constants_s[0] * factor))
analytical_z = result.gemitt_zeta[0] / (1 - result.Kz[-1] / (tw.damping_constants_s[2]))
print()
print("Emittance Constraint: Coupling")
print("Horizontal steady-state emittance:")
print("---------------------------------")
print(f"Analytical: {analytical_x}")
print(f"ODE: {result.eq_sr_ibs_gemitt_x}")
print("Vertical steady-state emittance:")
print("-------------------------------")
print(f"Analytical: {analytical_y}")
print(f"ODE: {result.eq_sr_ibs_gemitt_y}")
print("Longitudinal steady-state emittance:")
print("-----------------------------------")
print(f"Analytical: {analytical_z}")
print(f"ODE: {result.eq_sr_ibs_gemitt_zeta}")
# Emittance Constraint: Coupling
# Horizontal steady-state emittance:
# ---------------------------------
# Analytical: 8.450195693021898e-11
# ODE: 8.449194372183254e-11
# Vertical steady-state emittance:
# -------------------------------
# Analytical: 8.450195693021898e-11
# ODE: 8.449194372183254e-11
# Longitudinal steady-state emittance:
# -----------------------------------
# Analytical: 4.518939710729963e-06
# ODE: 4.5177482475985255e-06
# The results from the table can easily be plotted to view
# at the evolution of various parameters across time steps
# Complete source: xfields/examples/005_ibs/005_steady_state_emittances_coupling.py

This example, quite similar, shows how to do the same but with an excitation constraint between the transvserse planes. Notice how this time a specific factor between transverse emittances is respected through their evolution to the steady-state.
import matplotlib.pyplot as plt
import xtrack as xt
from scipy.constants import e
##########################
# Load xt.Line from file #
##########################
fname_line_particles = "../../../xtrack/test_data/bessy3/bessy3.json"
line = xt.Line.from_json(fname_line_particles) # has particle_ref
line.build_tracker()
########################
# Twiss with Radiation #
########################
# We need to Twiss with Synchrotron Radiation enabled to obtain
# the SR equilibrium emittances and damping constants
line.matrix_stability_tol = 1e-2
line.configure_radiation(model="mean")
line.compensate_radiation_energy_loss()
tw = line.twiss(eneloss_and_damping=True)
######################################
# Steady-State Emittance Calculation #
######################################
bunch_intensity = 1e-9 / e # 1C bunch charge
emittance_coupling_factor = 0.5 # for excitation this time
# One can provide specific values for starting emittances,
# but we need to ensure they respect the emittance coupling
# contraint we want to enforce
gemitt_x = 1.1 * tw.eq_gemitt_x # larger horizontal emittance
gemitt_y = emittance_coupling_factor * gemitt_x # enforce the constraint
# One can overwrite sigma_zeta / sigma_delta (larger
# values from potential well distortion for example)
overwrite_sigma_zeta = 1.2 * (tw.eq_gemitt_zeta * tw.bets0) ** 0.5 # larger sigma_zeta
overwrite_sigma_delta = 1.2 * (tw.eq_gemitt_zeta / tw.bets0) ** 0.5 # larger sigma_delta
# A specific time step or relative tolerance for convergence can also be provided.
result = tw.get_ibs_and_synrad_emittance_evolution(
formalism="nagaitsev", # can also be "bjorken-mtingwa"
total_beam_intensity=bunch_intensity,
gemitt_x=gemitt_x, # provided explicitely
gemitt_y=gemitt_y, # provided explicitely
overwrite_sigma_zeta=overwrite_sigma_zeta, # will recompute gemitt_zeta
overwrite_sigma_delta=overwrite_sigma_delta, # will recompute gemitt_zeta
emittance_coupling_factor=emittance_coupling_factor,
emittance_constraint="excitation",
)
# The returned object is a Table
print(result)
# Table: 989 rows, 9 cols
# time gemitt_x nemitt_x gemitt_y nemitt_y gemitt_zeta nemitt_zeta sigma_zeta sigma_delta Kx Ky Kz
# 1.2104374139405318e-06 1.07706e-10 5.26938e-07 5.3853e-11 2.63469e-07 4.05791e-06 0.069754 0.00413633 0.000981042 32.1477 -0.00949839 13.5932
# 9.306647341196751e-05 1.08146e-10 5.29091e-07 5.4073e-11 2.64546e-07 4.06402e-06 0.069859 0.00413944 0.00098178 32.1437 -0.00949715 13.5919
# 0.0001849225094099945 1.08574e-10 5.31185e-07 5.4287e-11 2.65592e-07 4.07004e-06 0.0699624 0.00414251 0.000982507 31.844 -0.00940374 13.4942
# ...
# 0.09057126193146717 1.22535e-10 5.99488e-07 6.12675e-11 2.99744e-07 4.70398e-06 0.0808596 0.00445345 0.00105626 21.8776 -0.00648841 9.10733
# 0.09066311796746519 1.22535e-10 5.99487e-07 6.12675e-11 2.99744e-07 4.70398e-06 0.0808597 0.00445345 0.00105626 21.8776 -0.00648841 9.10732
# 0.09075497400346322 1.22535e-10 5.99487e-07 6.12675e-11 2.99744e-07 4.70399e-06 0.0808598 0.00445345 0.00105626 21.8776 -0.00648841 9.10731
# The results from the table can easily be plotted . One notices
# the transverse coupling factor is respected at all steps.
# Complete source: xfields/examples/005_ibs/006_steady_state_emittances_excitation.py

IBS Kicks for Tracking
When trying to study the interplay of IBS effects with others such as space charge, e-cloud, beamb-beam etc. analytical growth rates are not enough and tracking becomes necessary. In Xfields beam elements are provided to model IBS tracking, which apply momenta kicks to particles. Two kick elements are available:
IBSAnalyticalKick
(based on R. Bruce) for kicks based on analytical growth rates;IBSKineticKick
(based on P. Zenkevich, adapted by M. Zampetakis) for kicks based on diffusion and friction terms from the kinetic theory of gases.
The following example illustrates how to create a kick element, inserting and configuring it for tracking. Refer to the Reference manual for the full list of parameters and their explanation, and to the Physics guide for full information.
See also: xtrack.Line.configure_intrabeam_scattering()
import xfields as xf
import xobjects as xo
import xpart as xp
import xtrack as xt
context = xo.ContextCpu(omp_num_threads="auto")
##########################
# Load xt.Line from file #
##########################
# This is SPS line with proton as particle ref
fname_line_particles = "../../../xtrack/test_data/sps_w_spacecharge/"\
"line_no_spacecharge_and_particle.json"
line: xt.Line = xt.Line.from_json(fname_line_particles)
line.build_tracker(_context=context)
#######################################
# Create and Install IBS Kick Element #
#######################################
# For the analytical kick formalism: kicks are computed based
# on the analytical growth rates (so it needs a formalism)
# ibs_kick = xf.IBSAnalyticalKick(formalism="nagaitsev", num_slices=50)
# For the kinetic formalism: kicks are computed based on the
# friction and diffusion terms of the kinetic theory of gases
ibs_kick = xf.IBSKineticKick(num_slices=50)
# By default the element is off until configuration. Let's install
# the kick at the end of the line and configure it. This internally
# provides the necessary information to the element
line.configure_intrabeam_scattering(
element=ibs_kick, name="ibskick", index=-1, update_every=50
)
############################################
# Define parameters and Generate Particles #
############################################
# Line is for SPS protons at injection
bunch_intensity: int = int(3.5e8)
nemitt_x: float = 2.5e-6
nemitt_y: float = 2.5e-6
sigma_delta: float = 9.56e-4
bunch_length: float = 8.98e-2
particles = xp.generate_matched_gaussian_bunch(
num_particles=10_000,
total_intensity_particles=bunch_intensity,
nemitt_x=nemitt_x,
nemitt_y=nemitt_y,
sigma_z=bunch_length,
line=line,
_context=context,
)
##############################################
# Track now applies an IBS kick at each turn #
##############################################
line.track(particles, num_turns=100, with_progress=5)
# Complete source: xfields/examples/005_ibs/003_tracking_with_kicks.py