Development of Beam Optics Code LEADS-v5

Jianqin Lü

State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing, China.

DOI: 10.4236/jamp.2022.104089                

Abstract

To calculate nonlinear transport of space charge dominated beam in 6D phase spaces, a computer code package LEADS-v5 (Linear and Electrostatic Accelerator Dynamics Simulations) has been developed. The codes calculate particle motions in the beam transport systems consisting of electrostatic and magnetic focusing lenses, ion analyzers, multipoles and RF accelerating structures. The nonlinear forces of external electric/magnetic fields are analyzed by the Lie algebraic method, and the space charge forces are obtained by the particle in cell (PIC) scheme. In the codes, Uniform and Gaussian particle distributions can be chosen to generate randomly the particle initial coordinates. The optimization procedures are provided to make the beam optics designs reasonable and fast. Graphically displays of calculated results are provided.

Keywords

Beam Transport, Nonlinear Trajectories, Lie Map, PIC Scheme, Computer Code

1. Introduction

Intense beam nonlinear transport is a very complicated problem, because the state of particle motion is dominated not only by the applied electromagnetic fields, but also by the beam induced electromagnetic fields (self fields). Moreover, the self fields are related to the beam dimensions and particle distributions. So, it is impossible to get the self-consistent solutions of particle motions analytically. To solve this problem, we use two methods: Lie algebraic method and PIC method.

In 1987 Dr. Dragt first introduced the Lie algebraic theory [1] [2] into accelerator physics. Since then the theory has been widely used by accelerator physicists in the research of charged particle nonlinear motions. The PIC method was originally developed for fluid dynamics simulations [3] [4]; however, its greatest application is in plasma physics [5] [6]. The method began to be used for particle accelerator simulations two decades ago. We combine the Lie algebraic method and the PIC method together to simulate intense beam nonlinear transport. With the Lie algebraic method we analyze the particle nonlinear trajectories in the applied electromagnetic fields, and with the PIC algorithm we calculate the effects of space charge forces of the beam. Based on the two methods, we have developed a computer program LEADS-v5.

Optimization plays an important role in the computational design for the particle accelerators and beam optical lines. It will not only save much computational time, but also make the design more reasonable. Many accelerator codes can perform optimization calculations. For example, program TRANSPORT by K. L. Brown [7] has powerful optimization ability. For this reason, we put some optimization subroutines in LEADS-5v, by using the Powell nonlinear optimization method.

The calculated results, such as the beam phase areas both in transverse and longitudinal directions as well as the beam envelopes can be displayed graphically on the computer monitor, so that it makes the calculated results simple and visual.

2. Particle Distributions

The particle distributions can be selected by the user. Two kinds of distributions are provided: uniform and Gaussian distributions in the 6D (x, x', y, y', ΔE, Δf) phase spaces as the initial particle coordinates which are generated randomly.

3. Beam Focusing, Analyzing and Accelerating Elements

The following focusing, analyzing and accelerating elements are provided by LEADS-v5:

3.1. Electromagnetic Focusing Elements, Multipoles (To Fifth Order Approximation)

· Electrostatic quadrupoles.

· Solenoid magnets.

· Magnetic quadrupoles.

· Sextupole magnets.

· Octupole magnets.

· Decapole magnets.

· Dodecapole magnets.

3.2. Ion Analyzing Elements (To Fifth Order Approximation)

· Dipole magnets (including fringe fields).

· Cylindrical electrostatic analyzers.

· Spherical electrostatic analyzers.

· E × B analyzers.

3.3. Axisymmetric Electrostatic Lenses (To Third Order Approximation) [8] [9] [10] [11]

· Two-cylinder accelerating gap lenses.

· Three-cylinder Einzel lenses.

· Three-aperture Einzel lenses.

· DC accelerator columns.

3.4. RF Accelerating Cavities (To First Order Approximation)

· Quarter wave resonators (QWRs).

· Split loop resonators (SLRs).

3.5. Other Elements

· Element rotation; Every physical element can be rotated about the z-axis in arbitrary angle (positive/negative).

· Arbitrary matrix. If you have known the transfer matrix of an optical element, you can put in the matrix to the code input data directly.

· Phase space plotting. The uses can plot the phase space diagrams of the beam in x-, y-, z-directions any where along the beam line.

· Beam envelope plotting. Beam envelope plotting is provided in the code.

4. Analysis Calculation Tools and Formulas

4.1. Lie Algebraic Tools

In the canonical phase space $\xi =\left(x,p_x,y,p_y,z,p_z\right)$, if one doesn’t take the space charge effects into count, the Hamiltonian of the particle motion is

