Numerical Methods in Accelerator Physics

Lecture Series by Dr. Adrian Oeftiger

Lecture 6

Run this notebook online!

Interact and run this jupyter notebook online:

via the public mybinder.org service:
mybinder.org logo
via the public gesis.org service:
gesis.org logo

Also find this lecture rendered as HTML slides on github $\nearrow$ along with the source repository $\nearrow$.

Run this first!

Imports and modules:

In [1]:
from config import (np, plt, plot_rfwave, tqdm, trange, 
                    beta, gamma, Machine, track_one_turn,
                    charge, mass, emittance, hamiltonian, 
                    plot_hamiltonian, plot_rf_overview,
                    plot_dist, plot_mp)
from scipy.constants import m_p, e, c
%matplotlib inline

If the progress bar by tqdm (trange) in this document does not work, run this:

In [2]:
!jupyter nbextension enable --py widgetsnbextension

Enabling notebook extension jupyter-js-widgets/extension...
      - Validating: OK

Refresher!

  • Hamiltonian for synchrotron motion (rf bucket, stationary case and small-amplitude approximation)
  • Monte-Carlo Techniques: sampling of distribution functions
  • uniformly distributed pseudo-random numbers (with Linear Congruential Generators)
  • Gaussian distributed pseudo-random numbers (Box-Muller method)
  • Equilibrium distributions
  • Equilibrium at small amplitudes: harmonic oscillation $\leftrightarrow$ Gaussian distribution
  • Rejection sampling

