Longitudinal dynamics of muon injection into a linac

A. Petrenko (CERN, 2018)

Consider longitudinal dynamics of muon 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 muon 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$:

In [1]:
import numpy as np
import holoviews as hv

hv.extension('bokeh')

Define some useful constants first

In [2]:
wavelength = 1.26 # mm
k = 2*np.pi/wavelength # 1/mm
mc = 0.105658 # GeV/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). $$
In [3]:
dim_xi = hv.Dimension('xi', label="ξ", unit='mm', range=(-0.5,+0.8))
dim_E  = hv.Dimension('E',  unit='GV/m', label='-E', range=(-2.1,+2.1))
In [4]:
xi_range = np.linspace(dim_xi.range[0],dim_xi.range[1], 30) # mm
In [5]:
%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) # GV/m
    xi_E = hv.Curve((xi_range,E), kdims=dim_xi, vdims=dim_E)*hv.HLine(0)
    return xi_E

xi_E_plot(-0.5)
Out[5]:
In [6]:
pz_max = 300 # GeV/c
pz = np.linspace(0,pz_max,100) # GeV/c
In [7]:
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))
In [8]:
%opts Scatter [width=500 height=300 show_grid=True]

dim_pz = hv.Dimension('pz', unit='GeV/c', range=(0, pz_max))
dim_sin_kxi = hv.Dimension('sin_kxi', label="sin(kξ)", range=(-1.2,+1.2))
In [9]:
p0 = 20 # GeV/c
xi0 = 0.0 # mm

Ea = -1.0 # GV/m

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

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

In [10]:
%opts Scatter.MARK (size=7, color='red')

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, vdims=dim_pz)

    xi0_p0 = hv.Scatter((xi0,p0), kdims=dim_xi, vdims=dim_pz, group="MARK")
    xi_E = xi_E_plot(Ea)

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

#plot_xi_pz(pz,xi0,p0,Ea)

Interactive version:

In [11]:
#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
In [12]:
xi0_points = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6] # mm
p0_points = [20,40,100] # GeV/c
Ea_points = [0.5,1.0,2.0] # GeV/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()
In [13]:
m
Out[13]:

As we can see velocity bunching is important for $\mu$-beam below 100 GeV/c.

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 muon momentum has the same effect on velocity bunching.

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

In [14]:
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 = 3000 # GeV/c

hv.Scatter((xi0,f_sin_kxi(pz,xi0,p0,Ea)), kdims=dim_xi, vdims=dim_sin_kxi)
Out[14]:
In [15]:
dim_xi_final = hv.Dimension('xi_f', label="ξ(high energy)", unit='mm', range=(-0.32,+0.32))
In [16]:
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, vdims=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)
Out[16]:
In [17]:
p0_points = [20,50,100] # GeV/c
Ea_points = [0.7,1.0,2.0] # GV/m

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()
In [18]:
m
Out[18]:
In [ ]: