HOME work papers html 2009fermi

4.3 Detailled description of the method

4.3.1 Frame definitions

**Figure 4.3:** Frame definition - laboratory - rest frame - bunch frame

The laboratory frame $\mathcal{R}$ , the rest frame $\mathcal{R}'$ . We will note respectivly

and

the canonical bases adapted to the laboratory frame, the rest frame and the main directions of the bunch.

4.3.2 Boosting the distribution to the rest frame

The Lorentz space transformation that should be applied to the particle has been studied in a previous section. For the numerical situation, the time dependence is removed so that, the matrix of equation 1.15 reduces to:

$\displaystyle A(\v{\beta}) =\left( \begin{array}{cccc} 1+ (\gamma-1)\frac{\beta... ...beta_z}{\beta^2} & 1+ (\gamma-1)\frac{\beta_z^2}{\beta^2}\\ \end{array} \right)$

(4.1)

Coordinates of a particle in the rest frame $\mathcal{R}'$ are obtained with the following multiplication:

$\displaystyle X\vert _{e'} = A \cdot X\vert _e$

(4.2)

4.3.3 Rotation of the bunch to its principal axis

To adapt the meshing to the distribution, one need to detect the distribution main axis. For this, we used the covariance matrix of the particles. Let being this 3 $\times$ 3-matrix. For this method, there should be no difference between using coordinates in the rest frame or in the laboratory frame. In our case, we chose the coordinates in the rest frame $\mathcal{R}'$ .

$\displaystyle C =$ Cov $\displaystyle \left\{x'_i , y'_i , z'_i \right\}$

(4.3)

is by definition a symmetric real matrix, hence diagonalizable to a matrix

. Once we find the eigen vectors of this matrix, we have the principal axis of the distribution. Thus diagonalizing the covariance matrix is an elegant method to further adapt the meshing on the distribution. The diagonalization process is described in the following section 4.3.4 and in annex B. Rotations and base transformation can be confusing, but not if one takes the time to write the proper formalism:

$\displaystyle P\vert _{e'}^{e'_b}$	$\displaystyle =$	Mat $\displaystyle \left( \text{Id}, e'_b, e' \right)$	(4.4)
$\displaystyle P\vert _{e'_b}^{e'}$	$\displaystyle =$	Mat $\displaystyle \left( \text{Id}, e', e'_b \right) = {P\vert _{e'}^{e'_b}}^{-1}$	(4.5)

Where $P\vert _{e'}^{e'_b}$ is the transformation matrix from the base

to the base

. We will immediatly note that these matrix are rotations, i.e. orthogonal matrix, such that $P^{-1}=P^{T}$ and $\det{P}=1$ . The diagonalization of

can be writen:

$\displaystyle D\vert _{e'_b} = P\vert _{e'_b}^{e'} \cdot C\vert _{e'}^{e'} \cdot P\vert _{e'}^{e'_b}$

(4.6)

Let us define $E'_{b,1}\vert _{e'}, E'_{b,2}\vert _{e'}, E'_{b,3}\vert _{e'}$ the three normalised eigen vectors of forming the base , which coordinates are expressed in the base . The 3 $\times$ 3-matrix formed by the three column vectors $E'_{b,i}$ is:

$\displaystyle P\vert _{e'}^{e'_b}=\left (E'_{b,1} \ E'_{b,2} \ E'_{b,3} \right)$

(4.7)

And to simplify notations we will further write it : According to the previous formalism we then have:

$\displaystyle R$	$\displaystyle =$	$\displaystyle P\vert _{e'}^{e'_b}$	(4.8)
$\displaystyle D$	$\displaystyle =$	$\displaystyle R^{T}\cdot C \cdot R$	(4.9)

And eventually, let's write down the fact that, for a given vector

expressed in the rest base

, its coordinaates in

are:

$\displaystyle X\vert _{e'_b} = P\vert _{e'_b}^{e'} \cdot X\vert _{e'} = R^{T} \cdot X\vert _{e'}$

(4.10)

Now that this algebra is clear, the implementation of the algorithm is straightforward

Compute the covariance matrix
Diagonalize it and store the eigenvectors as the rotation matrix
Multiply the coordinates of the particles in the rest frame by the transposed rotation matrix
Go on with the 3D poisson solver in this bunch-adapted base

4.3.4 Diagonalization algorithm

A choice between the several algorithms of diagonalization has been made. From our analysis, the Jacobi algorithm has been selected as it takes into account the fact that the matrix is symmetric. We developped in annex B the complete Jacobi method, and provide as well the source code I wrote in fortran to implement this algorithm. The other methods considered were:

The Givens reduction: which is a variant of the Jacobi method, which requires a finite number of iteration and return a tridiagonal matrix
The Householder algorithm: reduces a symmetric matrix A to tridiagonal form by orthogonal transformations.
The factorization method: another common technique consisting in decomposing the matrix into a left factor and a right factor. This method, is more subtle and does not exploit the symetric form of the matrix.

The Jacobi method has been selected for its simplicity of principle, and its rapidity to diagonalize a 3 $\times$ 3 symmetric matrix.

4.3.5 Space charge calculation

4.3.5.0.1 The 2D space charge algortihm of Astra

is a cylindrical symmetric algorithm. Each grid consist in a cylindrical ring in the radial direction forming a slice in the longitudinal direction. After ``boosting'' the grid in the rest frame, a static field calculation is performed, with the assumption that the charge is constant in each ring. Fields are calculated at the center of each grid and transformed back in the laboratory frame. To determine the fields in arbitrary location, a cubic spline interpolation is used. Outside of the grid, an

extrapolation is used.

4.3.6 The 3d space charge algorithm

follow the resolution method of Poisson's equation developped in chapter 5 of this document. Its main hypothesis is that the particle are at rest after boost, so that Poisson equation is valid. Indeed if the particle are at rest, no courant is present in the distribution and thus Maxwell equations reduces to Poisson equation. This also means that the magnetic field is null in the rest frame $\v{B'} = 0$ .

4.3.7 Solving the equations of motion

Astra uses a Runge-Kutta model to perfom the integration of the equations of motion. In these equations, both the contributions of the external and internal fields are taken into account. Nevertheless, internal forces being expected to evolve slowly, the space charge algorithm is not executed at every time step. Different criteria, such as the evolution of the aspect ratio of the bunch, determine whether the space charge algorithm should be called. If not, the previously calculated fields are scaled accordingly to the variation of the beam size.

Emmanuel Branlard
http://emmanuel.branlard.free.fr/