Today!

  1. Simulation of Full CERN PS Acceleration Ramp
  2. Initialisation of Phase Space (Cont'd)
  3. Emittance Preservation & Injection Errors

Part I: Simulation of Full CERN PS Acceleration Ramp

CERN PS Parameters (once more)

For reference, the CERN Proton Synchrotron (PS) parameters:

  • has a circumference of 2π·100m
  • takes protons from the PS Booster at a kinetic energy of 2GeV corresponding to $\gamma=3.13$
  • injects with 50kV of rf voltage, up to 200kV for ramp
  • runs at harmonic $h=7$
  • has a momentum compaction factor of $\alpha_c=0.027$
  • typical acceleration rate of (up to) $\dot{B}=2$ T/s, the bending radius is $\rho=70.08$ m

Exercise

In the following, let us properly simulate the full acceleration ramp (including crossing transition) with a bunch of small emittance (phase-space area)!

You as operator in the (simulation) control centre shall pay attention to preserve the longitudinal emittance throughout the acceleration!

$\implies$ Investigate the tracking results. Can you explain what happens? Think about what you need to change in the machine settings to achieve emittance preservation.

We start by instantiating the PS, Machine(...):

(hint: hit both Shift+Tab keys inside the parentheses () below to get info about the possible arguments to Machine)

In [3]:
m = Machine()

assert m.phi_s > 0, "machine is not accelerating...?"

If you want to check the rf systems setup, uncomment (remove #) the following plot command from the previous lecture:

In [4]:
# plot_rf_overview(m);

We initialise a Gaussian bunch distribution with very small rms bunch length $\sigma_z=1\,$m (so that the small-amplitude approximation holds):

In [5]:
sigma_z = 1

The "matched" rms momentum spread $\sigma_{\Delta p}$ (remember, $\sigma_{\Delta p}$ and $\sigma_z$ are linked via equal Hamiltonian values $\mathcal{H}_0$, the equilibrium condition):

In [6]:
sigma_deltap = np.sqrt(
    2 * m.p0() / -m.eta(0) * 
    charge * m.voltage * np.pi * m.harmonic / (beta(gamma(m.p0())) * c * m.circumference**2)
) * sigma_z

Generating Macro-particles via Box-Muller

(NB: with this small emittance / rms size in phase space, the Gaussian tails cannot reach outside the separatrix due to numerical reasons, so no rejection sampling is needed in this particular case.

Finite machine precision at FP64 is $\varepsilon\approx 2^{-53}$, therefore one can only generate pseudo-random numbers from the Gaussian distribution up to an amplitude of $x$ where $\exp(-x^2/2)=2^{-53}$, i.e.

In [7]:
np.sqrt(2 * -np.log(2**-53))
Out[7]:
8.571674348652905

$\implies$ given $\sigma_z=1\,$m, no particles can be generated outside of $z=8.6\,$m and the equivalent Hamiltonian contour in phase space)

In [8]:
N = 1000
In [9]:
np.random.seed(12345)

z = np.random.normal(loc=0, scale=sigma_z, size=N)
deltap = np.random.normal(loc=0, scale=sigma_deltap, size=N)
In [10]:
plot_hamiltonian(m)
plt.scatter(z, deltap / m.p0(), marker='.', s=1);

Compute Duration of Ramp

$(\Delta\gamma)_\mathrm{turn} = \cfrac{\Delta E_\mathrm{tot}}{m_0c^2} = \cfrac{qV\sin(\varphi_s)}{m_0c^2}$

such that accelerating from $\gamma_\mathrm{ref}=3.1$ to $\gamma_\mathrm{ref}=27.7$ takes as many turns as

$n_\mathrm{turns}=\cfrac{27.7-3.1}{(\Delta\gamma)_\mathrm{turn}}$

In [11]:
dgamma_per_turn = charge * m.voltage * np.sin(m.phi_s) / (mass * c**2)
n_turns = int(np.ceil((27.7 - 3.13) / dgamma_per_turn))
n_turns
Out[11]:
261756

Record longitudinal emittance during tracking:

In [12]:
epsn_z = np.zeros(n_turns, dtype=np.float64)

epsn_z[0] = emittance(z, deltap)

Let's go tracking! ($\approx 1\,$min total time)

In [13]:
for i_turn in trange(1, n_turns):
    z, deltap = track_one_turn(z, deltap, m)
    
    epsn_z[i_turn] = emittance(z, deltap)

  0%|          | 0/261755 [00:00<?, ?it/s]

We have reached the extraction energy of $\gamma=27.7$:

In [14]:
m.gamma_ref
Out[14]:
27.69995958340248

How did we do?

Check that we properly conserved the rms emittance:

In [15]:
plt.plot(np.arange(n_turns) / 100000, epsn_z / e)

plt.xlabel('Turns [100k]')
plt.ylabel('$\epsilon_z$ [eV.s]');

$\implies$ use the following phase-space plot besides the plotting functions plot_hamiltonian and plot_rf_overview above (from last lecture) as diagnostics. You can stop the tracking after a certain turn to investigate, e.g. by changing the value of n_turns inside the trange counter in the for loop.

In [16]:
plot_hamiltonian(m);
plt.scatter(z % 600, deltap / m.p0(), marker='.', s=1);

Part II: Initialisation of Phase Space (Cont'd)

Matching with Nonlinear Hamiltonians

Last lecture we looked at the thermal equilibrium distribution in small-amplitude approximation ($=$ Gaussian distribution) to generate a bunch of macro-particles in longitudinal phase space.

For the Gaussian distribution, the target rms bunch length $\hat{\sigma}_z$ directly enters the generator as a multiplication factor: a random variable $x$ following a Gaussian normal distribution has $\sigma_x=1$, then $z=\hat{\sigma}_z\cdot x$ will provide the Gaussian normal distributed random variable $z$ with $\sigma_z=\hat{\sigma}_z$.

For a general thermal equilibrium distribution in a nonlinear Hamiltonian, the PDF becomes

$$\psi=\psi(\mathcal{H})\propto \exp\left(\frac{\mathcal{H}}{\mathcal{H}_0}\right)$$

for a constant parameter $\mathcal{H}_0$. The value of $\mathcal{H}_0$ is determined through $\hat{\sigma}_z$ such that $\psi$ yields a distribution with corresponding rms bunch length. For nonlinear $\mathcal{H}$, the parameter $\mathcal{H}_0$ will generally have to be computed numerically. This process is referred to as "matching"!

PyHEADTAIL

In the following we will investigate the full nonlinear Hamiltonian case, $\psi=\psi(\mathcal{H})$.

Let's work with the python library PyHEADTAIL (simulation suite for beam dynamics). If PyHEADTAIL is not yet installed, please pip-install it:

$ pip install PyHEADTAIL==1.16.1

RF Buckets in PyHEADTAIL

PyHEADTAIL provides a class to represent rf buckets (for plotting as well as for matching):

In [17]:
from PyHEADTAIL.trackers.rf_bucket import RFBucket

PyHEADTAIL v1.16.1

We define a convenience function to provide RFBucket instances given a Machine instance:

In [18]:
def get_pyht_rfbucket(machine):
    m = machine
    deltap_per_turn = charge * m.voltage / (beta(gamma(m.p0())) * c) * np.sin(m.phi_s)
    rfb = RFBucket(m.circumference, m.gamma_ref, mass, charge, [m.alpha_c], deltap_per_turn, 
                   [m.harmonic], [m.voltage], [np.pi + m.phi_s])
    # PyHEADTAIL has a different convention for the phase and is offset by π compared to our lecture
    return rfb

Visualising the Distributions in the RF Bucket

In [19]:
m = Machine(gamma_ref=3.13)
In [20]:
rfb = get_pyht_rfbucket(m)
In [21]:
from PyHEADTAIL.particles.rfbucket_matching import (
    ThermalDistribution, WaterbagDistribution, ParabolicDistribution)
In [22]:
sigma_z = 8

Computing the (initial) guess for $\mathcal{H}_0$ based on the small-amplitude approximation:

In [23]:
H0 = rfb.guess_H0(sigma_z, from_variable='sigma')
H0
Out[23]:
153.62263838682642

Reminder: Equilibrium at Small Amplitudes

Small-amplitude approximation and a stationary rf bucket below transition, $\varphi_s=0$: particles follow harmonic oscillation with

$$\mathcal{H}_\mathrm{stat,small}(z,\Delta p) = \frac{1}{2}\frac{-\eta}{p_0} \Delta p{}^2 + \frac{qV}{\beta c}\cdot \cfrac{\pi h}{C^2}\cdot z^2 $$

The PDF of the thermal distribution becomes

$$\psi(\mathcal{H})\propto\exp\left(\frac{\mathcal{H}_{stat,small}}{\mathcal{H}_0}\right)=\exp\left(\frac{1}{2}\frac{\frac{qV\,2\pi h}{\beta c C^2}}{\mathcal{H}_0}\cdot z^2\right) \cdot \exp\left(\frac{1}{2}\frac{-\eta/p_0}{\mathcal{H}_0}\cdot \Delta p^2\right)$$

which is simply a bi-Gaussian distribution in $\Delta p$ and $z$!

Computing $\mathcal{H}_0$

$$\psi(\mathcal{H})\propto\exp\left(\frac{\mathcal{H}_{stat,small}}{\mathcal{H}_0}\right)=\exp\left(\frac{1}{2}\frac{\frac{qV\,2\pi h}{\beta c C^2}}{\mathcal{H}_0}\cdot z^2\right) \cdot \exp\left(\frac{1}{2}\frac{-\eta/p_0}{\mathcal{H}_0}\cdot \Delta p^2\right)$$

The constant $\mathcal{H}_0$ is determined by a choice of the rms bunch length $\sigma_z$ (or equivalently the rms momentum deviation $\sigma_\delta=\sigma_{\Delta p}/p_0$):

$$\mathcal{H}_0 = \mathcal{H}_\mathrm{stat,small}(\sigma_z,\Delta p=0) = \frac{qV}{\beta c}\cdot \cfrac{\pi h}{C^2}\cdot \sigma_z{}^2 \mathrel{\color{red}{\stackrel{!}{=}}} \mathcal{H}_\mathrm{stat,small}(z=0,\sigma_{\Delta p}) = \frac{1}{2}\frac{-\eta}{p_0} \sigma_{\Delta p}{}^2$$

$\implies$ as equilibrium or matching condition, $\sigma_z$ and $\sigma_\delta$ are linked to each other via equal Hamiltonian values!

In [24]:
plot_dist(ThermalDistribution, rfb, H0=0.05*H0);

$\implies$ adjust the $\sigma_z$ and observe the impact of the non-linearities on the matched $\psi$ (small $\sigma_z$ will make the distribution more Gaussian-like)!

Change the distribution type from ThermalDistribution to WaterbagDistribution (or e.g. ParabolicDistribution):

In [25]:
plot_dist(WaterbagDistribution, rfb, H0=H0*2);

Matching Algorithm

Given a target rms bunch length $\hat{\sigma}_z$, define the following "matching" algorithm:

  1. initial value for $\mathcal{H}_0$ given by small-amplitude approximation (the bi-Gaussian case)
  2. compute (numerically) the resulting $\sigma_z$ for full $\psi(\mathcal{H}/\mathcal{H}_0)$
  3. compare $\sigma_z$ with target $\hat{\sigma}_z$
  4. adjust $\mathcal{H}_0$ (for $\sigma_z<\hat{\sigma}_z$ increase $\mathcal{H}_0$ and vice versa), go to step 2 until converged

This matching algorithm is implemented in the RFBucketMatcher class in PyHEADTAIL:

In [26]:
from PyHEADTAIL.particles.generators import RFBucketMatcher
In [27]:
rfb_matcher = RFBucketMatcher(rfb, ThermalDistribution, sigma_z=sigma_z)
rfb_matcher.integrationmethod = 'cumtrapz' # better behaved numerical integration method

Calling the RFBucketMatcher.generate method will

(1.) iterate on $\mathcal{H}_0$ until the numerical integration of $\psi(\mathcal{H})$ for the bunch length converges to $\hat{\sigma}_z$, and then

(2.) sample this $\psi(\mathcal{H})$ by rejection sampling (see previous lecture) to generate the macro-particle phase-space coordinates $(z,\delta)$

In [28]:
z, delta, _, _ = rfb_matcher.generate(int(1e5))

*** Maximum RMS bunch length 8.888845463360669m.
... distance to target bunch length: -7.9156e+00
... distance to target bunch length: 6.5599e-01
... distance to target bunch length: 6.1413e-01
... distance to target bunch length: -3.2377e-01
... distance to target bunch length: 2.5291e-01
... distance to target bunch length: 6.5279e-02
... distance to target bunch length: -4.6398e-03
... distance to target bunch length: 2.6859e-04
... distance to target bunch length: 1.0364e-06
... distance to target bunch length: -8.7361e-06
--> Bunch length: 8.000001036447514
--> Emittance: 2.108543434891298

The relevant code in the RFBucketMatcher implementing the rejection sampling:

xmin, xmax = self.rfbucket.z_left, self.rfbucket.z_right
ymin = -self.rfbucket.dp_max(self.rfbucket.z_right)
ymax = -ymin

uniform = np.random.uniform
n_gen = macroparticlenumber
u = uniform(low=xmin, high=xmax, size=n_gen)
v = uniform(low=ymin, high=ymax, size=n_gen)
s = uniform(size=n_gen)

def mask_out(s, u, v):
    return s >= self.psi(u, v)

masked_out = mask_out(s, u, v)
while np.any(masked_out):
    masked_ids = np.where(masked_out)[0]
    n_gen = len(masked_ids)
    u[masked_out] = uniform(low=xmin, high=xmax, size=n_gen)
    v[masked_out] = uniform(low=ymin, high=ymax, size=n_gen)
    s[masked_out] = uniform(size=n_gen)
    masked_out[masked_ids] = mask_out(
        s[masked_out], u[masked_out], v[masked_out]
    )

Compare the converged value of $\mathcal{H}_0$...

In [29]:
rfb_matcher.psi_object.H0
Out[29]:
216.61619727806732

... to the small-amplitude approximation value:

In [30]:
rfb.guess_H0(0.1 * sigma_z, from_variable='sigma')
Out[30]:
1.5362263838682642

$\implies$ are the matched $\mathcal{H}_0$ and the small-amplitude approximation value closer for smaller target $\hat{\sigma}_z$?

Let's have a look at the generated macro-particle distribution:

In [31]:
plot_mp(z, delta, rfb);

Does the rms bunch length of the macro-particle distribution match the chosen target $\hat{\sigma}_z$?

In [32]:
np.std(z)
Out[32]:
8.00462076948197
In [33]:
sigma_z
Out[33]:
8

Matching Algorithm in More Detail: Root-Finding

Equilibrium distributions for nonlinear Hamiltonians

$\implies$ thermal distribution $\psi\propto \exp(\mathcal{H}/\mathcal{H}_0)$ where constant $\mathcal{H}_0$ is determined by rms bunch length (or momentum spread)

We discussed:

(1.) iterate on $\mathcal{H}_0$ until the numerical integration of $\psi(\mathcal{H})$ for the bunch length converges to $\hat{\sigma}_z$

(...)

$\implies$ this poses a root-finding problem!

Moments of Distribution $\psi$

The relevant moments of a generic distribution function $\psi$ are

  • monopole:
$$Q_z = \int\limits_{z_\text{left}}^{z_\text{right}} dz \int_\text{separatrix} d(\Delta p)\cdot\psi(z,\Delta p)$$
  • dipole (or centroid):
$$\langle z\rangle = \cfrac{1}{Q_z} \int\limits_{z_\text{left}}^{z_\text{right}} dz \int_\text{separatrix} d(\Delta p)\cdot\psi(z,\Delta p) \cdot z$$
  • quadrupole (or squared rms bunch length):
$$\sigma_z^2 = \cfrac{1}{Q_z}\int\limits_{z_\text{left}}^{z_\text{right}} dz \int_\text{separatrix} d(\Delta p)\cdot\psi(z,\Delta p) \cdot \left[z - \langle z \rangle \right]^2$$

where $z_\text{left}$ and $z_\text{right}$ denote the left and right boundary of the rf bucket, and the second integral integrates over $\Delta p$ from the lower to the upper separatrix at position $z$.

The Root-finding Problem

The equation for the rms bunch length expresses the root-finding problem.

Given a (non-normalised) thermal distribution $\psi(\mathcal{H}) = \exp(\mathcal{H}/\mathcal{H}_0)$ and a target rms bunch length $\hat{\sigma}_z$, then $\mathcal{H}_0$ is determined by the integral equation

$$0 = \hat{\sigma}_z{}^2 - \cfrac{1}{Q_z}\int\limits_{z_\text{left}}^{z_\text{right}} dz \int_\text{separatrix} d(\Delta p)\cdot\exp\left(\cfrac{\mathcal{H}(z,\Delta p)}{\color{red}{\mathcal{H}_0}}\right) \cdot \left[z - \langle z \rangle \right]^2$$

In the case of a linear harmonic oscillator Hamiltonian $\mathcal{H}_\text{lin}$ (typical result of small-amplitude approximation), the solution is given by

$$\mathcal{H}_{0,\text{lin}}(\hat{\sigma}_z)=\frac{qV}{\beta c}\cdot \frac{\pi h}{C^2}\cdot \hat{\sigma}_z{}^2$$

For nonlinear $\mathcal{H}(z,\Delta p)$, numerical root-finding algorithms can be employed (the integrals are solved numerically, too). Using as input the expression $\mathcal{H}_{0,\text{lin}}(\sigma_z)$, one varies $\sigma_z$ and obtains a converged value via

  • Brent-Dekker algorithm: interval for $\sigma_z$ with sign change of function, choose e.g. $[\hat{\sigma}_z / 100, C]$ for circumference $C$, or
  • secant method: starting point for $\sigma_z$, choose $\hat{\sigma}_z$ (and e.g. $0.9\hat{\sigma}_z$)

Let's Use scipy.optimize!¶

In [34]:
from scipy.optimize import brentq, newton

At what $x$ does $\cfrac{1}{\sqrt{x}+1}$ take the value 0.4?

In [35]:
def f(x):
    return 1 / (np.sqrt(x) + 1) - 0.4

Brent-Dekker algorithm with interval $x\in[0, 4]$:

In [36]:
brentq(f, 0, 4)
Out[36]:
2.250000000000207

... and Newton's Secant Method

Secant method with initial values $x_0=0$, $x_1=10^{-4}$:

In [37]:
newton(f, 0)
Out[37]:
2.2499999999999956

How does the function look like?

In [38]:
x = np.linspace(0, 4, 1000)
In [39]:
plt.plot(x, f(x))
plt.axhline(0, c='k', lw=2)
plt.xlabel('$x$')
plt.ylabel('$f(x)$');

Implementing the Secant Method

In [40]:
def secant_method(f, x0, x1, iterations, rtol=2e-12):
    """Return the root calculated using the secant method."""
    for i in range(iterations):
        x2 = x1 - f(x1) * (x1 - x0) / float(f(x1) - f(x0))
        if np.abs(x2 - x1) / x1 < rtol:
            break
        x0, x1 = x1, x2
    return x2
In [41]:
secant_method(f, 0, 1, 100)
Out[41]:
2.249999999999999

(The scipy implementation of the secant method in newton() computes the second guess $x_1$ as $x_1=x_0\cdot 1.0001\pm 0.0001$.)

Part III: Emittance Preservation & Injection Errors

Sources of Simulated Emittance Growth

Numerical:

  • distribution is not an exact equilibrium distribution (often the case for Gaussians)
  • bad (non-symplectic) numerical integrator

Physical:

  • magnet power supply noise, resonances
  • multi-particle / collective effects (intra-beam scattering, instabilities)
  • mismatch at injection into synchrotron

Section A: Numerical Emittance Growth

Based on the stationary rf bucket (no acceleration) case, we investigate the difference between initialising a bi-Gaussian and a matched thermal distribution based on the nonlinear Hamiltonian $\mathcal{H}$.

Start with the bi-Gaussian (simulation as previous lecture):

In [42]:
m = Machine(phi_s=0)
In [43]:
sigma_z = 13.5
In [44]:
def generate_gaussian_in_rfbucket(N, sigma_z, machine, seed=12345, margin=0.05):
    '''Generate a bi-Gaussian distribution with N macro-particles,
    rms bunch length sigma_z and a matched sigma_deltap via the
    machine settings.
    '''
    np.random.seed(seed)
    m = machine
    
    sigma_deltap = np.sqrt(
        2 * m.p0() / -m.eta(0) * 
        charge * m.voltage * np.pi * m.harmonic / (beta(gamma(m.p0())) * c * m.circumference**2)
    ) * sigma_z

    z_ini = np.random.normal(loc=0, scale=sigma_z, size=N)
    deltap_ini = np.random.normal(loc=0, scale=sigma_deltap, size=N)
    
    H_safetymargin = margin * hamiltonian(0, 0, m)
    H_values = hamiltonian(z_ini, deltap_ini, m) - H_safetymargin

    while any(H_values >= 0):
        mask_bad = H_values >= 0
        N_bad = np.sum(mask_bad)
        print (N_bad)
        # re-initialise bad particles:
        z_ini[mask_bad] = np.random.normal(loc=0, scale=sigma_z, size=N_bad)
        deltap_ini[mask_bad] = np.random.normal(loc=0, scale=sigma_deltap, size=N_bad)
        # re-evaluate rejection condition
        H_values = hamiltonian(z_ini, deltap_ini, m) - H_safetymargin
    
    return z_ini, deltap_ini
In [45]:
N = 10000
n_turns = 5000
In [46]:
z_ini, deltap_ini = generate_gaussian_in_rfbucket(N, sigma_z, m)

860
83
4
1
In [47]:
plot_mp(z_ini, deltap_ini / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=20);

Track the bi-Gaussian distribution...

In [48]:
z = np.zeros((n_turns, N), dtype=np.float64)
deltap = np.zeros_like(z)

z[0] = z_ini
deltap[0] = deltap_ini
In [49]:
for i_turn in trange(1, n_turns):
    z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)

  0%|          | 0/4999 [00:00<?, ?it/s]

The rms emittance evolution:

In [50]:
epsn_z = np.array([emittance(z_i, deltap_i) for z_i, deltap_i in zip(z, deltap)])

ylim_m = np.median(4 * np.pi * epsn_z / e)
ylim_d = 1.1 * np.max(np.abs(ylim_m - 4 * np.pi * epsn_z / e))
In [51]:
plt.plot(4 * np.pi * epsn_z / e)

plt.ylim(ylim_m - ylim_d, ylim_m + ylim_d)

plt.xlabel('Turns')
plt.ylabel('$4\pi\epsilon_z$ [eV.s]');

$\leadsto$ the Gaussian particle distribution is not exactly in equilibrium for sufficiently large rms values in the nonlinear potential, the particles filament and the rms emittance grows (a little)!

$\implies$ compare to using full nonlinear Hamiltonian to construct PDF $\psi(\mathcal{H})\propto\exp\left(\cfrac{\mathcal{H}}{\mathcal{H}_0}\right)$

In [52]:
rfb = get_pyht_rfbucket(m)

rfb_matcher = RFBucketMatcher(rfb, ThermalDistribution, sigma_z=sigma_z)
rfb_matcher.integrationmethod = 'cumtrapz'
In [53]:
z_ini, delta_ini, _, _ = rfb_matcher.generate(N)

*** Maximum RMS bunch length 14.146533693849044m.
... distance to target bunch length: -1.3365e+01
... distance to target bunch length: 2.5641e-01
... distance to target bunch length: 2.4078e-01
... distance to target bunch length: -1.8484e-01
... distance to target bunch length: 5.7289e-02
... distance to target bunch length: 1.0673e-02
... distance to target bunch length: -1.5085e-04
... distance to target bunch length: 1.9092e-06
... distance to target bunch length: -4.9149e-06
--> Bunch length: 13.500001909209015
--> Emittance: 5.835458508892962
In [54]:
deltap_ini = delta_ini * m.p0()
In [55]:
plot_mp(z_ini, delta_ini, rfb, n_bins=20);

Track the matched thermal distribution...

In [56]:
z = np.zeros((n_turns, N), dtype=np.float64)
deltap = np.zeros_like(z)

z[0] = z_ini
deltap[0] = deltap_ini
In [57]:
for i_turn in trange(1, n_turns):
    z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)

  0%|          | 0/4999 [00:00<?, ?it/s]

The rms emittance evolution:

In [58]:
epsn_z = np.array([emittance(z_i, deltap_i) for z_i, deltap_i in zip(z, deltap)])

ylim_m = np.median(4 * np.pi * epsn_z / e)
In [59]:
plt.plot(4 * np.pi * epsn_z / e)

plt.ylim(ylim_m - ylim_d, ylim_m + ylim_d)

plt.xlabel('Turns')
plt.ylabel('$4\pi\epsilon_z$ [eV.s]');

$\implies$ this result shows that the nonlinearly matched thermal distribution is in equilibrium from the start (up to macro-particle noise, the fluctuations reduce with $1/\sqrt{N}$)!

Section B: Physical Emittance Growth – Dipole Injection Mismatch

In [60]:
m = Machine(phi_s=0)

sigma_z = 8

z_ini, deltap_ini = generate_gaussian_in_rfbucket(N, sigma_z, m, margin=0.15)

24

Simulate a dipole injection mismatch (e.g. when the rf phase is not well synchronised between the injector and the synchrotron):

In [61]:
z_ini -= 0.5 * sigma_z

4 meter mismatch in $z$ correspond to a phase mismatch of 16 degree:

In [62]:
4 / (m.circumference / m.harmonic) * 360
Out[62]:
16.04281826366305
In [63]:
plot_mp(z_ini, deltap_ini / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=20);

$\implies$ note the offset towards negative $z$, the contours of the macro-particle density are no longer matched to the Hamiltonian contours.

The safety margin inside the separatrix (where no particles are generated in generate_gaussian_in_rfbucket) should be chosen large enough such that no particles are located outside the rf bucket after the mismatch:

In [64]:
assert all(hamiltonian(z_ini, deltap_ini, m) < 0), 'particles have been generated outside the rf bucket!'

Tracking the mismatched distribution of macro-particles:

In [65]:
z = np.zeros((n_turns, N), dtype=np.float64)
deltap = np.zeros_like(z)

z[0] = z_ini
deltap[0] = deltap_ini
In [66]:
for i_turn in trange(1, n_turns):
    z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)

  0%|          | 0/4999 [00:00<?, ?it/s]

