Linear Transverse Motion

In this chapter,we will discuss the particle motion in perpendicular to the particle ideal trajectory in linear magnets. Linear magnets refer to the dipole and quadrupole magnets that produce linear magnetic field. The force of the particle is given by the Lorentz force:

\begin{equation} \mathbf{F}=q\mathbf{v}\times \mathbf{B} \end{equation}

The force is always perpendicular to its velocity, therefore the particle can not be accelerated by the magnetic field.

Design Trajectory

The geometry of the accelerator is determined by the dipole magnet. Let us assume that the magnetic field is constant in certain region and along the vertical direction $\mathbf{B}=B\hat{y}$. The particle travels through this region initially in $\hat{s}$ direction in the $x$-$s$ plane, as shown below.

The equation of motion is

\begin{equation} \frac{d\mathbf{P}}{dt}=q\mathbf{v}\times\mathbf{B} \end{equation}

The momentum change will be only in the $x$-$s$ plane, the velocity is always perpendicular to the magnetic field and the energy of the beam does not change, therefore:

\begin{equation} m\gamma \frac{d\mathbf{v}}{dt}=qB\mathbf{v}\times\hat{y} \end{equation}

Using the fact that the amplitude of the velocity will not change therefore:

\begin{equation} m\gamma \frac{v^2}{\rho}=qvB \end{equation}

The bending radius is given by:

\begin{equation} \rho=\frac{\left|\mathbf{P}\right|}{qB} \end{equation}

Here, we define the term 'rigidity' $B\rho=P/q$, to quantify 'how hard to bend the beam'. And naturally it is only function of the particle momentum amplitude. It is a very important normalization factor of the transverse motion, as will shown later. After normalization, most our calculation becomes particle energy independent.

The dipoles that defines the geometry of the accelerator are sometimes called main dipoles. They define a ideal trajectory. A charge particle with correct energy, once inject with correct position and angle, will exactly follow the ideal trajectory.

First look: Beam motion in Quadrupole field

The quadrupole field, from the previous section, has

\begin{align} B_x&=Gy\\ B_y&=Gx \end{align}

where the gradient is $G=B_0 b_1$.

Quadrupole is usually positioned so that the zero field axis is located on the ideal trajectory. The ideal particle will not experience any field from the Quadrupole.

From the Lorentz force, we have $$\gamma m \frac{d\mathbf{v}}{dt} = q \mathbf{v}\times \mathbf{B}$$

Expand the vector production and get

\begin{align} \gamma m\frac{d^{2}x}{dt^{2}}&=-qGxv_{z}\\ \gamma m\frac{d^{2}y}{dt^{2}}&=qGyv_{z} \end{align}

Here we need to make following approximation:

  • paraxial approximation, indicate that the particle will only deviate the design direction by small angles, i.e. if z is the forward direction, the two transverse velocity satisfy: $v_x,v_y \ll \left|\mathbf{v}\right|$
  • In such approximation,
\begin{align} v_z&=\sqrt{v^2-v_x^2-v_y^2}\nonumber\\ &=v\left(1-\frac{v_x^2+v_y^2}{2v^2}\right)\sim v \end{align}
  • Using such approximation, we have $$\frac{d}{dt}=v_s\frac{d}{ds}\sim v\frac{d}{ds}\;\frac{d^2}{dt^2}\sim v^2\frac{d^2}{ds^2}$$

Then the equation of motion in quadrupole is simplified as

\begin{align} \gamma m v\frac{d^{2}x}{ds^{2}}&=-qGx\\ \gamma m v\frac{d^{2}y}{ds^{2}}&=qGy \end{align}

Using the definition of the rigidity $B\rho$, the EOM can be rewritten as :

\begin{align} \frac{d^{2}x}{ds^{2}}&=-\frac{G}{B\rho}x\\ \frac{d^{2}y}{ds^{2}}&=\frac{G}{B\rho}y \end{align}

Then, we define the normalized quadrupole strength $k=G/(B\rho)$, the EOM is simply:

\begin{align} \frac{d^{2}x}{ds^{2}}+kx&=0\\ \frac{d^{2}y}{ds^{2}}-ky&=0 \end{align}

Depend on the sign of $k$, it is easy to solve the differential equation:

\begin{equation} x\left(s\right)= \begin{cases} a \cos \left(\sqrt{k}s\right)+b \sin \left(\sqrt{k}s\right)&k>0\\ a s + b &k=0\\ a \cosh \left(\sqrt{-k}s\right)+b\sinh \left(\sqrt{-k}s\right)&k<0 \end{cases} \end{equation}

The three cases corresponds to a focusing quadrupole, drift space and a defocusing quadrupole respectively. The parameter $a$, $b$ are determined by the initial particle position $x$ and its derivative $x'$. We also learn that the quadrupole will always focus the particles in one transverse direction and defocus in the other.

In general, the focusing strength can be also function of the longitudinal direction and the EOM of transverse motion (using horizontal direction $\hat{x}$ as an example) gives:

\begin{equation} \frac{d^{2}x}{ds^{2}}+k\left(s\right)x=0 \end{equation}

This ODE is named Hill's equation. It can be rewritten as:

\begin{align} \frac{dx}{ds}&=x'\\ \frac{dx'}{ds}&=-k(s)x \end{align}

They are just the form of Hamiltonian equations of the Hamiltonian:

\begin{align} H=\frac{1}{2}x'^2+\frac{1}{2}k(s)x^2 \end{align}

Latter we will use alternative method to rederive this equation and try to solve it in general.

Curvilinear Coordinate System

The ideal trajectory is usually used as reference for the dynamics of the charged particles. This trajectory can be determined from a vector $\mathbf{r_0}\left(t\right)$. The length along the curve is given by

\begin{align} s(t)=\int_0^t \left|\mathbf{r'_0}\left(t^\prime\right)\right| dt^\prime \end{align}

Therefore our study is usually performed in the curvilinear coordinate system, named Frenet-Serret coordinate system, as illustrated in the figure:

In this representation, the independent variable is the location on the reference orbit $s$. The reference point on the reference orbit is defined by the vector $\mathbf{r_0}\left(s\right)$. The unit vector of the tangent direction is given by:

\begin{equation} \hat{s}\left(s\right)=\frac{d\mathbf{r_0}\left(s\right)}{ds} \end{equation}

The unit vector of the tangent plane and the normal plane is:

\begin{equation}\hat{x}\left(s\right)=-\rho\left(s\right)\frac{d\hat{s}\left(s\right)}{ds} \end{equation}

Here $\rho$ denotes the local radius, which is the recipocal of the curvature $\kappa$. At last, the binormal vector is simply defined by:

\begin{equation} \hat{y}\left(s\right)=\hat{s}\left(s\right)\times\hat{x}\left(s\right) \end{equation}

We defining The derivative of the unit vectors are:

\begin{align} \frac{d\hat{s}\left(s\right)}{ds} &= -\kappa\left(s\right) \hat{x}\left(s\right)\\ \frac{d\hat{x}\left(s\right)}{ds} &= \kappa\left(s\right) \hat{s}\left(s\right)+\tau\left(s\right)\hat{y}\left(s\right)\\ \frac{d\hat{y}\left(s\right)}{ds} &= -\tau\left(s\right)\hat{x}\left(s\right) \end{align}

Here, $\tau$ is the local torsion of the curve. In matrix form, the relation can be expressed as:

\begin{equation} \frac{{d}}{ds}\left(\begin{array}{c} \hat{{s}}\\ \hat{{x}}\\ \hat{y} \end{array}\right)=\left(\begin{array}{ccc} 0 & -\kappa & 0\\ \kappa & 0 & \tau\\ 0 & -\tau & 0 \end{array}\right)\left(\begin{array}{c} \hat{{s}}\\ \hat{{x}}\\ \hat{y} \end{array}\right) \end{equation}

In most cases, we focus our discussion on the motion in plane geometry, therefore the torsion is zero.

Then, the coordinate of particle can be expressed as:

\begin{equation} \mathbf{r}=\mathbf{r_0}\left(s\right)+x\hat{x}\left(s\right)+y\hat{y}\left(s\right) \end{equation}

Compared with the Cartesian coordinate system, $\mathbf{r}=x \hat{x}+y \hat{y}+z \hat{z} $, the differentials give: $\partial \mathbf{r}/\partial x = \hat{x}$, $\partial \mathbf{r}/\partial y = \hat{y}$, $\partial \mathbf{r}/\partial z = \hat{z}$; the differentials of the Frenet-Serret Coordinate yields:

\begin{align} \frac{\partial \mathbf{r}}{\partial s} &= \hat{s}\left(1+ \kappa x\right) \\ \frac{\partial \mathbf{r}}{\partial x} &= \hat{x} \\ \frac{\partial \mathbf{r}}{\partial y} &= \hat{y} \end{align}

The Lame coefficient of the Frenet-Serret coordinate system is $h_s=1+\kappa x$ and $h_x=h_y=1$.

In this coordinate system, the vector operation has different form, for example:

\begin{equation} \nabla \phi = \frac{1}{1+\kappa x}\frac{\partial \phi}{\partial s}\hat{s}+\frac{\partial \phi}{\partial x}\hat{x}+\frac{\partial \phi}{\partial y}\hat{y} \end{equation}

And the curl of the vector gives:

\begin{equation} \nabla\times A=\frac{1}{1+\kappa x}\left|\begin{array}{ccc} \left(1+\kappa x\right)\hat{{s}} & \hat{{x}} & \hat{y}\\ \frac{\partial}{\partial s} & \frac{\partial}{\partial x} & \frac{\partial}{\partial y}\\ \left(1+\kappa x\right)A_{s} & A_{x} & A_{y} \end{array}\right| \end{equation}

Hamiltonian Approach

Lagrangian to Hamiltonian

We consider the approach from Lagrangian of a relativistic charged particle in a E-M field.

\begin{equation} L\left(\mathbf{x}, \mathbf{v}, t\right)= -mc^2\sqrt{1-\frac{v^2}{c^2}} +e\mathbf{v}\cdot \mathbf{A}\left(\mathbf{x},t\right) -e\phi \left(\mathbf{x},t\right) \end{equation}

Here, $\mathbf{x}$ is the position and the velocity is $\mathbf{v}=\dot{\mathbf{x}} =d\mathbf{x}/dt$. The vector potential $\mathbf{A}$ and scalar potential $\phi$ is the function of the position and time. The fields is calculated as:

\begin{align} \mathbf{E}&=-\nabla\phi-\frac{\partial\mathbf{A}}{\partial t}\\ \mathbf{B}&=\nabla\times\mathbf{A} \end{align}

The equation of motion can be derived from the equation

\begin{equation} \frac{\partial L}{x_i}=\frac{d}{dt}\frac{\partial L}{\dot{x}_i} \end{equation}

Using the Lagrangian above will get the EOM with Lorentz force.

We are interested in the generalized momentum $P_i$ of the $i^{th}$ coordinate $x_i$ which is given by:

\begin{align} P_i&=\frac{\partial L}{\partial \dot{x}_i} \nonumber\\ &=\frac{m\dot{x}_i}{\sqrt{1-\frac{v^2}{c^2}}}+eA_i \nonumber\\ &=p_i+eA_i \end{align}

Here $p_i$ is the momentum of the particle. The above components can be combined to get a vector form, $\mathbf{P}=\mathbf{p}+e\mathbf{A}$, between the generalized (or conjugate) momentum $\mathbf{P}$, the momentum $\mathbf{p}$ and the vector potential $\mathbf{A}$.

Then the Hamiltonian is calculated as:

\begin{align} H&=\mathbf{v}\cdot\mathbf{P}-L \nonumber\\ &=\mathbf{v}\cdot\mathbf{p}+e\mathbf{v}\cdot\mathbf{A}+mc^2\sqrt{1-\frac{v^2}{c^2}} -e\mathbf{v}\cdot \mathbf{A}\left(\mathbf{x},t\right) +e\phi \left(\mathbf{x},t\right) \nonumber\\ &=\gamma m v^2 +mc^2\sqrt{1-\frac{v^2}{c^2}}+e\phi \left(\mathbf{x},t\right) \nonumber\\ &=(\gamma\beta^2+1/\gamma)mc^2+e\phi \left(\mathbf{x},t\right) \nonumber\\ &=\gamma mc^2 +e\phi \left(\mathbf{x},t\right) \end{align}

With no surprise, the Hamiltonian reflects the total energy of the particle and the potential energy in the E-M field. For further analysis, we express the Hamiltonian as function of the conjugate variable pair:

\begin{align} H(\mathbf{x},\mathbf{p})&=H\left(x_1,\cdots,x_n, P_1,\cdots,P_n, t\right)\\ &=\sqrt{\left(m c^2\right)^2 +\left(\mathbf{P}-e\mathbf{A}\left(\mathbf{x},t\right)\right)^2c^2} +e\phi\left(\mathbf{x},t\right) \end{align}

Later we will use a series of canonical transformations to change the coordinates and independent variable. A generating function is handy to create the canonical transformations from old coordinate $(q,p)$ to a new set (Q,P). There are four types of generating functions:

\begin{align} F_1=F_1(\mathbf{q},\mathbf{Q},t) \quad & \mathbf{p}=\frac{\partial F_1}{\partial \mathbf{q}}; \mathbf{P}=-\frac{\partial F_1}{\partial \mathbf{Q}}\\ F_2=F_2(\mathbf{q},\mathbf{P},t) \quad & \mathbf{p}=\frac{\partial F_2}{\partial \mathbf{q}}; \mathbf{Q}=\frac{\partial F_2}{\partial \mathbf{P}}\\ F_3=F_3(\mathbf{p},\mathbf{Q},t) \quad & \mathbf{q}=-\frac{\partial F_3}{\partial \mathbf{p}}; \mathbf{P}=-\frac{\partial F_3}{\partial \mathbf{Q}}\\ F_4=F_4(\mathbf{p},\mathbf{P},t) \quad & \mathbf{q}=-\frac{\partial F_4}{\partial \mathbf{P}}; \mathbf{Q}=\frac{\partial F_4}{\partial \mathbf{P}} \end{align}

The right figure shows a way to remember the relative relations of the generating functions. The new Hamiltonian $K(Q,P,t)$ is calculated by:

\begin{align} K(\mathbf{Q},\mathbf{P},t)=H\left(\mathbf{q}(\mathbf{Q},\mathbf{P}),\mathbf{p}(\mathbf{Q},\mathbf{P})\right)+\frac{\partial F}{\partial t} \end{align}

Transformation to Frenet-Serret coordinate system

The first necessary transformation that we need to make is to change the coordinate system to the Frenet-Serret coordinate system, i.e. from Cartesian coordinate $\mathbf{r}=x\hat{x}+y\hat{y}+z\hat{z}$ to $\mathbf{r}=\mathbf{r}_0(s)+x\hat{x}(s)+y\hat{y}(s)$.
We use the generating function of the $3^{rd}$ type to perform this transformation:

\begin{equation} F_3=-\mathbf{P}\cdot\left(\mathbf{r}_0(s)+x\hat{x}(s)+y\hat{y}(s)\right) \end{equation}

The new conjugate momentum, $\boldsymbol{\Pi}$ can be calculated as:

\begin{align} \Pi_x&=-\frac{\partial F_3}{\partial x}=\mathbf{P}\cdot\hat{x}\\ \Pi_y&=-\frac{\partial F_3}{\partial y}=\mathbf{P}\cdot\hat{y}\\ \Pi_s&=-\frac{\partial F_3}{\partial s}=\mathbf{P}\cdot\hat{s}\left(1+\frac{x}{\rho}\right)\\ \end{align}

Therefore, in the Frenet-Serret coordinate system:

\begin{equation} H=c\sqrt{ \left(m c\right)^2 +\left(\Pi_x-eA_x\right)^2+\left(\Pi_y-eA_y\right)^2 +\left(\frac{\Pi_s}{1+x/\rho}-eA_s\right)^2 } +e\phi \end{equation}

The vector potential now also take its components from the Frenet-Serret coordinate system as $A_x=\mathbf{A}\cdot\hat{x}$, $A_y=\mathbf{A}\cdot\hat{y}$ and $A_s=\mathbf{A}\cdot\hat{s}$. Sometimes, the longitudinal component of the vector potential also can be define as $A_s=\mathbf{A}\cdot \frac{\partial \mathbf{r}}{\partial s}= (1+x/\rho) \mathbf{A}\cdot\hat{s}$, which will NOT be used in the latter contents.

Use $s$ as independent variable

Next, we need charge the independent coordinate from time to '$s$', the position on the reference orbit. The Hamiltonian above is function of $(x,\Pi_x, y, \Pi_y, s, \Pi_s,t)$ and the Hamiltonian equation are:

\begin{align} \frac{\partial H}{\partial x}&=-\dot{\Pi}_x & \frac{\partial H}{\partial \Pi_x}&=\dot{x} \\ \frac{\partial H}{\partial y}&=-\dot{\Pi}_y & \frac{\partial H}{\partial \Pi_y}&=\dot{y} \\ \frac{\partial H}{\partial s}&=-\dot{\Pi}_s & \frac{\partial H}{\partial \Pi_s}&=\dot{s} \end{align}

To change the independent coordinate to $s$, we first look at the derivative:

\begin{equation} x^\prime\equiv\frac{dx}{ds}=\frac{dx}{dt}\left(\frac{ds}{dt}\right)^{-1} =\frac{\partial H}{\partial \Pi_x}\left(\frac{\partial H}{\partial \Pi_s}\right)^{-1} \end{equation}

Using the differential chain rule of the implicit equation $H(x,\Pi_x, y, \Pi_y, s, \Pi_s,t)-f(t)=0$, we have:

\begin{equation} x^\prime\equiv\frac{dx}{ds}=-\frac{\partial \Pi_s}{\partial \Pi_x} \end{equation}

Use the similar trick, we will find that

\begin{align} {\Pi_x}^\prime&\equiv\frac{d\Pi_x}{ds}=\frac{\partial \Pi_s}{\partial x}\\ y^\prime&\equiv\frac{dy}{ds}=-\frac{\partial \Pi_s}{\partial \Pi_y}\\ {\Pi_y}^\prime&\equiv\frac{d\Pi_y}{ds}=\frac{\partial \Pi_s}{\partial y}\\ \end{align}

We also find the other two relationship easily:

\begin{equation} t^\prime\equiv\frac{dt}{ds}=\left(\frac{\partial H}{\partial \Pi_s}\right)^{-1}=\frac{\partial \Pi_s}{\partial H} \end{equation}\begin{equation} H^\prime\equiv\frac{dH}{ds}=-\frac{\partial \Pi_s}{\partial t} \end{equation}

Therefore we can choose the new Hamiltonian $K=-\Pi_s(x,\Pi_x, y, \Pi_y, H, t,s)$, which will satisfy the Hamiltonian relation. However, it is awkward to set $H$ as the 'coordinate' and $t$ as 'general momentum'. To fix this, we can set the conjugate pair as $(t,-H)$. Therefore the new Hamiltonian reads:

\begin{align} K\left(x,\Pi_x, y, \Pi_y, t, -H,s\right) = -\Pi_s \nonumber\\ =-\left(1+\frac{x}{\rho}\right)\sqrt{\left(\frac{H-e\phi}{c}\right)^2-(mc)^2-\left(\Pi_x-eA_x\right)^2-\left(\Pi_y-eA_y\right)^2}\nonumber\\-eA_s\left(1+\frac{x}{\rho}\right) \end{align}

Simplifications of Hamiltonian

Since the new Hamiltonian has a very complicated form, we will use the following simplifications.

  • First, we have seen that for 2-D magnetic field, only the longitudinal vector potential is needed. We can set the both transverse potential to be zero: $A_x=A_y=0$. Then, the generalized momentum reduce to the particle's momentum $\Pi_{x/y}=p_{x/y}$.

  • Second, let us assume the potential $\phi$ is not time dependent, then we can make canonical transformation $(t, -H)$ to $(t, -h)$ with $h=H-e\phi$. $h$ is just the total energy of the particle and the first two terms give the momentum of the particle $p^2=\left(h/c\right)^2-(mc)^2$.

  • Third, Let's define a reference momentum of the particle $p_0$, and each particle will have a deviation from the reference, $p=p_0(1+\delta)$. $\delta$ is called fractional momentum deviation. We scale the Hamiltonian $K$ by factor of $p_0$ and get new $\tilde{K}$.

\begin{equation} \frac{K}{p_0} \equiv\tilde{K}\left(x, \tilde{p_x}=\frac{p_x}{p_0}, y, \tilde{p_y}=\frac{p_y}{p_0}, t,-h=\frac{-H}{p_0}\right)=\\ -\left(1+\frac{x}{\rho}\right) \sqrt{\left(\frac{h}{c}\right)^2-\left(\frac{mc^2}{p_0}\right)^2-\tilde{p_x}^2-\tilde{p_y}^2} -e\frac{A_s}{p_0}\left(1+\frac{x}{\rho}\right) \end{equation}
  • Last, we will make a canonical transformation from the pair $(t, -h)$ to $(z,\delta)$ with generating function $F_2=f(\delta)t$ and
\begin{equation} -h=\frac{\partial F_2}{\partial t}=f(\delta) \end{equation}

The new coordinate $z$ is calculated as

\begin{align} z&=\frac{\partial F_2}{\partial \delta} \nonumber \\ &=\left(-\sqrt{c^2(1+\delta)^2+m^2c^4}\right)'t\nonumber \\ &=-\beta ct \end{align}

the Hamiltonian becomes:

\begin{equation} \tilde{K}\left(x, \tilde{p_x}, y, \tilde{p_y}, z=-\beta ct, \delta,s\right)=-\left(1+\frac{x}{\rho}\right)\sqrt{(1+\delta)^2-\tilde{p_x}^2-\tilde{p_y}^2} -e\frac{A_s}{p_0}\left(1+\frac{x}{\rho}\right) \end{equation}

The new longitudinal coordinate $z$ has the meaning: the distance of the particle travels at current speed.

For later convenience, we remove the tilde of the above Hamiltonian and it looks like:

\begin{equation} \tilde{K}\left(x, \tilde{p_x}, y, \tilde{p_y}, z=-\beta ct, \delta,s\right)=-\left(1+\frac{x}{\rho}\right)\sqrt{(1+\delta)^2-\tilde{p_x}^2-\tilde{p_y}^2} -e\frac{A_s}{p_0}\left(1+\frac{x}{\rho}\right)\label{eq:Hamiltonian_popular_form} \end{equation}

Further assumptions

The above Hamiltonian can be consider as exact in most accelerator cases. Next, we will make assumptions and expand the Hamiltonian up to 2nd order. The assumptions includes:

  • Momentum spread is very small. For any particle, $\delta \ll 1$
  • Paraxial approximation, the longitudinal momentum is close to the reference momentum and $p_x,p_y \ll 1$

Then we can expand the Hamiltonian up to 2nd order in $p_x$ and $p_y$ and have:

\begin{align} K\left(x, p_x, y, p_y, z, \delta,s\right)&=-(1+\delta)\left(1+\frac{x}{\rho}\right)\left(1-\frac{p_x^2}{2(1+\delta)^2}-\frac{p_y^2}{2(1+\delta)^2}\right) -e\frac{A_s}{p_0}\left(1+\frac{x}{\rho}\right) \nonumber\\ &=-\left(1+\frac{x}{\rho}\right)\left(1+\delta-\frac{p_x^2}{2(1+\delta)}-\frac{p_y^2}{2(1+\delta)}\right) -e\frac{A_s}{p_0}\left(1+\frac{x}{\rho}\right) \end{align}

Recall that $$p_x=\frac{\mathbf{p}\cdot\hat{x}}{p_0} \;;\; p_y=\frac{\mathbf{p}\cdot\hat{y}}{p_0}$$

is different from the orbit slope:

\begin{equation} x^\prime\equiv\frac{dx}{ds}=\frac{\mathbf{p}\cdot\hat{x}}{p_s}\sim\frac{\mathbf{p}\cdot\hat{x}}{p_0(1+\delta)}=\frac{p_x}{1+\delta} \end{equation}\begin{equation} y^\prime\equiv\frac{dy}{ds}=\frac{\mathbf{p}\cdot\hat{y}}{p_s}\sim\frac{\mathbf{p}\cdot\hat{y}}{p_0(1+\delta)}=\frac{p_y}{1+\delta} \end{equation}

Therefore, we only consider the transverse motion without energy dependence, the slope $x^\prime$ and $y^\prime$ can be used as conjugate momentum in paraxial approximation. Otherwise using the slope instead of $p_x$ and $p_y$ will violate Hamiltonian dynamics.

In 4-Dimensional motion for particle with reference momentum, the Hamiltonian is

\begin{equation} K\left(x, p_x, y, p_y,s\right)=-\left(1+\frac{x}{\rho}\right)\left(1-\frac{p_x^2}{2}-\frac{p_y^2}{2}\right) -e\frac{A_s}{p_0}\left(1+\frac{x}{\rho}\right) \end{equation}

Vector potentials for linear elements

Now, its time to consider the real form for the vector potential $A_s$. We will limit our discussion within linear regime.

For a main dipole, where the reference orbit bend inside with radius $\rho$ or its reciprocal, curvature $\kappa=1/\rho$, we need to find $A_s$ to get vertical magnetic field $\mathbf{B}=B_0(s)\hat{y}$ in the Frenet-Serret coordinate system:

\begin{equation} B_0(s)=-\frac{1}{1+\kappa x}\frac{\partial\left((1+\kappa x)A_s\right)}{\partial x} \end{equation}

It seems exact solution for this differential equation doesnot exist. However, with the assumption of $\left|x\right|\kappa \ll 1$, we can get approximation of $A_s$ up to second order. Let's assume $A_s=-B_0(s)(a x + bx^2)$. Then we get:

\begin{equation} 1+\kappa x=\frac{\partial\left((1+\kappa x)(ax+bx^2)\right)}{\partial x} \end{equation}

By solve the index of the $0^{th}$ and $1^{th}$ terms, we get $a=1$ and $b=-\kappa/2$. Therefore in a constant field main dipole:

\begin{equation} A_{s,d}\approx-B_0\left(x - \frac{x^2}{2\rho}\right) \end{equation}

For quadrupole, the radius is infinity. The vector potential will be the same as its form in a Cartesian coordinate. From the Beth's representation:

\begin{equation} A_{s,q}=-B_0 b_1(s) \left(x^2-y^2\right)/2=-G(s)\left(x^2-y^2\right)/2 \end{equation}

Here $G(s)$ is the field gradient in unit of [Tesla/m].

We get the Hamiltonian and truncate down to second order, and recall that the rigidity $B_0 \rho=p_0/e$:

\begin{align} K\left(x, p_x, y, p_y,s\right)&=-\left(1+\frac{x}{\rho}\right)\left(1-\frac{p_x^2}{2}-\frac{p_y^2}{2}\right) +\frac{eB_0}{p_0}\left(1+\frac{x}{\rho}\right)\left(x-\frac{x^2}{2\rho}\right)\nonumber\\ +\frac{eB_0b_1}{2p_0}\left(x^2-y^2\right) \nonumber\\ &=-\left(1+\frac{x}{\rho}\right)\left(1-\frac{p_x^2}{2}-\frac{p_y^2}{2}\right) +\frac{1}{\rho}\left(1+\frac{x}{\rho}\right)\left(x-\frac{x^2}{2\rho}\right)+\frac{G(s)}{2B_0\rho}\left(x^2-y^2\right) \nonumber\\ &\approx -1+\frac{p_x^2}{2}+\frac{p_y^2}{2}+\frac{x^2}{2\rho^2}+k(s)\left(x^2-y^2\right)/2 \end{align}

Here we use normalized quadrupole strength $k(s)=G(s)/B_0\rho$. It is the most useful quantity in the simulation to characterize the strength of a quadrupole.

In the approximate Hamiltonian of dipole and normal quadrupole, we can first ignore the constant, then split it in to two separate Hamiltonian, one for each transverse direction. $$K\left(x, p_x, y, p_y,s\right)=K_x\left(x, p_x,s\right)+K_y\left( y, p_y,s\right)$$ This indicate that the motion of two directions are de-coupled and can be treat separately.

\begin{align} K_x\left(x,p_x,s\right)&=\frac{p_x^2}{2}+\frac{x^2}{2\rho^2}+\frac{k(s)x^2}{2} \\ K_y\left(y,p_y,s\right)&=\frac{p_y^2}{2}-\frac{k(s)y^2}{2} \end{align}

The equations of motion can be found from Hamiltonian equation:

\begin{align} & x''+\left(\frac{1}{\rho^2}+k(s)\right)x=0 \\ & y''-k(s)y=0 \end{align}

We returned to the Hill's equation including both the main dipole and quadrupole in the Frenet-Serret coordinate system.

Combine the effect of dipole and quad, there is an chance to have focusing effect in both plane:

\begin{align} -\frac{1}{\rho^2}<k(s)<0 \end{align}

If a machine can have this condition through out, the motion of both transverse plane will be stable. This is called weak focusing.

Transfer Matrix for Linear Elements

It is convenient to write the solution of the Hill's equation using matrix form. Such matrix is named betatron transfer matrix. The motion in each transverse direction can be described by $2\times 2$ matrix.

Drift space

The EOM for drift space is simply $x''=0$ and $y''=0$. Using horizontal direction as instance, if the initial condition at the starting point $s_0$ are $x(s=s_0)=x_0$ and $x'(s=s_0)={x_0'}$. The solution is:

\begin{equation} x(s=s_0+l)=x_0+x'_0l \end{equation}

Or in matrix form

\begin{equation} \left(\begin{array}{c} x\\ x' \end{array}\right)_{s=s_0+l}=\left(\begin{array}{cc} 1 & l\\ 0 & 1 \end{array}\right)\left(\begin{array}{c} x\\ x \end{array}\right)_{s=s_0} \end{equation}

We write the transfer matrix for drift space as:

\begin{equation}M\left(s_0+s,s_0\right)=\left(\begin{array}{cc} 1 & s\\ 0 & 1 \end{array}\right) \end{equation}

Dipole (main)

For the dipole that bend the beam with the radius $\rho$ in horizontal plane, the $x$ direction follows an EOM as a harmonic oscillator, while in $y$ direction, the motion is just as in a drift space.

The matrix form in the dipole can be easily get as:

\begin{equation} \left(\begin{array}{c} x\\ x' \end{array}\right)_{s=s_0+l}=\left(\begin{array}{cc} \cos(l/\rho) & \rho\sin(l/\rho)\\ -\sin(l/\rho)/\rho & \cos(l/\rho) \end{array}\right)\left(\begin{array}{c} x\\ x \end{array}\right)_{s=s_0} \end{equation}

The transfer matrix for dipole in horizontal space is:

\begin{equation} M\left(s_0+l,s_0\right)=\left(\begin{array}{cc} \cos\left(\frac{l}{\rho}\right) & \rho\sin\left(\frac{l}{\rho}\right)\\ -\frac{1}{\rho}\sin\left(\frac{l}{\rho}\right) & \cos\left(\frac{l}{\rho}\right) \end{array}\right) \end{equation}

In a small angle approximation $l/\rho \ll1$, we can expand the matrix up to first order of $l/\rho$. The matrix reduce to a drift space:

\begin{equation} M\left(s_0+l,s_0\right)=\left(\begin{array}{cc} 1 & l\\ 0 & 1 \end{array}\right) \end{equation}

Quadrupole (Normal)

In the example before, we already derived the solution for particle in the quadrupole. It can be rewritten in the matrix form as

\begin{equation} M\left(s_0+s,s_0\right)= \begin{cases} \left(\begin{array}{cc} \cos\left(\sqrt{k}s\right) & \sin\left(\sqrt{k}s\right)/\sqrt{k}\\ -\sqrt{k}\sin\left(\sqrt{k}s\right) & \cos\left(\sqrt{k}s\right) \end{array}\right)&k>0, \text{Focusing Quad}\\ \\ \left(\begin{array}{cc} \cosh\left(\sqrt{k}s\right) & \sinh\left(\sqrt{k}s\right)/\sqrt{k}\\ \sqrt{k}\sinh\left(\sqrt{k}s\right) & \cosh\left(\sqrt{k}s\right) \end{array}\right)&k<0, \text{Defocusing Quad}\\ \end{cases} \end{equation}

Usually, the quadrupole has short length and strong field gradient. We may model such short quad as a 'lens with zero length', and define the focal length $f$ as:

\begin{equation} f=\lim_{l\to 0}\frac{1}{|k|l} \end{equation}

In such limit the transfer map for a thin-length quad is

\begin{equation} M_{\text{quad}}= \begin{cases} \left(\begin{array}{cc} 1 & 0\\ -1/f & 1 \end{array}\right)& \text{Focusing Quad}\\ \\ \left(\begin{array}{cc} 1 & 0\\ 1/f & 1 \end{array}\right)& \text{Defocusing Quad}\\ \end{cases} \end{equation}

Chain of linear elements

In accelerators, we have chain of magnets, as illustrated above. Each magnet has a transfer matrix, the transfer matrix for the chain can be found by matrix multiplications:

\begin{equation} M\left(s_n,s_0\right)=M\left(s_n,s_{n-1}\right)\cdots M\left(s_2,s_1\right)M\left(s_1,s_0\right) \end{equation}

Note that the order of the multiplication is the reverse order of the placement of the magnet.

Properties of the betatron transfer matrix

Symplectic condition

Let's return to the most general form of the matrix $M$:

\begin{equation} M(s)=\left(\begin{array}{cc} a(s) & b(s)\\ c(s) & d(s) \end{array}\right) \end{equation}

The matrix represent the solution:

\begin{align} x(s)&=x_0a(s)+x'_0b(s)\\ x'(s)&=x_0c(s)+x'_0d(s) \end{align}

that satisfy the Hill's equation with general focusing term $K(s)$:

\begin{equation} x''+K(s)x=0 \end{equation}

It can be constructed from the Hamiltonian, which reads:

\begin{equation} H(x,x')=x'^2+K(s)x^2 \end{equation}

The implicit requirement raised from fact that the system is described by a Hamiltonian is called symplecticity. The symplectic property for a 1-D system is the simplest case. To review this property, we use the concept of Wronskian of differential equation. For 1-D second order ODE, there aways exist 2 independent solution. We name it $x_1(s)$ and $x_2(s)$. Wronskian is defined as:

\begin{equation} W\equiv x_1x'_2-x_2x'_1 \end{equation}

And

\begin{align} \frac{dW}{ds}&=x'_1x'_2+x_1x''_2-x'_2x'_1-x_2x''_1 \nonumber\\ &=-K(s)x_1x_2+K(s)x_2x_1 \nonumber\\ &=0 \end{align}

Therefore Wronskian is the constant of motion and not zero(the two solution are indepedent). Then from the map:

\begin{equation} \left(\begin{array}{cc} x_{1} & x_{2}\\ x'_{1} & x'_{2} \end{array}\right)_{s=s_{0}+s}=\left(\begin{array}{cc} a(s) & b(s)\\ c(s) & d(s) \end{array}\right)\left(\begin{array}{cc} x_{1} & x_{2}\\ x'_{1} & x'_{2} \end{array}\right)_{s=s_{0}} \end{equation}

Take the determinant on both sides, we have $$W(s_0+s)=det(M(s_0+s,s_0))W(s)$$

Since $W(s)$ is constant of motion, therefore the determinant of transfer matrix is always 1. This is a sufficient and necessary condition of symplecticity for 1-D system. We will learn later, in multi-dimension problem, this is still a necessary condition but not sufficient anymore.

Long-term stability condition

The symplecticity ensures the area phase space is constant. However, in most times, we have stronger requirements for a repeating chain of elements.

In accelerator, the beam usually will pass through a chain of elements (called cell) many times. A long lattice usually consist of repeating cells. In a ring accelerator, the particle will travels through the lattice billions of turns. If the transfer matrix for one cell is $M$, we are interested in the behavior of repeating $k$ cells $$M_{total}=M\cdot M\cdots M\cdot M=M^k$$ where k is a very large integer or infinity.

To avoid the either $x$ or $x'$ to be very large (or infinity) after k cells. We required that the absolute value of the eigenvalues of $M$ is no larger than 1. The eigenvalue is found by:

\begin{equation} \left|\begin{array}{cc} \lambda-a & -b\\ -c & \lambda-d \end{array}\right|=0 \end{equation}

or

\begin{equation} \lambda^2-(a+d)\lambda+1=0 \end{equation}

To satisfy $|\lambda|\le 1$, we require

\begin{equation} \left|Tr(M)\right|=\left|a+d\right|\le 2 \end{equation}

This is the stable condition for one cell in periodical structure.

Revisit, the 'weak focusing'

Let's write the matrix term for the weak focusing. The weak focusing satisfy:

\begin{align} x''+k_x(s)x&=0 \\ y''+k_y(s)y&=0 \end{align}

Here, $k_x(s)$ and $k_y(s)$ are piece wise constant and always positive, and $k_x(s)+k_y(s)=1/\rho ^2$.

Each constant focusing piece, we have:

\begin{align} M_x\left(s_0+s,s_0\right) &=\left(\begin{array}{cc} \cos\left(\sqrt{k_x}s\right) & \sin\left(\sqrt{k_x}s\right)/\sqrt{k_x}\\ -\sqrt{k_x}\sin\left(\sqrt{k_x}s\right) & \cos\left(\sqrt{k_x}s\right) \end{array}\right)\\ M_y\left(s_0+s,s_0\right) &=\left(\begin{array}{cc} \cos\left(\sqrt{k_y}s\right) & \sin\left(\sqrt{k_y}s\right)/\sqrt{k_y}\\ -\sqrt{k_y}\sin\left(\sqrt{k_y}s\right) & \cos\left(\sqrt{k_y}s\right) \end{array}\right) \end{align}

The matrix indicates that except a scaling factor, the transfer matrix is just a rotation of the phase space. If we define $P_x=x'/\sqrt{k_x}$ and $P_y=y'/\sqrt{k_y}$, the transfer matrix for $(x,P_x)$, $(y,P_y)$ are just rotation matrix with angle $\sqrt{k_{x/y}}s$ respectively. Then we call the angle betatron phase advance.

If we have an entire ring, made of weak focusing structure as shown above, the matrix for entire ring gives:

\begin{align} M_x\left(s_0+C,s_0\right) &=\left(\begin{array}{cc} \cos\left(\sqrt{k_x}2\pi\rho\right) & \sin\left(\sqrt{k_x}2\pi\rho\right)/\sqrt{k_x}\\ -\sqrt{k_x}\sin\left(\sqrt{k_x}2\pi\rho\right) & \cos\left(\sqrt{k_x}2\pi\rho\right) \end{array}\right)\\ M_y\left(s_0+C,s_0\right) &=\left(\begin{array}{cc} \cos\left(\sqrt{k_y}2\pi\rho\right) & \sin\left(\sqrt{k_y}2\pi\rho\right)/\sqrt{k_y}\\ -\sqrt{k_y}\sin\left(\sqrt{k_y}2\pi\rho\right) & \cos\left(\sqrt{k_y}2\pi\rho\right) \end{array}\right) \end{align}

It means, after one turn the phase advance of the two transverse plane is just $2\pi\sqrt{k_{x/y}}\rho$. We define the transverse tune $\nu_{x/y}$ is the one turn phase advance in unit of $2\pi$. Therefore:

\begin{align} \nu_x&=\sqrt{k_{x}}\rho=\sqrt{1+\rho^2 k}\\ \nu_y&=\sqrt{k_{y}}\rho=\sqrt{-k}\rho \end{align}

Therefore for continuous weak-focusing ring, we always have:

\begin{align} \nu_x^2+\nu_y^2=1 \end{align}

We may also look at the weak focusing ring is made of several sections with small drift space $L$ in between. The segment matrix is represented as:

\begin{align} M_{x,seg} &=\left(\begin{array}{cc} 1 & L/2\\ 0 & 1 \end{array}\right) \left(\begin{array}{cc} \cos\left(\sqrt{k_x}s\right) & \sin\left(\sqrt{k_x}s\right)/\sqrt{k_x}\\ -\sqrt{k_x}\sin\left(\sqrt{k_x}s\right) & \cos\left(\sqrt{k_x}s\right) \end{array}\right) \left(\begin{array}{cc} 1 & L/2\\ 0 & 1 \end{array}\right)\\ \end{align}
In [1]:
import sympy
sympy.init_printing(use_unicode=True)
In [2]:
sqkx,s,l=sympy.symbols("\sqrt{k_x},s,L")
def mat_f(sqk,s):
    return sympy.Matrix([[sympy.cos(sqk*s),sympy.sin(sqk*s)/sqk],
                         [-sympy.sin(sqk*s)*sqk,sympy.cos(sqk*s)]])
def mat_d(d):
    return sympy.Matrix([[1,d],[0,1]])

mat_seg=sympy.simplify(mat_d(l/2)*mat_f(sqkx,s)*mat_d(l/2))
mat_seg
Out[2]:
$$\left[\begin{matrix}- \frac{L \sqrt{k_x}}{2} \sin{\left (\sqrt{k_x} s \right )} + \cos{\left (\sqrt{k_x} s \right )} & - \frac{L^{2} \sqrt{k_x}}{4} \sin{\left (\sqrt{k_x} s \right )} + L \cos{\left (\sqrt{k_x} s \right )} + \frac{1}{\sqrt{k_x}} \sin{\left (\sqrt{k_x} s \right )}\\- \sqrt{k_x} \sin{\left (\sqrt{k_x} s \right )} & - \frac{L \sqrt{k_x}}{2} \sin{\left (\sqrt{k_x} s \right )} + \cos{\left (\sqrt{k_x} s \right )}\end{matrix}\right]$$

When $L\sqrt{k_x}\ll1$, the diagnal terms can be written as:

\begin{align} M_{11}=M_{22} &=\cos\left(\sqrt{k_x}s + \Delta\phi\right)\nonumber\\ &\approx\cos\left(\sqrt{k_x}\rho\theta +L\sqrt{k_x}/2\right) \end{align}

where $\Delta \phi=\arcsin(L\sqrt{k_x}/2)$. Here we ignore any terms beyond $L^2 k(x)$. And we know that $L^2 k(x) \le L^2/\rho^2$

Therefore for one turn, which has $N$ such segments, the phase advance is :

\begin{align} \phi(L)&=N\left(\sqrt{k_x}\rho \frac{2\pi}{N} +L\sqrt{k_x}/2\right)\nonumber\\ &\approx 2\pi\rho\sqrt{k_x}\left(1+\frac{NL}{4\pi\rho}\right)\nonumber\\ &=\phi(L=0)\left(1+\frac{NL}{4\pi\rho}\right) \end{align}

And the scaling factor for change the $(x,x')$ to $(x,P_x)$ also changes by an extra factor $\left(1+\frac{NL}{4\pi\rho}\right)$.

The weak-focusing scheme has large limitation since the bending gradient has to be continuous. Later on the theory of strong focusing enables smaller beam size.

In [ ]: