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$:
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.105658 # GeV/c
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))
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) # 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)
pz_max = 300 # GeV/c
pz = np.linspace(0,pz_max,100) # GeV/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='GeV/c', range=(0, pz_max))
dim_sin_kxi = hv.Dimension('sin_kxi', label="sin(kξ)", range=(-1.2,+1.2))
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)
Now define the function producing the plot for a given set of parameters.
%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:
#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 = [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()
m
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.$
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)
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, 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)
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()
m