Longitudinal dynamics of electron injection into a linac

A. Petrenko (CERN, 2018)

Consider longitudinal dynamics of electron injection into a linear accelerator (e.g. plasma wakefield) with electric field

$$E_z(z,t) = E_a \cos (k z -\omega t),$$

where wavevector $k = \omega / c = 2\pi/\lambda$. Defining $\xi = z - ct$,

$$E_z(\xi,t) = E_a \cos k \xi.$$

Assuming only longitudinal movement of electron the equation of its motion then can be written as

$$\begin{cases} \frac{dp_z}{dt} = -eE_a \cos k \xi \\ \frac{d\xi}{dt} = v_z - c \\ p_z = \gamma m v_z \end{cases}$$

($e>0$ is assumed here). Since $\frac{v_z^2}{c^2} = \left(1 + \frac{m^2c^2}{p_z^2}\right)^{-1}$ therefore

$$\begin{cases} \frac{dp_z}{cdt} = -\frac{eE_a}{c} \cos k \xi \\ \frac{d\xi}{cdt} = \left(1 + \frac{m^2c^2}{p_z^2} \right)^{-1/2} - 1. \end{cases}$$

The dependence of $\xi(p_z)$, i.e. velosity bunching, can be found if we divide one equation by another:

$$\frac{eE_a}{c} \cos k \xi d\xi = dp_z -\left(1 + \frac{m^2c^2}{p_z^2} \right)^{-1/2}dp_z.$$

Both sides are easy to integrate:

$$\frac{eE_a}{kc} d(\sin k \xi) = dp_z - mc d\sqrt{1 + \frac{p_z^2}{m^2c^2}}.$$

And the solution (exact!) is

$$\frac{eE_a}{kmc^2} ( \sin k \xi - \sin k\xi_0 ) = \frac{p_z - p_0}{mc} - \sqrt{1 + \frac{p_z^2}{m^2c^2}} + \sqrt{1 + \frac{p_0^2}{m^2c^2}}.$$

Let's plot this solution for different initial conditions $(\xi_0, p_0)$ and $E_a$:

import numpy as np
import holoviews as hv

hv.extension('bokeh')

Define some useful constants first

wavelength = 1.26 # mm
k = 2*np.pi/wavelength # 1/mm
mc = 0.511 # MeV/c

$$\sin k \xi = \sin k\xi_0 + \frac{k mc^2}{e E_a} \left( \frac{p_z - p_0}{mc} - \sqrt{1 + \frac{p_z^2}{m^2c^2}} + \sqrt{1 + \frac{p_0^2}{m^2c^2}} \right).$$

dim_xi = hv.Dimension('xi', label="ξ", unit='mm', range=(-0.5,+0.8))
dim_E  = hv.Dimension('E',  unit='MV/m', label='-E', range=(-950,+950))

xi_range = np.linspace(dim_xi.range[0],dim_xi.range[1], 30) # mm

%opts Curve [width=500 height=150 show_grid=True]
%opts HLine (color='black', line_width=0.5)

def xi_E_plot(Ea):
    E = -Ea*np.cos(xi_range*k) # MV/m
    xi_E = hv.Curve((xi_range,E), kdims=dim_xi, vdims=dim_E)*hv.HLine(0)
    return xi_E

xi_E_plot(700)

pz_max = 50 # MeV/c
pz = np.linspace(0,pz_max,100) # MeV/c

def f_sin_kxi(pz,xi0,p0,Ea):
    return np.sin(k*xi0) + (mc*k/(Ea*1e-3))*((pz-p0)/mc - np.sqrt(1 + pz*pz/mc/mc) + np.sqrt(1+p0*p0/mc/mc))

%opts Scatter [width=500 height=300 show_grid=True]

dim_pz = hv.Dimension('pz', unit='MeV/c', range=(0, pz_max))
dim_sin_kxi = hv.Dimension('sin_kxi', label="sin(kξ)", range=(-1.2,+1.2))

p0 = 20 # MeV/c
xi0 = 0.0 # mm

Ea = -800.0 # MV/m

hv.Scatter((f_sin_kxi(pz,xi0,p0,Ea),pz), kdims=[dim_sin_kxi, dim_pz])

Now define the function producing the plot for a given set of parameters.

def plot_xi_pz(pz,xi0,p0,Ea):
    sin_kxi = f_sin_kxi(pz,xi0,p0,Ea)
    pz      =      pz[np.abs(sin_kxi)<=1]
    sin_kxi = sin_kxi[np.abs(sin_kxi)<=1]
    xi = np.arcsin(sin_kxi)/k
    xi_pz  = hv.Scatter((xi,pz), kdims=[dim_xi, dim_pz])

    xi0_p0 = hv.Scatter(([xi0],[p0]), kdims=[dim_xi, dim_pz]).options(size=7)
    xi_E = xi_E_plot(Ea)

    return (xi_E + xi_pz*xi0_p0).cols(1)

#plot_xi_pz(pz,xi0,p0,Ea)

Interactive version:

#xi0_points = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6] # mm
#items = [(xi0, plot_xi_pz(pz,xi0,p0,Ea)) for xi0 in xi0_points]

#m = hv.HoloMap(items, kdims = ['ξ0']).collate()
#m

xi0_points = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6] # mm
p0_points = [3,5,10,20] # MeV/c
Ea_points = [200,500,800] # MeV/c

dim_p0 = dim_pz.clone()
dim_p0.label = 'p0'

items = [( (xi0,p0,Ea), plot_xi_pz(pz,xi0,p0,-Ea) ) for xi0 in xi0_points  for p0 in p0_points  for Ea in Ea_points]

m = hv.HoloMap(items, kdims = [dim_xi, dim_p0, dim_E]).collate()

m

ξ: 0 mm, p0: 3 MeV/c, -E: 200 MV/m

ξ (mm):

p0 (MeV/c):

-E (MV/m):

As we can see for a high-amplitude wakefield velocity bunching can be used only starting with a very low-energy e-beam (few MeV).

The effect of velocity bunching is defined by the dependense of $\xi(\xi_0, p_z \gg p_0)$. Since $x -\sqrt{x^2 + 1} \rightarrow 0$, if $x \rightarrow \infty$,

$$\sin k \xi_{p_z \rightarrow \infty} = \sin k\xi_0 + \frac{k mc^2}{e E_a} \left(\sqrt{1 + \frac{p_0^2}{m^2c^2}} - \frac{p_0}{mc}\right) \approx \sin k\xi_0 + \frac{k mc^2}{e E_a} \cdot \frac {mc}{2 p_0}.$$

Increasing accelerating field and initial electron momentum has the same effect on velocity bunching.

Now let's plot $\xi_{p_z \rightarrow \infty}$ vs $\xi_0$ for different $p_0.$

xi0 = np.linspace(0, wavelength/2, 100) # region of focusing plasma wakefield

# We can use the same function f_sin_kxi,
# but now xi0 is a vector instead of pz which is set to some large value:
pz = 1000 # MeV/c

hv.Scatter((xi0,f_sin_kxi(pz,xi0,p0,Ea)), kdims=[dim_xi, dim_sin_kxi])

dim_xi_final = hv.Dimension('xi_f', label="ξ(high energy)", unit='mm', range=(-0.32,+0.32))

def plot_xi0_xi(pz,xi0,p0,Ea):
    sin_kxi = f_sin_kxi(pz,xi0,p0,Ea)
    
    xi0      =    xi0[np.abs(sin_kxi)<=1.0]
    sin_kxi = sin_kxi[np.abs(sin_kxi)<=1.0]
    
    xi = np.arcsin(sin_kxi)/k
    
    xi0_xi  = hv.Scatter((xi0,xi), kdims=[dim_xi, dim_xi_final])

    xi_E = xi_E_plot(Ea)

    return (xi_E + xi0_xi*hv.HLine(0)).cols(1)

#plot_xi0_xi(pz,xi0,p0,Ea)

p0_points = [0.5,1,2,3,4,5,10,20,50] # MeV/c
Ea_points = [200,500,800] # MeV/c

items = [( (p0,Ea), plot_xi0_xi(pz,xi0,p0,-Ea) ) for p0 in p0_points for Ea in Ea_points]

m = hv.HoloMap(items, kdims = [dim_p0, dim_E]).collate()

m

p0: 0.5 MeV/c, -E: 200 MV/m

p0 (MeV/c):

-E (MV/m):