$H_t=-{\left[m_0^2{c}^4+c^2{\left(p_x-q{A}_x\right)}^2+c^2{\left(p_y-q{A}_y\right)}^2+c^2{\left(p_z-q{A}_z\right)}^2\right]}^{\frac{1}{2}}+q\psi$ (1)

where x, y and z are the particle coordinates in the real space; p_x, p_y and p_z are the particle canonical momentum components; A_x, A_y and A_z are the magnetic vector components; y is the electric potential; q is the particle charge, m₀ is the particle rest mass. Now, take $\varsigma =\left(x,x^{\prime },y,y^{\prime },\tau ,p_{\tau }\right)$ as a new canonical phase space, where $\tau =t-z/v_0$ is the time difference between the arbitrary particle and the reference particle, $p_{\tau }=p_t-p_t^0$, $p_t=-H_t$, $p_t^0$ is the value of $p_t$ for the reference particle. From phase space x to phase space z, transformation x ® z is canonical. In the new phase space z, the Hamiltonian is

$\begin{array}{l}H=-{\left[-{\left(p_x-q{A}_x\right)}^2-{\left(p_y-q{A}_y\right)}^2-{\left(p_t+p_t^0+q\psi \right)}^2/c^2-m_0^2{c}^2\right]}^{1/2}\\ +q{A}_z-\frac{p_t+p_t^0}{v_0}\end{array}$ (2)

and in the phase space $\zeta =\left(x,x^{\prime },y,y^{\prime },\tau ,p_{\tau }\right)$ the solutions of particle trajectories can be expressed as

$\begin{array}{l}{\zeta }_f{}_i={\displaystyle \underset{j=1}{\overset{6}{\sum }}{M}_{i,j}{\zeta }_j}+{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\displaystyle \underset{k=1}{\overset{6}{\sum }}{S}_{i,j,k}{\zeta }_j{\zeta }_k}}+{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\displaystyle \underset{k=1}{\overset{6}{\sum }}{\displaystyle \underset{l=1}{\overset{6}{\sum }}{T}_{i,j,k,l}{\zeta }_j{\zeta }_k{\zeta }_l}}}\\ +{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\displaystyle \underset{k=1}{\overset{6}{\sum }}{\displaystyle \underset{l=1}{\overset{6}{\sum }}{\displaystyle \underset{m=1}{\overset{6}{\sum }}{U}_{i,j,k,l,m}{\zeta }_j{\zeta }_k{\zeta }_l}}}{\zeta }_m}\\ +{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\displaystyle \underset{k=1}{\overset{6}{\sum }}{\displaystyle \underset{l=1}{\overset{6}{\sum }}{\displaystyle \underset{m=1}{\overset{6}{\sum }}{\displaystyle \underset{n=1}{\overset{6}{\sum }}{W}_{i,j,k,l,m,n}{\zeta }_j{\zeta }_k{\zeta }_l}}}{\zeta }_m{\zeta }_n}}+\cdots +\Delta {\zeta }_i,i=1,\cdots ,6\end{array}$ (3)

where z_f is the particle final coordinate, z the initial coordinate; $M_{i,j}$, $S_{i,j,k}$, $T_{i,j,k,l}$, $U_{i,j,k,l,m}$ and $W_{i,j,k,l,m,n}$ are the first order, second order and to the fifth order coefficients of the particle trajectories contributed by the applied electromagnetic fields, $\Delta \zeta$ is the contributions of the space charge forces to the particle trajectories.

As mentioned before, the particle trajectories in the external fields are analyzed with the Lie algebraic method, and the effects of space charge forces are calculated with the PIC scheme.

Expanding the Hamiltonian H into power series gives

$H=H_0+H_1+H_2+H_3+H_4+H_5+H_6+\cdots$ (4)

In the phase space, the final coordinate z_f and the initial coordinate z of a particle are related by a map M:

${\varsigma }_{{}_f}=\text{M}\varsigma$ (5)

According to the references [1] and [2], M is

$\text{M}=\cdots {\text{M}}_{\text{4}}{\text{M}}_{\text{3}}{\text{M}}_{\text{2}}$ (6)

where

$\begin{array}{l}{\text{M}}_{\text{2}}=\text{exp}\left(:f_{\text{2}}:\right),{\text{M}}_{\text{3}}=\text{exp}\left(:f_{\text{3}}:\right),{\text{M}}_{\text{4}}=\text{exp}\left(:f_{\text{4}}:\right),\\ {\text{M}}_{\text{5}}=\text{exp}\left(:f_{\text{5}}:\right),{\text{M}}_{\text{6}}=\text{exp}\left(:f_{\text{6}}:\right),\cdots \end{array}$ (7)