Centroid Results

In [67]:
plt.plot(np.mean(z, axis=1))

plt.xlabel('Turns')
plt.ylabel(r'$\langle z \rangle$ [m]');

$\implies$ exponential decay of the initial offset (due to the non-linearity of the rf bucket)

RMS Bunch Length Results

In [68]:
plt.plot(np.std(z, axis=1))

plt.xlabel('Turns')
plt.ylabel(r'$\sigma_z$ [m]');

$\implies$ saturation of the rms bunch length growth

RMS Emittance Results

In [69]:
epsn_z = np.array([emittance(z_i, deltap_i) for z_i, deltap_i in zip(z, deltap)])
In [70]:
plt.plot(epsn_z / e)

plt.xlabel('Turns')
plt.ylabel('$\epsilon_z$ [eV.s]');

$\implies$ in this example, 10% emittance growth as a result of the 4 meter injection offset.

In [71]:
plot_mp(z[-1], deltap[-1] / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=40);

$\implies$ the filamentation of the macro-particle distribution is clearly visible!

Section C: Physical Emittance Growth – Quadrupole Injection Mismatch

In [72]:
m = Machine(phi_s=0)

sigma_z = 8

z_ini, deltap_ini = generate_gaussian_in_rfbucket(N, sigma_z, m, margin=0.15)

24

Simulate a quadrupole injection mismatch (e.g. when the rf voltage (rf bucket height) is not matched between the injector and the synchrotron):

In [73]:
deltap_ini *= 0.5
In [74]:
plot_mp(z_ini, deltap_ini / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=20);

$\implies$ note the squeezed rms momentum spread, the contours of the macro-particle density are no longer matched to the Hamiltonian contours.

Tracking the mismatched distribution of macro-particles:

In [75]:
z = np.zeros((n_turns, N), dtype=np.float64)
deltap = np.zeros_like(z)

z[0] = z_ini
deltap[0] = deltap_ini
In [76]:
for i_turn in trange(1, n_turns):
    z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)

  0%|          | 0/4999 [00:00<?, ?it/s]

Centroid Results

In [77]:
plt.plot(np.mean(z, axis=1))

plt.xlabel('Turns')
plt.ylabel(r'$\langle z \rangle$ [m]');

$\implies$ only residual centroid fluctuations (due to macro-particle noise), note the amplitude of the oscillation in comparison to the rf bucket length!

RMS Bunch Length Results

In [78]:
plt.plot(np.std(z, axis=1))

plt.xlabel('Turns')
plt.ylabel(r'$\sigma_z$ [m]');

$\implies$ exponential decay of the initial momentum mismatch (due to the non-linearity of the rf bucket)

