Релятивисткая разностная схема для расчета динамики частиц в сложных электрических и магнитных полях

First, we write the main quantities in the difference scheme, then we consider the algorithm.

Basic values

At the initial time \(t\) set Cartesian coordinates of the \(i -\)particles \(\overrightarrow{r_{i}} = (x_i, y_i, z_i)\)nd the initial pulses \(\overrightarrow{p_{i}} = (p_{x_i}, p_{y_i}, p_{z_i})\). Minimum time step is set \(\delta t\).

For convenience, made dimensionless physical quantities (with tilde): \begin{equation} \widetilde{\overrightarrow{r_i}} = \frac{\overrightarrow{r_i}}{c \delta t}, \end{equation} \begin{equation} \widetilde{\overrightarrow{E_i}} = \frac{\overrightarrow{E_i}e}{mc} \delta t, \end{equation} \begin{equation} \widetilde{\overrightarrow{H_i}} = \frac{\overrightarrow{H_i}e}{mc} \delta t, \end{equation} where \(m\) is the speed of light, \(e\) is the particle charge and \(m\) is the particle mass.

Algorithm

Consider the relativistic difference scheme algorithm step by step:

• Calculation of electric fields \(\overrightarrow{E_i}\) at the points where the particles are.

• Pulse increment: \begin{equation} \overrightarrow{p_i} = \overrightarrow{p_i} + 2\cdot\overrightarrow{E_i}, \end{equation} where the coefficient means that the increment of the pulse is made for the entire time interval \(2\delta t\).

• New speeds are calculated by new impulses: \begin{equation} \overrightarrow{v_i} = \frac{\overrightarrow{p_i}}{\gamma_i}, \end{equation} where \(\gamma_i = \sqrt{1 + \overrightarrow{p_i}^2}.\)

• At new speeds, new coordinates are calculated: \begin{equation} \overrightarrow{r_i} = \overrightarrow{r_i} + 1\cdot\overrightarrow{v_i}, \end{equation} where the coefficient means that the increment of the coordinate is made for the time interval \(\delta t.\)

• Calculation of magnetic fields \(\overrightarrow{H_i}\) at the new points where the particles are.

• Calculated velocity values (after rotation in a magnetic field):

  \[ b_1 = 1 - \frac{H_i^2}{\gamma_i}, b_2 = 1 + \frac{H_i^2}{\gamma_i}, b_3 = 2\cdot \frac{\overrightarrow{v_i}\cdot\overrightarrow{H_i}}{\gamma_i}, \]

  \[ \overrightarrow{f_i} = 2\cdot \frac{\overrightarrow{v_i}\times \overrightarrow{H_i}}{\gamma_i}, \]

  \begin{equation} \overrightarrow{v_i} = \frac{\overrightarrow{v_i}b_1 + \overrightarrow{f_i} + \frac{\overrightarrow{H_i}}{\gamma_i}b_3}{b_2}. \end{equation}

• At new speeds, new coordinates are calculated: \begin{equation} \overrightarrow{r_i} = \overrightarrow{r_i} + 1\cdot\overrightarrow{v_i}, \end{equation} where the coefficient means that the increment of the coordinate is made for the time interval $\delta t.$

• New impulses are calculated according to the new rates: \begin{equation} \overrightarrow{p_i} = \overrightarrow{v_i}\gamma_i. \end{equation}

One cycle is performed in a time interval \(2\delta t.\)

Python

import redpic as rp
import numpy as np
import holoviews as hv
hv.extension('matplotlib')

import warnings
warnings.filterwarnings('ignore')

rp.__version__

Some plotting options

%output size=100 backend='matplotlib' fig='png' dpi=300
%opts Curve Scatter [aspect=3 show_grid=True]
%opts Curve (linewidth=1 alpha=0.7 color='blue')
%opts Scatter (alpha=0.7 s=0.5)

Define accelerator beamline parameters

acc = rp.accelerator.Accelerator(0.7, 5, 0.01)

#              Unique name,  z-position [m],  Ez [MV/m],  Ez(z) profile
acc.add_accel('Acc. 1',      4.096,          -1.1,         'Ez.dat')
acc.add_accel('Acc. 2',      5.944,          -1.1,         'Ez.dat')

#                 Unique name,  z-position [m],  Bz [T],  Bz(z) profile
acc.add_solenoid('Sol. 1',      0.450,          -0.0580,   'Bz.dat')
acc.add_solenoid('Sol. 2',      0.957,           0.0390,   'Bz.dat')
acc.add_solenoid('Sol. 3',      2.107,           0.0250,   'Bz.dat')
acc.add_solenoid('Sol. 4',      2.907,           0.0440,   'Bz.dat')
acc.add_solenoid('Sol. 5',      3.670,           0.0400,   'Bz.dat')
acc.add_solenoid('Sol. 6',      4.570,           0.0595,   'Bz.dat')
acc.add_solenoid('Sol. 7',      5.470,           0.0590,   'Bz.dat')

acc.compile()

dim_z  = hv.Dimension('z',  unit='m')
dim_Ez = hv.Dimension('Ez', unit='MV/m', label='$E_z$')
dim_Bz = hv.Dimension('Bz', unit='Gs', label='$B_z$')

z  = acc.z
z_Ez = hv.Curve((z, acc.Ez(z)), kdims=dim_z, vdims=dim_Ez)
z_Bz = hv.Curve((z, acc.Bz(z)*1e4), kdims=dim_z, vdims=dim_Bz)

(z_Ez + z_Bz).cols(1)

print(acc)

Accelerator structure.
	Solenoids:
	[ 0.45 m, -0.058 T, Bz.dat, Sol. 1, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 0.957 m, 0.039 T, Bz.dat, Sol. 2, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 2.107 m, 0.025 T, Bz.dat, Sol. 3, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 2.907 m, 0.044 T, Bz.dat, Sol. 4, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 3.67 m, 0.04 T, Bz.dat, Sol. 5, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 4.57 m, 0.0595 T, Bz.dat, Sol. 6, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 5.47 m, 0.059 T, Bz.dat, Sol. 7, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	Accelerating modules:
	[ 4.096 m, -1.1 T, Ez.dat, Acc. 1, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 5.944 m, -1.1 T, Ez.dat, Acc. 2, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	Quadrupoles:
	Correctors x:
	Correctors y:

Define beam parameters

beam = rp.beam.RadialUniformBeam(
    type=rp.constants.electron, 
    energy = 1.32,          # MeV
    current = 0.5e3,  # A
    radius_x = 48e-3, # initial r (m)
    radius_y = 48e-3, # initial r (m)
    radius_z = 3.5,
    radius_xp = 2*35.0e-3,     # initial r' (rad)
    radius_yp = 2*35.0e-3,     # initial r' (rad)
    x  = 0.0e-3,   # horizontal centroid position (m)
    xp = 0.0e-3,     # horizontal centroid angle (rad)
    y = 0,          # vertical centroid position (m)
    normalized_emittance = 200e-6 # m*rad
)

beam.generate(10_000)

2023-08-31 21:35:33,649 - redpic.beam.base - INFO - Generate a radial-uniform beam with 10000 particles

beam.df

dim_x = hv.Dimension('x', unit='m', range=(-0.1, 0.1))
dim_y = hv.Dimension('y', unit='m', range=(-0.1, 0.1))
dim_z = hv.Dimension('z', unit='m', range=(acc.z_start, acc.z_stop))
dim_px = hv.Dimension('px', unit='MeV/c', label='$p_x$')
dim_py = hv.Dimension('py', unit='MeV/c', label='$p_y$')

beam_x_y = hv.Scatter(beam.df, kdims=[dim_x, dim_y])
beam_z_x = hv.Scatter(beam.df, kdims=[dim_z, dim_x])
beam_x_px = hv.Scatter(beam.df, kdims=[dim_x, dim_px])
beam_y_py = hv.Scatter(beam.df, kdims=[dim_y, dim_py])

WARNING:param.Scatter: Chart elements should only be supplied a single kdim
2023-08-31 21:35:33,772 - param.Scatter - WARNING - Chart elements should only be supplied a single kdim
WARNING:param.Scatter: Chart elements should only be supplied a single kdim
2023-08-31 21:35:33,789 - param.Scatter - WARNING - Chart elements should only be supplied a single kdim
WARNING:param.Scatter: Chart elements should only be supplied a single kdim
2023-08-31 21:35:33,791 - param.Scatter - WARNING - Chart elements should only be supplied a single kdim
WARNING:param.Scatter: Chart elements should only be supplied a single kdim
2023-08-31 21:35:33,801 - param.Scatter - WARNING - Chart elements should only be supplied a single kdim

(beam_x_y + beam_z_x + beam_x_px + beam_y_py).cols(2)

print(beam)

Beam parameters:
            Type	electron
            Distribution	radial-uniform
            Particles	10000
            Current	500.0 A
            Energy	1.32 MeV
            Total momentum	1.7582483378931433 MeV/c
            Rel. factor	3.583175582013202
            Radius x	48.0 mm
            Radius y	48.0 mm
            Radius z	3.5 m
            Radius x prime	70.0 mrad
            Radius y prime	70.0 mrad
            Horizontal centroid position	0.0 mm
            Vertical centroid position	0.0 mm
            Horizontal centroid angle	0.0 mrad
            Vertical centroid angle	0.0 mrad
            Normalized emittance x	200.0 mm*mrad
            Normalized emittance y	200.0 mm*mrad

Run simulation

rp_sim = rp.solver.Simulation(beam, acc)
rp_sim.track()

z = 4.98 m (99.5 %) 

Plot the simulation results

def plot(i):
    df = rp_sim.result[i]
    rp_z_x = hv.Scatter(df, kdims=[dim_z, dim_x], label='redpic')
    return rp_z_x

items = [(i, plot(i)) for i in list(rp_sim.result.keys())]

holomap = hv.HoloMap(items, kdims = ['z'])

hv.output(holomap, widget_location='bottom')