and

$f_2=-{\displaystyle {\int }_0^{\mathcal{l}}{H}_2}\left(z\right)\text{d}z$ (8)

$f_3=-{\displaystyle {\int }_0^{\mathcal{l}}{h}_3^{int}}\left(z\right)\text{d}z$ (9)

$f_4=-{\displaystyle {\int }_0^{\mathcal{l}}{h}_4^{int}}\left(z\right)\text{d}z$ (10)

$\begin{array}{c}{f}_5=-{\displaystyle {\int }_0^{\mathcal{l}}{h}_5^{int}}\left(z_1\right)\text{d}{z}_1+{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1}{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}\left[-h_3^{\mathrm{int}}\left(z_2\right),-h_4^{\mathrm{int}}\left(z_1\right)\right]\\ +\frac{1}{3}{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1}{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}{\displaystyle {\int }_{z_0}^{z_2}\text{d}{z}_3}(\left[-h_3^{\mathrm{int}}\left(z_3\right),\left[-h_3^{\mathrm{int}}\left(z_2\right),-h_3^{\mathrm{int}}\left(z_1\right)\right]\right],\\ +\left[-h_3^{\mathrm{int}}\left(z_2\right),\left[-h_3^{\mathrm{int}}\left(z_3\right),-h_3^{\mathrm{int}}\left(z_1\right)\right]\right])\end{array}$ (11)

$\begin{array}{l}{f}_6=-{\displaystyle {\int }_0^{\mathcal{l}}{h}_6^{int}}\left(z_1\right)\text{d}{z}_1+\frac{1}{4}{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1}{\displaystyle {\int }_{z_0}^{z{}_1}\text{d}{z}_2}\left[-h_3^{\mathrm{int}}\left(z_2\right),-h_5^{\mathrm{int}}\left(z_1\right)\right]\\ +\frac{1}{2}{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1}{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}\left[-h_4^{\mathrm{int}}\left(z_2\right),-h_4^{\mathrm{int}}\left(z_1\right)\right]\\ +\frac{1}{4}({\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1}{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}{\displaystyle {\int }_{z_0}^{z_2}\text{d}{z}_3}(\left[-h_4^{\mathrm{int}}\left(z_3\right),\left[-h_3^{\mathrm{int}}\left(z_2\right),-h_3^{\mathrm{int}}\left(z_1\right)\right]\right]\\ \end{array}$

$\begin{array}{l}+{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}{\displaystyle {\int }_{z_0}^{z_2}\text{d}{z}_3}(\left[-h_4^{\mathrm{int}}\left(z_2\right),\left[-h_3^{\mathrm{int}}\left(z_3\right),-h_3^{\mathrm{int}}\left(z_1\right)\right]\right]}\\ +3{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}{\displaystyle {\int }_{z_0}^{z_2}\text{d}{z}_3}(\left[-h_3^{\mathrm{int}}\left(z_3\right),\left[-h_3^{\mathrm{int}}\left(z_2\right),-h_4^{\mathrm{int}}\left(z_1\right)\right]\right]}\end{array}$

$\begin{array}{l}+{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}{\displaystyle {\int }_{z_0}^{z_2}\text{d}{z}_3}(\left[-h_3^{\mathrm{int}}\left(z_2\right),\left[-h_3^{\mathrm{int}}\left(z_3\right),-h_4^{\mathrm{int}}\left(z_1\right)\right]\right]})\\ +\frac{1}{4}({\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1}{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}{\displaystyle {\int }_{z_0}^{z_2}\text{d}{z}_3}{\displaystyle {\int }_{z_0}^{z_3}(\left[-h_3^{\mathrm{int}}\left(z_3\right),\left[-h_3^{\mathrm{int}}\left(z_2\right),\left[-h_3^{\mathrm{int}}\left(z_4\right),-h_3^{\mathrm{int}}\left(z_1\right)\right]\right]\right]}\\ +{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1}{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}{\displaystyle {\int }_{z_0}^{z_2}\text{d}{z}_3}{\displaystyle {\int }_{z_0}^{z_3}(\left[-h_3^{\mathrm{int}}\left(z_4\right),\left[-h_3^{\mathrm{int}}\left(z_2\right),\left[-h_3^{\mathrm{int}}\left(z_3\right),-h_3^{\mathrm{int}}\left(z_1\right)\right]\right]\right]}\\ +{\displaystyle {\int }_{z_0}^{\mathcal{l}}\text{d}{z}_1}{\displaystyle {\int }_{z_0}^{z_1}\text{d}{z}_2}{\displaystyle {\int }_{z_0}^{z_2}\text{d}{z}_3}{\displaystyle {\int }_{z_0}^{z_3}(\left[-h_3^{\mathrm{int}}\left(z_4\right),\left[-h_3^{\mathrm{int}}\left(z_3\right),\left[-h_3^{\mathrm{int}}\left(z_2\right),-h_3^{\mathrm{int}}\left(z_1\right)\right]\right]\right]})\end{array}$ (12)

where ℓ is the length of the element, and

$h_n^{int}\left(\varsigma ,z\right)={\text{M}}_{\text{2}}{H}_n$ (13)

To simplify the expressions we define the multiple integrations as follows:

$\langle \begin{array}{c}i\\ \alpha \end{array},\begin{array}{c}j\\ \beta \end{array},\cdots ,\begin{array}{c}m\\ \epsilon \end{array},\rangle ={\displaystyle {\int }_{t^{in}}^t\text{d}{t}_i}{\displaystyle {\int }_{t^{in}}^{t_i}\text{d}{t}_j}\cdots {\displaystyle {\int }_{t^{in}}^{t_l}\text{d}{t}_m}\left[:-H_{\alpha }^{\mathrm{int}}\left(t_i\right)::-H_{\beta }^{\mathrm{int}}\left(t_j\right):\cdots :-H_{\epsilon }^{\mathrm{int}}\left(t_m\right):\right]$ (14)

Then, the Equations (9)-(12) can be rewritten as

$f_3=\left[\begin{array}{c}1\\ 3\end{array}\right]$ (15)

$f_4=\left[\begin{array}{c}1\\ 4\end{array}\right]+\frac{1}{2}\left[\begin{array}{c}2\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]$ (16)

$f_5=\left[\begin{array}{c}1\\ 5\end{array}\right]+\left[\begin{array}{c}2\\ 3\end{array},\begin{array}{c}1\\ 4\end{array}\right]+\frac{1}{3}\left[\begin{array}{c}3\\ 3\end{array},\left[\begin{array}{c}2\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]\right]+\frac{1}{3}\left[\begin{array}{c}2\\ 3\end{array},\left[\begin{array}{c}3\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]\right]$ (17)

$\begin{array}{c}{f}_6=\left[\begin{array}{c}1\\ 6\end{array}\right]+\left[\begin{array}{c}2\\ 3\end{array},\begin{array}{c}1\\ 5\end{array}\right]+\frac{1}{2}\left[\begin{array}{c}2\\ 4\end{array},\begin{array}{c}1\\ 4\end{array}\right]\\ +\frac{1}{4}\left(\left[\begin{array}{c}3\\ 4\end{array},\left[\begin{array}{c}2\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]\right]+\left[\begin{array}{c}2\\ 4\end{array},\left[\begin{array}{c}3\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]\right]+3\left[\begin{array}{c}3\\ 3\end{array},\left[\begin{array}{c}2\\ 3\end{array},\begin{array}{c}1\\ 4\end{array}\right]\right]+\left[\begin{array}{c}2\\ 3\end{array},\left[\begin{array}{c}3\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]\right]\right)\\ +\frac{1}{4}\left(\left[\begin{array}{c}3\\ 3\end{array},\left[\begin{array}{c}2\\ 3\end{array},\left[\begin{array}{c}4\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]\right]\right]+\left[\begin{array}{c}4\\ 3\end{array},\left[\begin{array}{c}2\\ 3\end{array},\left[\begin{array}{c}3\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]\right]\right]+\left[\begin{array}{c}4\\ 3\end{array},\left[\begin{array}{c}3\\ 3\end{array},\left[\begin{array}{c}2\\ 3\end{array},\begin{array}{c}1\\ 3\end{array}\right]\right]\right]\right)\end{array}$ (18)

The first order, second order, and to fifth order terms of the orbit solutions ( ${\zeta }^i,i=\text{1},\text{2},\cdots$ ) are expressed as

${\varsigma }^1=f_2\varsigma$ (19)

${\varsigma }^2=f_3{\varsigma }^1$ (20)

${\varsigma }^3=\left(:f_4:+\frac{1}{2}:f_3{:}^2\right){\varsigma }^1$ (21)

${\varsigma }^4=\left(:f_5:+:f_4::f_3:+\frac{1}{6}:f_3{:}^3\right){\varsigma }^1$ (22)

${\varsigma }^5=\left(:f_6:+:f_5::f_3:+\frac{1}{2}:f_4{:}^2+\frac{1}{2}:f_4::f_3{:}^2+\frac{1}{24}:f_3{:}^4\right){\varsigma }^1$ (23)

where:*: stands for the Poisson operation. For example, for the given two functions f an g,:f: g = [f, g].

4.2. The PIC Algorithms [5] [6]

The beam self-fields are calculated with the PIC scheme. The initial particle distributions are generated randomly.

Generally speaking, to simulate the particle motion in the beam self-field with the PIC method, the following steps would be taken:

Figure 1. Mesh generation.

Figure 2. A mesh cell.

1) Randomly generate the particle initial distributions.

2) Suppose the beam moves along the straight line (ignoring the curvilinear orbit due to the short time step Δt); Divide the beam into cubic mesh cell (see Figure 1 and Figure 2); Leth_x,h_y, and h_z be the grid widths in the x, y and z directions respectively; the grid node numbers are N_x, N_y and N_z in the x, y and z directions respectively.

3) Assign the charge q of each macroparticle (MP) in the cell to the cell nodes: the charge of every MP is distributed to the eight cell nodes. The portion of q assigned to the cell nodes is determined by the position of the MP within the cell in accordance with the following relationship:

${\rho }_{l,m,n}={\displaystyle \underset{p=1}{\overset{N_c}{\sum }}{\rho }_{\text{MP}}}\left(1-\frac{\left|x_p-x_l\right|}{h_x}\right)\left(1-\frac{\left|y_p-y_m\right|}{h_y}\right)\left(1-\frac{\left|z_p-z_n\right|}{h_z}\right)$ (24)

where l = i, i + 1, m = j, j + 1, n = k, k + 1 are integer grid indices; N_c is MP number in this cell; p is the sequence number of MP in the cell; (x_p, y_p, z_p) are the coordinates in the cell, (x_l, y_m, z_n) are the coordinates of a node; and

${\rho }_{MP}=\frac{q}{h_x{h}_y{h}_z}$. (25)

From Equation (24) we see that the MP is not a point, but has a finite volume h_x× h_y× h_z, and r_MP is just the charge density of a MP.

4) Solve the Poisson equation ${\nabla }^2\varphi =-\rho \left(x,y,z\right)/{\epsilon }_0$ (ε₀ is the vacuum permittivity): After the charge of all of the MPs has been assigned to the mesh nodes, we use the Green’s function to calculate the potentials f_i.j,k at the mesh nodes. Usually, the beam transverse dimensions are much smaller than the vacuum pipe, so we can use the opening boundary conditions, and take ${\varphi }_{r\to \infty }=0$.

The solutions of the Poisson equation are expressed as:

$\varphi \left(x,y,z\right)=-\frac{1}{4\pi {\epsilon }_0}{\displaystyle \int G\left(x,x_0,y,y_0,z,z_0\right)\rho \left(x_0,y_0,z_0\right)\text{d}{x}_0\text{d}{y}_0\text{d}{z}_0}$ (26)

here $\left(x,y,z\right)$ is field point, $\left(x_0,y_0,z_0\right)$ is source point, and the Green’s function is

$G\left(x,x_0,y,y_0,z,z_0\right)=\frac{1}{\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2+{\left(z-z_0\right)}^2}}$ (27)

Rewrite Equation (27) to discrete form, we get the potentials at the mesh nodes:

$\varphi \left(x_i,y_j,z_k\right)=-\frac{h_x{h}_y{h}_z}{4\pi {\epsilon }_0}{\displaystyle \underset{i^{\prime }=1}{\overset{N_x}{\sum }}{\displaystyle \underset{j^{\prime }=1}{\overset{N_y}{\sum }}{\displaystyle \underset{k^{\prime }=1}{\overset{N_z}{\sum }}G\left(x_i-x_{i^{\prime }},y_j-y_{j^{\prime }},z_k-z_{k^{\prime }}\right)\rho \left(x_{i^{\prime }},y_{j^{\prime }},z_{k^{\prime }}\right)}}}$ (28)

5) The electric fields at the mesh nodes are calculated with the central interpolation method:

$E_x=-\frac{{\varphi }_{i+1,j,k}-{\varphi }_{i-1,j,k}}{2h_x}$, $E_y=-\frac{{\varphi }_{i,j+1,k}-{\varphi }_{i,j-1,k}}{2h_y}$, $E_z=-\frac{{\varphi }_{i,j,k+1}-{\varphi }_{i,,j,k-1}}{2h_z}$. (29)

The electric fields at all MP positions are calculated by interpolating from the electric field at the mesh nodes using the same weighting scheme as what for the charge assignment.