RMS Emittance Results

In [79]:
epsn_z = np.array([emittance(z_i, deltap_i) for z_i, deltap_i in zip(z, deltap)])
In [80]:
plt.plot(epsn_z / e)

plt.xlabel('Turns')
plt.ylabel('$\epsilon_z$ [eV.s]');

$\implies$ in this example, 20% emittance growth as a result of the 50% momentum spread mismatch.

In [81]:
plot_mp(z[-1], deltap[-1] / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=40);

$\implies$ again, the filamentation of the macro-particle distribution is clearly visible!

Summary

  • full simulation of acceleration ramp, incl. transition crossing (and necessary synchronous phase adjustments!)
  • equilibrium distributions with nonlinear Hamiltonian
  • matching algorithm to find $\mathcal{H}_0$ (given target rms bunch length $\hat{\sigma}_z$)
  • root-finding problem (Brent-Dekker or secant algorithm to find root)
  • emittance growth mechanisms (numerical + physical: dipole and quadrupole injection mismatch)

PS: Solving the Exercise in Part I

To avoid the emittance growth in the simulation of the acceleration ramp in the CERN Proton Synchrotron, the operators should switch the synchronous phase at transition energy in order to restore phase focusing. (The previous lectures provide the conceptual background to understand what is going on.)

The switch in the synchronous phase can be implemented in the tracking loop as follows:

  • compute the new synchronous phase for the relativistic regime from the initial m.phi_s (which determines the energy gain) and store the value in a variable, e.g. phi_s_1, before the tracking starts
  • define the condition where transition energy is crossed
  • if this condition is met, set the synchronous phase m.phi_s to the correct new value phi_s_1
phi_s_1 =  # compute here the correct new synchronous phase above transition
# by using the one from below transition which is m.phi_s

for i_turn in trange(1, n_turns):
    z, deltap = track_one_turn(z, deltap, m)

    epsn_z[i_turn] = emittance(z, deltap)

    condition =  # write here X > Y after the "="
    if condition:
        m.phi_s = phi_s_1