\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,amsmath,amssymb,bm,graphicx}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdds}[2]{\frac{\partial^2 #1}{\partial{#2}^2}}
\newcommand{\pddd}[3]{\frac{\partial^2 #1}{\partial #2\partial #3}}

\newcommand{\rowfour}[4]{\left[\begin{array}{cccc}#1&#2&#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\newcommand{\mxtwo}[4]{\left[\begin{array}{cc}#1&#2\\#3&#4\end{array}\right]}

\newcommand{\sq}[1]{\left(#1\right)^2}
\newcommand{\floor}[1]{\left\lfloor#1\right\rfloor}
\newcommand{\abs}[1]{\left|#1\right|}

\newcommand{\rem}[1]{}

\newcommand{\A}{\mathbf A}
\newcommand{\E}{\mathbf E}
\newcommand{\B}{\mathbf B}
\newcommand{\F}{\mathbf F}
\newcommand{\J}{\mathbf J}

\newcommand{\x}{\mathbf x}
\newcommand{\y}{\mathbf y}
\newcommand{\z}{\mathbf z}

\newcommand{\p}{\mathbf p}
\newcommand{\q}{\mathbf q}
\newcommand{\s}{\mathbf s}
\newcommand{\rv}{\mathbf r}

\newcommand{\uv}{\mathbf u}
\newcommand{\vv}{\mathbf v}

\newcommand{\e}{\mathbf e}
\newcommand{\f}{\mathbf f}
\newcommand{\g}{\mathbf g}

\newcommand{\de}{\,\mathrm d}

\newcommand{\st}{$^\mathrm{st}$ }
\newcommand{\nd}{$^\mathrm{nd}$ }
\newcommand{\rd}{$^\mathrm{rd}$ }
\newcommand{\thh}{$^\mathrm{th}$ }

\begin{document}

\title{First Order Transport of a Cold Uniform Ellipsoid of Charge}
\author{Stephen Brooks}
\maketitle

\section{Definitions}
A single particle has a 6-element phase space vector $\s=(\x,\gamma\vv)$ where $\gamma\vv=\frac{\p}{m}$ is proportional to the relativistic momentum.

The ellipsoid is represented by a 21-element vector $(\s,C,V)$ where $\s$ is the phase space location of the centroid.  $C$ is the $3\times3$ covariance matrix defined by $C_{ij}=\langle \delta x_i \delta x_j \rangle$, where $\delta\x$ is the distance of a particle from the centroid $\x$.  $C$ is symmetric so only counts as 6 elements of the overall vector.  $V$ is a general $3\times3$ matrix that determines the velocity distribution in the bunch by $\delta(\gamma\vv)=V\delta\x$.  In total this makes $6+6+9=21$ elements.

\begin{figure}[!htb]
   \centering
   \includegraphics*[width=\columnwidth]{bunchmodel}
   %\caption{21-variable bunch model.}
   %\label{fig:bunchmodel}
\end{figure}

The ellipsoid is uniformly filled with total charge $Q$ and has zero temperature, so occupies only three dimensions of the six dimensional phase space.  It can be seen as the image of the unit ball $\{\uv:|\uv|\le1\}$ under the mapping $\uv\mapsto\s+(X\uv,VX\uv)$, where $X$ is some $3\times3$ matrix that maps the unit ball to an ellipsoid with covariance matrix $C$.  The two are related by 
\[ C = \langle\delta\x\delta\x^T\rangle = \langle X\uv\uv^TX^T \rangle = X\langle\uv\uv^T\rangle X^T = X\frac15IX^T = \frac15XX^T, \]
where the fact $\langle\uv\uv^T\rangle=\frac15I$ is related to the RMS of a coordinate in a unit ball being $\frac1{\sqrt{5}}$.

\section{Time Derivative of $\s$}
The time derivative of position is $\dot\x=\vv=\frac1{\gamma}\gamma\vv$, so a way to calculate $\frac1{\gamma}$ from $\gamma\vv$ is needed.  Relativistic formulae give $|\gamma\vv|^2=(\beta\gamma)^2c^2=(\gamma^2-1)c^2$, so  $c^2+|\gamma\vv|^2=(\gamma c)^2$ and therefore $\frac1{\gamma}=\frac{c}{\sqrt{c^2+|\gamma\vv|^2}}$.

The time derivative of $\gamma\vv$ is determined by the Lorentz force law $\F=q(\E+\vv\times\B)$ in the fields $\E(\x,t)$ and $\B(\x,t)$.  Note that $\dot{(\gamma\vv)}=\frac{\dot\p}{m}=\frac{\F}{m}$, so $\dot{(\gamma\vv)}=\frac{q}{m}(\E+\frac1\gamma\gamma\vv\times\B)$.

\section{Time Derivative of $C$}
Differentiating a previous formula, $\dot C = \frac15(\dot X X^T + X \dot X^T)$, the second term being the transpose of the first.  To evaluate $\dot X$, remember that $\delta\x=X\uv$ and $\delta(\gamma\vv)=VX\uv$ and use the phase space flow
\[ \delta\dot\x = \pdd{\dot\x}{\x}\delta\x + \pdd{\dot\x}{(\gamma\vv)}\delta(\gamma\vv), \]
since the time derivative of the phase space position is defined entirely in terms of the phase space position.  The Jacobian matrix is defined as $\left(\pdd{\f}{\g}\right)_{ij}=\pdd{f_i}{g_j}$ so that $\pdd{\f}{\g}\delta\g \simeq \delta\f$ for small displacements.  For this case, $\pdd{\dot\x}{\x}=0$ as $\dot\x$ only depends on $\gamma\vv$.  The vector $\uv$ may be considered a constant that `labels' the particle, so $\delta\dot\x=\dot X \uv$ and
\[ \dot X \uv = \pdd{\dot\x}{(\gamma\vv)}VX\uv \]
for any $\uv$, so $\dot X = \pdd{\dot\x}{(\gamma\vv)}VX$.  This gives  $\dot C = \frac15\pdd{\dot\x}{(\gamma\vv)}VX X^T + \mbox{transpose} = \pdd{\dot\x}{(\gamma\vv)}VC + \mbox{transpose}$.

\section{Time Derivative of $V$}
Take the other part of the phase space flow
\[ \delta\dot{(\gamma\vv)} = \pdd{\dot{(\gamma\vv)}}{\x}\delta\x + \pdd{\dot{(\gamma\vv)}}{(\gamma\vv)}\delta(\gamma\vv), \]
substituting in the matrix definitions (and removing $\uv$) gives
\[ \dot V X + V \dot X = \pdd{\dot{(\gamma\vv)}}{\x}X + \pdd{\dot{(\gamma\vv)}}{(\gamma\vv)}VX. \]
Subtracting $V \dot X$ both sides and postmultiplying by $X^{-1}$ gives
\[ \dot V = \pdd{\dot{(\gamma\vv)}}{\x} + \pdd{\dot{(\gamma\vv)}}{(\gamma\vv)}V - V \dot X X^{-1}. \]
Using the formula for $\dot X$ from the previous section,
\[ \dot V = \pdd{\dot{(\gamma\vv)}}{\x} + \pdd{\dot{(\gamma\vv)}}{(\gamma\vv)}V - V \pdd{\dot\x}{(\gamma\vv)}V. \]

\section{Derivatives of Phase Space Flow}
The term $\pdd{\dot\x}{(\gamma\vv)}$ would be equal to $I$ in a non-relativistic setting, so it encodes relativistic kinematics.
\[ \pdd{\dot\x}{(\gamma\vv)} = \frac{-1}{c^2\gamma^3}(\gamma\vv)(\gamma\vv)^T + \frac1{\gamma}I. \]

The term $\pdd{\dot{(\gamma\vv)}}{\x}$ contains contributions from the field gradients:
\[ \pdd{\dot{(\gamma\vv)}}{x_i} = \frac{q}{m}\left(\pdd{\E}{x_i}+\frac1\gamma\gamma\vv\times\pdd{\B}{x_i}\right). \]

The term $\pdd{\dot{(\gamma\vv)}}{(\gamma\vv)}$ comes from the magnetic force's dependence on the velocity $\dot\x=\frac1\gamma\gamma\vv$:
\[ \pdd{\dot{(\gamma\vv)}}{(\gamma v_i)} = \frac{q}{m}\pdd{\dot\x}{(\gamma v_i)}\times\B = \frac{q}{m}\left( \frac{-1}{c^2\gamma^3}(\gamma\vv)\gamma v_i + \frac1{\gamma}\e_i \right)\times\B. \]

\section{Electrostatic Repulsion}
A derivation using elliptic coordinates \cite{CaiEllipsoid,KelloggPotentialTheory} gives the interior potential of an axis-aligned uniformly charged ellipsoid as
\[ U(\x) = \pi abc\rho \int_0^\infty \left( 1-\frac{x^2}{a^2+v}-\frac{y^2}{b^2+v}-\frac{z^2}{c^2+v} \right) \frac1{\sqrt{(a^2+v)(b^2+v)(c^2+v)}} \de v, \]
where $a,b,c$ are the radii in $x,y,z$ and $\rho=Q/V=Q/(\frac43\pi abc)$ is the charge density.  This integral must be done numerically.  The corresponding electric field is
\[ E_i = -\dd{U}{x_i} = \pi abc\rho \int_0^\infty \frac{2x_i}{a_i^2+v} \frac1{\sqrt{(a^2+v)(b^2+v)(c^2+v)}} \de v \]
\[ = \frac32Q\left( \int_0^\infty \frac{1}{(a_i^2+v)\sqrt{(a^2+v)(b^2+v)(c^2+v)}} \de v \right) x_i, \]
where $a_i$ is the radius in the $x_i$ axis.  This field varies linearly with $\x$ and gives a gradient matrix $\pdd{\E}{\x}$ that is constant and diagonal.  For non-axis-aligned ellipsoids, the principal axes must be found (for example by SVD of $C$) and the field rotated.  The rotated electric field will remain linear and the gradient matrix will remain constant but not diagonal.  This self electric field can be added onto any external field and gradients being applied.

\subsection{Moving Reference Frame}
The above derivation assumes that the charges are stationary for all time and thus does not generate a $\B$ field.  For moving bunches this is not valid and a much better approximation is to assume the charges are stationary in the average rest frame of the bunch (this is still not entirely accurate if the bunch changes shape or accelerates).

Define $\bm\beta=\vv/c$, $\beta=|\vv|/c$ and the $3\times3$ spatial part of the Lorentz transformation matrix
\[ \Lambda=I+\frac{\gamma-1}{\beta^2}\bm\beta\bm\beta^T. \]
The covariance matrix in the bunch's rest frame is $C_\mathrm{bunch}=\Lambda C\Lambda$.  This may be used in the method in the previous section to give  $\pdd{\E_\mathrm{bunch}}{\x_\mathrm{bunch}}$.  The transformation of $\E_\mathrm{bunch}$ back to the laboratory frame also generates a $\B$ field:
\[ \E = \gamma\E_\mathrm{bunch}-(\gamma-1)\E_\mathrm{bunch}^\parallel, \qquad \B=\frac1c \gamma \bm\beta\times\E_\mathrm{bunch}, \]
where $\E_\mathrm{bunch}^\parallel = \frac{(\E_\mathrm{bunch}\cdot\bm\beta)\bm\beta}{\beta^2}$ is the electric field parallel to the bunch velocity.  Finally, the spatial derivatives can be transformed back to the laboratory frame with
\[ \pdd{\E}{\x} = \pdd{\E}{\x_\mathrm{bunch}}\Lambda, \qquad \pdd{\B}{\x} = \pdd{\B}{\x_\mathrm{bunch}}\Lambda. \]

\section{Numerical Time Integration}
Now that all the time derivatives of the variables $(\s,C,V)$ have been defined, any numerical integration method such as 4\thh order Runge--Kutta may be used on the 21-element vector.

This model is called ``first order transport'' because it is first order in space.  The change in shape of the ellipsoid is represented by linear transformations (matrices) and no higher-order distortions are taken into account.  


\begin{thebibliography}{9}

\bibitem{CaiEllipsoid} \emph{Potential Field of a Uniformly Charged Ellipsoid}, Wei Cai (Department of Mechanical Engineering, Stanford University), May 28, 2007.

\bibitem{KelloggPotentialTheory} \emph{Foundations of Potential Theory}, O. D. Kellogg, (Dover, New York, 1953).

\end{thebibliography}

\end{document}