6) From the Newton’s law, the contributions of the self-fields to the particle trajectories are:

$\Delta x=\frac{1}{2}q{E}_x\Delta t^2/m$, $\Delta y=\frac{1}{2}q{E}_y\Delta t^2/m$, $\Delta z=\frac{1}{2}q{E}_z\Delta t^2/m$. (30)

where m is the particle mass, Δt = Δz/v₀, Δz is the step length along the reference orbit, v₀ is the velocity of the reference particle. The particle velocity changes are

$\Delta v_x=q{E}_x\Delta t/m$, $\Delta v_y=q{E}_y\Delta t/m$, $\Delta v_z=q{E}_z\Delta t/m$. (31)

Combining the results obtained both from the PIC scheme and the Lie algebraic method, we get the particle trajectories in the optical elements.

4.3. RF Gap Formulas

4.3.1. Single RF Gap

The code provides the beam transport calculations for the RF accelerating structures consisting of QWR or SLR cavities, see Figure 3 and Figure 4.

Figure 3. Quarter wave resonator.

Figure 4. Split loop resonator.

For each RF gap of QWR or SLR the transfer matrix is

$M=M\left(i,j\right),i,j=1,6$ (32)

where

$\begin{array}{l}{M}_{11}=1-\left(\alpha /2\right)\left(T_{k}k/{\gamma }^3+T/\gamma \right)\cos \phi \\ M_{12}=-\left(\alpha /2\right)\left(T_{kk}k/{\gamma }^3\right)\sin \phi \\ M_{21}=-\left(\alpha kT/\left(2{\gamma }^3\right)\right)\sin \phi \\ M_{22}=1-\left(\alpha /\gamma \right)\left(1/{\gamma }^2+{\beta }^2\right)T\cos \phi \end{array}$

$\begin{array}{l}{M}_{33}=M_{11},M_{34}=M_{12}\\ M_{43}=M_{21},M_{44}=M_{22}\\ M_{55}=1+\left(\alpha k{T}_k/{\gamma }^3\right)\cos \phi \\ M_{65}=-QTV\sin \phi \\ M_{66}=1\end{array}$ (33)

The matrix entries not listed above are all zero, and Q is the number of particle charge state, V_m is the voltage amplitude across the RF gap, α = QV/(2E), E the energy of reference particle, k = ω/v, ω is the frequency of the RF field, v is the velocity of the reference particle, T, T_k and T_kk are the time transit factor and its first and second deviations to k, g is the relativistic factor.

4.3.2. Quarter Wave Resonators

To calculate the transfer matrix at each gap the energy of reference particle should be known first. For the QWRs, let f_s be the phase angle of the reference particle at the middle way of a QWR, the energy gain of reference particle at the first gap is:

$\Delta E=Q\cdot V_m\cdot T_{tf}\cdot \cos \left({\varphi }_s+\pi /2-\omega \mathcal{l}/v+{\varphi }_1\right)$ (34)

where T = sin(ωg/2v)/(ωg/2v) is the time transit factor; 2ℓ is the length between the two centers of the two gaps; f₁is the initial phase angle of the reference particle when arriving at the first gap. In the above equation the particle energy depends on its phase angle, and the phase angle is related to the energy (in other word, the velocity v). So, (34) is a transcendental equation and can be solved iteratively. From the gap 1 to gap 2 the particle drift distance 2ℓ, the energy keeps unchanged, but the phase angle increases by 2ℓω/v. The energy gain of reference particle in gap 2 is

$\Delta E=Q\cdot V_m\cdot T_{tf}\cdot \cos \left({\phi }_s-\pi /2+\omega \mathcal{l}/v+{\phi }_2\right)$ (35)

where φ₂ is the phase angle of reference particle at gap 2. At this time the velocity v is a constant, it is not necessary to calculate the energy gain iteratively.

4.3.3. Split Loop Resonator

For the SLR element, at the first gap the phase angle and energy gain of the reference particle can be calculated iteratively according to Equation (34). At the second gap the phase angle is f_s, the energy gain is

$\Delta E=Q\cdot V_m\cdot T_{tf}\cdot \cos {\phi }_s$ (36)

When the particle arrives at the third gap, its phase angle and energy gain can be obtained from Equation (35) directly.

5. Optimization Procedure [12]

We have developed the optimization subroutines which automatically determine the optimal parameters of the optical elements to match a beam from a given initial state to a prescribed final state.

Several optical conditions can be prescribed, such as forming an image, making a beam waist, chromaticity, etc. If the transfer matrix of an element or a part of the beam line is M (i, j), i, j = 1, 6, the beam matrix is σ(i, j), i, j = 1, 6, the following optical conditions could be inserted into the input data file:

