Ballistic Compression Using a Short RF Element

A 20 MeV electron beam propagates through a short RF element near zero-crossing, inducing a head-tail energy correlation. This is followed by ballistic motion in a drift, which is used to compress the rms bunch length from 16 ps to 10 ps (compression of 5/3).

The beam is not exactly on-crest (phase = -89.5 deg), so there is an energy gain of 4.5 MeV.

The transverse emittance is sufficiently small that the horizontal and verticle beam size are essentially unchanged. Due to RF curvature, there is some growth of the longitudinal emittance.

In this test, the initial and final values of \(\sigma_x\), \(\sigma_y\), \(\sigma_t\), \(\epsilon_x\), \(\epsilon_y\), and \(\epsilon_t\) must agree with nominal values.

Run

This example can be run either as:

  • Python script: python3 run_compression.py or

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

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

Listing 42 You can copy this file from examples/compression/run_compression.py.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Marco Garten, 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.particle_shape = 2  # B-spline order
sim.space_charge = False
# sim.diagnostics = False  # benchmarking
sim.slice_step_diagnostics = True

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

# load a 250 MeV proton beam with an initial
# unnormalized rms emittance of 1 mm-mrad in all
# three phase planes
kin_energy_MeV = 20.0  # reference energy
bunch_charge_C = 1.0e-9  # used with space charge
npart = 10000  # 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.Waterbag(
    lambdaX=0.5e-3,
    lambdaY=0.5e-3,
    lambdaT=5.0e-3,
    lambdaPx=1.0e-5,
    lambdaPy=1.0e-5,
    lambdaPt=4.0e-6,
)
sim.add_particles(bunch_charge_C, distr, npart)

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

# design the accelerator lattice
sim.lattice.append(monitor)
#   Short RF cavity element
shortrf1 = elements.ShortRF(V=1000.0, freq=1.3e9, phase=-89.5)
#   Drift element
drift1 = elements.Drift(ds=1.7)

sim.lattice.extend([shortrf1, drift1])

sim.lattice.append(monitor)

# run simulation
sim.evolve()

# clean shutdown
sim.finalize()
Listing 43 You can copy this file from examples/compression/input_compression.in.
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 20.0
beam.charge = 1.0e-9
beam.particle = electron
beam.distribution = waterbag
beam.lambdaX = 0.5e-3
beam.lambdaY = 0.5e-3
beam.lambdaT = 5.0e-3
beam.lambdaPx = 1.0e-5
beam.lambdaPy = 1.0e-5
beam.lambdaPt = 4.0e-6
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0


###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor shortrf1 drift1 monitor

monitor.type = beam_monitor
monitor.backend = h5

shortrf1.type = shortrf
shortrf1.V = 1000.0
shortrf1.freq = 1.3e9
shortrf1.phase = -89.5

drift1.type = drift
drift1.ds = 1.7


###############################################################################
# Algorithms
###############################################################################
algo.particle_shape = 2
algo.space_charge = false

Analyze

We run the following script to analyze correctness:

Script analysis_compression.py
Listing 44 You can copy this file from examples/compression/analysis_compression.py.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#

import numpy as np
import openpmd_api as io
from scipy.stats import moment


def get_moments(beam):
    """Calculate standard deviations of beam position & momenta
    and emittance values

    Returns
    -------
    sigx, sigy, sigt, emittance_x, emittance_y, emittance_t
    """
    sigx = moment(beam["position_x"], moment=2) ** 0.5  # variance -> std dev.
    sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
    sigy = moment(beam["position_y"], moment=2) ** 0.5
    sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
    sigt = moment(beam["position_t"], moment=2) ** 0.5
    sigpt = moment(beam["momentum_t"], moment=2) ** 0.5

    epstrms = beam.cov(ddof=0)
    emittance_x = (sigx**2 * sigpx**2 - epstrms["position_x"]["momentum_x"] ** 2) ** 0.5
    emittance_y = (sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2) ** 0.5
    emittance_t = (sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2) ** 0.5

    return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)


# openPMD data series at the beam monitors
series = io.Series("diags/openPMD/monitor.h5", io.Access.read_only)

# first and last step
final_step = list(series.iterations)[-1]
first_it = series.iterations[1]
final_it = series.iterations[final_step]

# initial beam & reference particle gamma
initial = first_it.particles["beam"].to_df()
initial_gamma_ref = first_it.particles["beam"].get_attribute("gamma_ref")

# final beam & reference particle gamma
final = final_it.particles["beam"].to_df()
final_gamma_ref = final_it.particles["beam"].get_attribute("gamma_ref")

# compare number of particles
num_particles = 10000
assert num_particles == len(initial)
assert num_particles == len(final)

print("Initial Beam:")
sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(initial)
print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
print(
    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
)
print(f"  gamma={initial_gamma_ref:e}")

atol = 0.0  # ignored
rtol = 1.9 * num_particles**-0.5  # from random sampling of a smooth distribution
print(f"  rtol={rtol} (ignored: atol~={atol})")

assert np.allclose(
    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t, initial_gamma_ref],
    [
        5.0e-04,
        5.0e-04,
        5.0e-03,
        4.952764e-09,
        5.028325e-09,
        1.997821e-08,
        40.1389432485322889,
    ],
    rtol=rtol,
    atol=atol,
)


print("")
print("Final Beam:")
sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(final)
print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
print(
    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
)
print(f"  gamma={final_gamma_ref:e}")


atol = 0.0  # ignored
rtol = 1.9 * num_particles**-0.5  # from random sampling of a smooth distribution
print(f"  rtol={rtol} (ignored: atol~={atol})")

assert np.allclose(
    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t, final_gamma_ref],
    [
        5.004995e-04,
        5.005865e-04,
        3.033949e-03,
        4.067876e-09,
        4.129937e-09,
        6.432081e-05,
        48.8654787469061860,
    ],
    rtol=rtol,
    atol=atol,
)