A Single Bend with ISR

This test illustrates the effects of incoherent synchrotron radiation (ISR) on a high-energy electron bunch in a bending dipole. The effects of CSR are turned off.

An initially cold (zero emittance), 100 GeV electron bunch with a 0.1 mm rms beam size propagates in an ideal sector bend, through a region of uniform magnetic field.

Due to the stochastic emission of synchrotron radiation, the bunch experiences a mean loss of energy as well as a growth in rms energy spread. Due to the presence of nonzero dispersion, this also results in a growth of the bunch rms emittance.

The resulting quantities are compared against the expressions given in:

N. Yampolsky and B. Carlsten, “Beam debunching due to ISR-induced energy diffusion,” Nucl. Instrum. and Methods in Phys. Res. A, 870, 156-162 (2017), equations (17-20).

In this test, the values of \(\langle{p_x\rangle}\), \(\langle{p_y\rangle}\), \(\langle{p_t\rangle}\), \(\sigma_{p_x}\), \(\sigma_{p_y}\), and \(\sigma_{p_t}\) must agree with nominal values.

In addition, the final values of \(\langle{p_t\rangle}\) and \(\sigma_{p_t}\) must agree with predicted values.

Run

This example can be run either as:

• Python script: python3 run_bend_isr.py or

• ImpactX executable using an input file: impactx input_bend_isr.in

For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.

Python: Script

Listing 187 You can copy this file from examples/incoherent_synchrotron/run_bend_isr.py.

#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-

from impactx import ImpactX, distribution, elements

sim = ImpactX()

# set numerical parameters and IO control
sim.space_charge = False
sim.isr = True
sim.isr_order = 1
sim.slice_step_diagnostics = True

# domain decomposition & space charge mesh
sim.init_grids()

# load a 2 GeV electron beam with an initial
# unnormalized rms emittance of 2 nm
kin_energy_MeV = 100.0e3  # reference energy
bunch_charge_C = 1.0e-12  # used with space charge
npart = 100000  # number of macro particles

#   reference particle
ref = sim.particle_container().ref_particle()
ref.set_charge_qe(-1.0).set_mass_MeV(0.510998950).set_kin_energy_MeV(kin_energy_MeV)

#   particle bunch
distr = distribution.Gaussian(
    lambdaX=1.0e-4,
    lambdaY=1.0e-4,
    lambdaT=1.0e-4,
    lambdaPx=0.0,
    lambdaPy=0.0,
    lambdaPt=0.0,
    muxpx=0.0,
    muypy=0.0,
    mutpt=0.0,
)
sim.add_particles(bunch_charge_C, distr, npart)

# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")

# design the accelerator lattice)
ns = 25  # number of slices per ds in the element

bend = elements.ExactSbend(
    name="sbend_exact", ds=0.1, phi=0.5729577951308232, nslice=ns
)

lattice_line = [monitor, bend, monitor]

# define the lattice
sim.lattice.extend(lattice_line)

# run simulation
sim.track_particles()

# clean shutdown
sim.finalize()

Executable: Input File

Listing 188 You can copy this file from examples/incoherent_synchrotron/input_bend_isr.in.

###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 100000
beam.units = static
beam.kin_energy = 100.0e3  # 100 GeV nominal energy
beam.charge = 1.0e-12  # 1 pC charge
beam.particle = electron
beam.distribution = gaussian
beam.lambdaX = 1.0e-4
beam.lambdaY = beam.lambdaX
beam.lambdaT = 1.0e-4
beam.lambdaPx = 0.0
beam.lambdaPy = 0.0
beam.lambdaPt = 0.0
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0

###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor bend monitor
lattice.nslice = 25

monitor.type = beam_monitor
monitor.backend = h5

bend.type = sbend_exact
bend.ds = 0.1  #used only to parameterize reference arc length
bend.phi = 0.5729577951308232  #deg

#bend.type = sbend
#bend.ds = 0.1  #used only to parameterize reference arc length
#bend.rc = 10.0  #0.5729577951308232 deg

###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
algo.isr = true
algo.isr_order = 1

###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
diag.eigenemittances = true

Analyze

We run the following script to analyze correctness:

Script analysis_bend_isr.py

Listing 189 You can copy this file from examples/incoherent_synchrotron/analysis_bend_isr.py.

#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#

import glob
import re

import numpy as np
import pandas as pd


def read_file(file_pattern):
    for filename in glob.glob(file_pattern):
        df = pd.read_csv(filename, delimiter=r"\s+")
        if "step" not in df.columns:
            step = int(re.findall(r"[0-9]+", filename)[0])
            df["step"] = step
        yield df


def read_time_series(file_pattern):
    """Read in all CSV files from each MPI rank (and potentially OpenMP
    thread). Concatenate into one Pandas dataframe.

    Returns
    -------
    pandas.DataFrame
    """
    return pd.concat(
        read_file(file_pattern),
        axis=0,
        ignore_index=True,
    )  # .set_index('id')


# read reduced diagnostics
rbc = read_time_series("diags/reduced_beam_characteristics.*")

s = rbc["s"]
px_mean = rbc["mean_px"]
py_mean = rbc["mean_py"]
pt_mean = rbc["mean_pt"]
sigma_px = rbc["sigma_px"]
sigma_py = rbc["sigma_py"]
sigma_pt = rbc["sigma_pt"]

px_meani = px_mean.iloc[0]
py_meani = py_mean.iloc[0]
pt_meani = pt_mean.iloc[0]
sig_pxi = sigma_px.iloc[0]
sig_pyi = sigma_py.iloc[0]
sig_pti = sigma_pt.iloc[0]

length = len(s) - 1

sf = s.iloc[length]

px_meanf = px_mean.iloc[length]
py_meanf = py_mean.iloc[length]
pt_meanf = pt_mean.iloc[length]
sig_pxf = sigma_px.iloc[length]
sig_pyf = sigma_py.iloc[length]
sig_ptf = sigma_pt.iloc[length]

print("Initial Beam:")
print(f"  px_mean={px_meani:e} py_mean={py_meani:e} pt_mean={pt_meani:e}")
print(f"  sig_px={sig_pxi:e} sig_py={sig_pyi:e} sig_pt={sig_pti:e}")


atol = 1.0e-6
print(f"  atol={atol}")
assert np.allclose(
    [px_meani, py_meani, pt_meani, sig_pxi, sig_pyi, sig_pti],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    atol=atol,
)

# Physical constants:
re_classical = 2.8179403205e-15  # classical electron radius
lambda_compton_reduced = 3.8615926744e-13  # reduced Compton wavelength

# Problem parameters:
ds = 0.1
rc = 10.0
gamma = 195696.117901
num_particles = 10000

# Characteristic length of energy loss:
length_isr = 3.0 / 2.0 * rc**2 / (re_classical * gamma**3)

print("")
print("Length and characteristic length for energy loss [m]:")
print(f" Length={ds:e} Length_ISR={length_isr:e}")

# Predicted energy loss and energy spread:
dpt = 2.0 / 3.0 * re_classical * ds / (rc**2) * gamma**3

dsigpt2 = (
    55.0
    / (24 * np.sqrt(3.0))
    * lambda_compton_reduced
    * re_classical
    * ds
    / rc**3
    * gamma**5
)
dsigpt = np.sqrt(dsigpt2)

print("")
print("Final Beam:")
print(f" pt_mean={pt_meanf:e} sig_pt={sig_ptf}")
print("Predicted (assuming that Length << Length_isr):")
print(f" pt_mean={dpt:e} sig_pt={dsigpt:e}")
print("")

rtol = 10.0 * num_particles**-0.5  # from random sampling of a smooth distribution
assert np.allclose(
    [pt_meanf, sig_ptf],
    [
        dpt,
        dsigpt,
    ],
    rtol=rtol,
    atol=atol,
)

A Single Bend with ISR, Reference Energy Loss

This is identical to the preceding test, except for the flag isr_on_ref_part option being enabled. In this test, the reference particle experiences radiative energy loss. For the beam particles, whose coordinates and momenta are measured relative to the reference particle, the primary effect of ISR is to induce an increase in energy spread. Little effect is visible on the beam centroid.

Run

This example can be run either as:

• Python script: python3 run_bend_isr_ref.py or

• ImpactX executable using an input file: impactx input_bend_isr_ref.in

For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.

Python: Script

Listing 190 You can copy this file from examples/incoherent_synchrotron/run_bend_isr_ref.py.

#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-

from impactx import ImpactX, distribution, elements

sim = ImpactX()

# set numerical parameters and IO control
sim.space_charge = False
sim.isr = True
sim.isr_order = 1
sim.isr_on_ref_part = True
sim.slice_step_diagnostics = True

# domain decomposition & space charge mesh
sim.init_grids()

# load a 2 GeV electron beam with an initial
# unnormalized rms emittance of 2 nm
kin_energy_MeV = 100.0e3  # reference energy
bunch_charge_C = 1.0e-12  # used with space charge
npart = 100000  # number of macro particles

#   reference particle
ref = sim.particle_container().ref_particle()
ref.set_charge_qe(-1.0).set_mass_MeV(0.510998950).set_kin_energy_MeV(kin_energy_MeV)

#   particle bunch
distr = distribution.Gaussian(
    lambdaX=1.0e-4,
    lambdaY=1.0e-4,
    lambdaT=1.0e-4,
    lambdaPx=0.0,
    lambdaPy=0.0,
    lambdaPt=0.0,
    muxpx=0.0,
    muypy=0.0,
    mutpt=0.0,
)
sim.add_particles(bunch_charge_C, distr, npart)

# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")

# design the accelerator lattice)
ns = 25  # number of slices per ds in the element

bend = elements.ExactSbend(
    name="sbend_exact", ds=0.1, phi=0.5729577951308232, nslice=ns
)

lattice_line = [monitor, bend, monitor]

# define the lattice
sim.lattice.extend(lattice_line)

# run simulation
sim.track_particles()

# clean shutdown
sim.finalize()

Executable: Input File

Listing 191 You can copy this file from examples/incoherent_synchrotron/input_bend_isr_ref.in.

###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 100000
beam.units = static
beam.kin_energy = 100.0e3  # 100 GeV nominal energy
beam.charge = 1.0e-12  # 1 pC charge
beam.particle = electron
beam.distribution = gaussian
beam.lambdaX = 1.0e-4
beam.lambdaY = beam.lambdaX
beam.lambdaT = 1.0e-4
beam.lambdaPx = 0.0
beam.lambdaPy = 0.0
beam.lambdaPt = 0.0
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0

###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor bend monitor
lattice.nslice = 25

monitor.type = beam_monitor
monitor.backend = h5

bend.type = sbend_exact
bend.ds = 0.1  #used only to parameterize reference arc length
bend.phi = 0.5729577951308232  #deg

#bend.type = sbend
#bend.ds = 0.1  #used only to parameterize reference arc length
#bend.rc = 10.0  #0.5729577951308232 deg

###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false
algo.isr = true
algo.isr_order = 1
algo.isr_on_ref_part = true

###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
diag.eigenemittances = true

Analyze

We run the following script to analyze correctness:

Script analysis_bend_isr_ref.py

Listing 192 You can copy this file from examples/incoherent_synchrotron/analysis_bend_isr_ref.py.

#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#

import glob
import re

import numpy as np
import pandas as pd


def read_file(file_pattern):
    for filename in glob.glob(file_pattern):
        df = pd.read_csv(filename, delimiter=r"\s+")
        if "step" not in df.columns:
            step = int(re.findall(r"[0-9]+", filename)[0])
            df["step"] = step
        yield df


def read_time_series(file_pattern):
    """Read in all CSV files from each MPI rank (and potentially OpenMP
    thread). Concatenate into one Pandas dataframe.

    Returns
    -------
    pandas.DataFrame
    """
    return pd.concat(
        read_file(file_pattern),
        axis=0,
        ignore_index=True,
    )  # .set_index('id')


# read reduced diagnostics
rbc = read_time_series("diags/reduced_beam_characteristics.*")
ref = read_time_series("diags/ref_particle.*")

s = rbc["s"]
px_mean = rbc["mean_px"]
py_mean = rbc["mean_py"]
pt_mean = rbc["mean_pt"]
sigma_px = rbc["sigma_px"]
sigma_py = rbc["sigma_py"]
sigma_pt = rbc["sigma_pt"]
ref_gamma = ref["gamma"]

px_meani = px_mean.iloc[0]
py_meani = py_mean.iloc[0]
pt_meani = pt_mean.iloc[0]
sig_pxi = sigma_px.iloc[0]
sig_pyi = sigma_py.iloc[0]
sig_pti = sigma_pt.iloc[0]
ref_gammai = ref_gamma.iloc[0]

length = len(s) - 1

sf = s.iloc[length]

px_meanf = px_mean.iloc[length]
py_meanf = py_mean.iloc[length]
pt_meanf = pt_mean.iloc[length]
sig_pxf = sigma_px.iloc[length]
sig_pyf = sigma_py.iloc[length]
sig_ptf = sigma_pt.iloc[length]
ref_gammaf = ref_gamma.iloc[-1]

print("Initial Beam:")
print(f"  px_mean={px_meani:e} py_mean={py_meani:e} pt_mean={pt_meani:e}")
print(f"  sig_px={sig_pxi:e} sig_py={sig_pyi:e} sig_pt={sig_pti:e}")
print(f" reference gamma={ref_gammai:e}")

atol = 1.0e-6
print(f"  atol={atol}")
assert np.allclose(
    [px_meani, py_meani, pt_meani, sig_pxi, sig_pyi, sig_pti],
    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    atol=atol,
)

# Physical constants:
re_classical = 2.8179403205e-15  # classical electron radius
lambda_compton_reduced = 3.8615926744e-13  # reduced Compton wavelength

# Problem parameters:
ds = 0.1
rc = 10.0
gamma = 195696.117901
num_particles = 10000

# Characteristic length of energy loss:
length_isr = 3.0 / 2.0 * rc**2 / (re_classical * gamma**3)

print("")
print("Length and characteristic length for energy loss [m]:")
print(f" Length={ds:e} Length_ISR={length_isr:e}")

# Predicted energy loss and energy spread:
dpt = 2.0 / 3.0 * re_classical * ds / (rc**2) * gamma**3

dsigpt2 = (
    55.0
    / (24 * np.sqrt(3.0))
    * lambda_compton_reduced
    * re_classical
    * ds
    / rc**3
    * gamma**5
)
dsigpt = np.sqrt(dsigpt2)

# Relative (total) energy change in the reference particle:
dpt_ref = -(ref_gammaf - ref_gammai) / ref_gammai

print("")
print("Final Beam:")
print(f" pt_mean={pt_meanf:e} sig_pt={sig_ptf}")
print("Predicted (assuming that Length << Length_isr):")
print(f" pt_mean={dpt:e} sig_pt={dsigpt:e}")
print(f" reference gamma={ref_gammaf:e}")
print("")

print(f" relative reference energy loss={dpt_ref:e}")

rtol = 10.0 * num_particles**-0.5  # from random sampling of a smooth distribution
assert np.allclose(
    [dpt_ref, sig_ptf],
    [
        dpt,
        dsigpt,
    ],
    rtol=rtol,
    atol=atol,
)