image in x-plane: M(1, 2) = 0;

image in y-plane: M(3, 4) = 0;

chromaticity: M(1, 6) = 0, or M(2, 6) = 0;

Parallel—Focusing image: M(1, 1) = 0;

Focusing—Parallel image: M(2, 2) = 0;

Parallel—Parallel image (telescope system): M(2, 1) = 0;

form a waist in x-plane: σ(1, 2) = 0;

beam size limit in x-plane: σ(1, 1) = given value;

form a waist in y-plane: σ(3, 4) = 0;

beam size limit in y-plane: σ(3, 3) = given value;

beam waist in the longitudinal direction: σ(5, 6) = 0.

Powell nonlinear optimization subroutines [4] have been incorporated in the codes. The goal of the optimization calculations is to find out the minimum values of the following object function

$F={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[\left(f_i\left(x_1,x_2,\cdots ,x_m\right)/{\epsilon }_i\right)\right]}^2}=\min$ (37)

where $f_i\left(i=\text{1},\text{2},\cdots ,n\right)$ are the required optical conditions; fin are the given value for this conditions; $x_j\left(j=\text{1},\text{2},\cdots ,m\right)$ are the variable parameters, such as magnetic field, voltage, element length and so on; ε_i is the tolerance for each condition (weight factor).

6. Periodical Structure Calculation [8]

In the linear particle accelerators consisting of QWRs or SLRs, or in some periodically arranged beam lines, the particle beams will pass through these periodic structures. In order to keep the particle motion stable, the program automatically adjusts the magnetic quadrupole fields to fit the following stability condition:

$\left|\cos \mu \right|=\left|0.5Tr\left(M\right)\right|\le 1$, (38)

where μ is the phase shift per period, $M$ is the Twiss matrix shown as the follows:

$M=\left(\begin{array}{cc}\cos \mu +\alpha \sin \mu & \beta \sin \mu \\ -\gamma \sin \mu & \cos \mu -\alpha \sin \mu \end{array}\right)$. (39)

The periodic structures could be combined with the magnetic quadrupoles, drift spaces, as well as QWRs and SLRs.

7. Simulation Examples

In this section we present four example applications of the code. The first example is a low energy beam transport system (LEBT) of the RFQ accelerator. By this example we just want to illustrate the space charge effects to the beam transport. The second one is the optical system of the 400 keV high voltage accelerator. Most of elements of the accelerator are Axisymmetric electrostatic lenses. The third one is a particle distribution uniformization system, which exhibits the applications of multipoles. The last one is an RF linear accelerator system combined with QWRs and quadrupoles.

7.1. The LEBT System of an RFQ Accelerator

With the code, we calculated the LEBT system which delivers the D ⁺ beam to the RFQ accelerating cavity (see Figure 5).

Figure 5. Layout of the LEBT.

The LEBT system consists of two solenoid lenses and some drift spaces. The radius of the solenoids are 50 mm, the lengths of them are all 200 mm, the longitudinal magnetic fields in the two solenoids are all 4.77 kG.

The D ⁺ beam energy is 50 keV, the average bunched beam current is 5 mA, the initial parameters of the beam in the phase space are 1 ´ 40 mm-mrad both in the x- and y-directions. The pulse repetition is 9 MHz, and the initial longitudinal phase space (δE× δf) is 1 keV × 60˚. The calculated beam envelopes of the system are shown in Figure 6.

We can see from Figure 6 that the beam envelopes of 0 mA beam current are much different from the envelopes of 5 mA beam current.

Figure 6. Beam envelopes in the LEBT system. (a) Beam envelopes with 0 mA current; (b) Beam envelopes with 5 mA current.

7.2. The 400 keV High Voltage Accelerator

The second example is an optical system of the 400 keV high voltage accelerator (Figure 7). The initial proton beam parameters are: x = y = ±1 mm, x¢ = y¢ = ±40 mrad, δE= ±150 eV, δf= ±60˚. The Einzel lens 1 focuses the beam to the center of the gap lens 2 which matches the beam to the high voltage accelerating columns. Usually, the beam emerging from the accelerating columns is divergent, so a quadrupole doublet 4 focuses the beam again to the target. Figure 8 shows the beam envelopes in the optical system.

Figure 7. Layout of the 400 keV high voltage accelerator.

Figure 8. Beam envelopes in the 400 keV high voltage accelerator.

7.3. Particle Uniform Distribution System

The third example of the applications of the code is a particle uniform distribution system. It is made up of two quadrupole triplets, one quadrupole doublet and two octupoles as shown in Figure 9. The first octupole locates at a beam waist in x-plane; another octupole locates at the beam peak.

Figure 9. Particle distribution uniformization system.

The initial proton beam energy is 400 KeV. The initial particle distribution in the x-y transverse plane is Gaussian. The initial beam parameters in the 6D phase space are (x, x¢) = (±1 mm, ±20 mrad), (y, y¢) = (±1 mm, ±20 mrad), (δE, δf) = (±150 eV, ±60˚). The beam envelopes are shown in Figure 10, the initial and final transverse distributions are shown in Figure 11.

7.4. The RF Linear Accelerator

The fourth example is a proton RF linear accelerator consisting of some QWRs. The accelerator consists of an array of 24 QWRs distributed in 6 cryostat modules (Figure 12). Each module contains 4 QWRs of the same structure and same physical parameters. The designed β of the linac is 0.1945, synchrotron phase

Figure 10. Beam envelopes in the uniform distribution system.

Figure 11. Particle distributions in transverse plane (x, y). (a) Initial, Gaussian distribution; (b) Final, approximately uniform and square.

1, 2, 3, 4: QWRs, 5: Quadrupole doublet

Figure 12. Structure of a module.

angle −26.0˚. The accelerating electric field of each QWR is 3.0 MV/m, working frequency 150 MHz. Figure 13 is the layout of the linac.

1, 8: Quadrupole triplets; 2, 3, 4, 5, 6, 7: Modules

Figure 13. Layout of the linac.

The initial proton beam energy is 18.13 MeV. The initial beam parameters in the 6D phase space are (x, x¢) = (±1 mm, ±5 mrad), (y, y¢) = (±1 mm, ±5 mrad), (δE, δf) = (±300 eV, ±60˚). The beam envelopes in the linac are shown in Figure 14.

Figure 14. Beam envelopes in the QWR linac.

8. Conclusion

A code package LEADS-v5 for nonlinear beam transport study with space charge effects has been developed. The code calculates high voltage accelerators which consist of axisymmetric electrostatic lenses and accelerating columns. Beam transport lines and RF linacs of QWR/SLR structures can be simulated also. The code provides multipole magnets which allow the user to study the particle distribution uniformization systems. The optimization procedures provide the possibility to obtain the prescribed optical condition automatically. Many calculated examples show that the code is reliable.

Founding

Work supported by the National Natural Science Foundation of China, Contract No. 1057009, and by SRFDP (20070001001).

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Dragt, A.J. (1982) Lectures on Nonlinear Orbit Dynamics. In: Carrigan, R.A., et al., Eds., AIP Conference Proceedings No. 87, American Institute of Physics, New York, 147-269.
[2]	Dragt, A.J. and Finn, J.M. (1976) Lie Series and Invariant Functions for Analytic Symplectic Maps. Journal of Mathematical Physics, 17, 2215. https://doi.org/10.1063/1.522868
[3]	Buneman, O. (1959) Dissipation of Currents in Ionized Media. Physical Review, 115, 503-517. https://doi.org/10.1103/PhysRev.115.503
[4]	Dawson, J.M. (1962) One-Dimensional Plasma Model. The Physics of Fluids, 5, 445-459. https://doi.org/10.1063/1.1706638
[5]	Birdsall, C.K. (1997) Clouds-in-Clouds, Clouds-in-Cells Physics for Many-Body Plasma Simulation. Journal of Computational Physics, 135, 141-148. https://doi.org/10.1006/jcph.1997.5723
[6]	Hockney, R.W. and Eastwood, J.W. (1987) Computer Simulation Using Particles. CRC Press, Boca Raton. https://doi.org/10.1201/9781439822050
[7]	Brown, K.L., Rothacker, F., Carey, D.C. and Iselin, Ch. (1974) Report SLAC-91.
[8]	Adams, A. and Read, F.H. (1972) Electrostatic Cylinder Lenses II: Three Element Einzel Lenses. Journal of Physics E: Scientific Instruments, 5, 150. https://doi.org/10.1088/0022-3735/5/2/019
[9]	Grivet, P. (1965) Electron Optics. Pergamon, New York.
[10]	Harting, E. and Read, F.H. (1976) Electrostatic Lenses. Elsevier, Amsterdam, Oxford, New York.
[11]	Courant E.D. and Snyder M.S. (1976) Electrostatic Lenses. Elsevier Scientific Pub. Co., Amsterdam, New York.
[12]	Himmelblau, D.M. (1972) Applied Nonlinear Programming. McGraw-Hill, New York.