\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,ordinals}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\ddh}[2]{\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{\vecthree}[3]{\left[\begin{array}{c}#1\\#2\\#3\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{\E}{\mathbf E}
\newcommand{\e}{\mathrm e}
\newcommand{\ii}{\mathrm i}
\newcommand{\F}{\mathbf F}
\newcommand{\x}{\mathbf x}
\newcommand{\p}{\mathbf p}
\newcommand{\vv}{\mathbf v}
\newcommand{\uv}{\mathbf u}

\begin{document}

\title{Arbitrary Beams in Uniform Pipes}
\author{Stephen Brooks}
\maketitle

%\tableofcontents\clearpage

\newcommand{\rhot}{\rho_\mathrm{2D}}
\newcommand{\rhol}{\rho_\mathrm{1D}}
\newcommand{\Vt}{V_\mathrm{2D}}
\newcommand{\rb}{r_\mathrm{beam}}
\newcommand{\lap}{\nabla^2_{x,y}}

\section{Beam Slice Problem}
Before considering the whole beam, a simpler problem is a slice (at rest) with charge density \[\rho=\rhot(x,y)\delta(z).\]  The electric field is given by a potential $\E=-\nabla V$ where $V=0$ on a perfectly conducting pipe whose shape does not change in $z$.  Since $\nabla\cdot\E=\rho/\epsilon_0$, the equation to solve is $\nabla^2V=-\rho/\epsilon_0$.

For $z\ne0$, $\rho=0$ so $\nabla^2V=0$ and the problem can be restated as a `time evolution' problem
\[ \dds{V}{z}=-\dds{V}{x}-\dds{V}{y} \qquad \Rightarrow \qquad V''=-\lap V, \]
where prime denotes differentiation by $z$ and $\lap$ is the 2D Laplacian.

\subsection{Eigenmode Decomposition}
Suppose there is a set of 2D eigenpotentials $V_\lambda(x,y)$ for various $\lambda$ such that
\[ \lap V_\lambda = \lambda V_\lambda \qquad \mbox{and} \qquad V_\lambda=0 \quad \mbox{on pipe}. \]
Assuming the eigenpotentials span the space of all 2D potentials within the pipe, decompose the 3D potential at each $z$ via $V=\sum_\lambda a_\lambda(z)V_\lambda$.  Now,
\[ V'' = \sum_\lambda a''_\lambda V_\lambda \qquad \mbox{and} \qquad V'' = -\lap V = -\sum_\lambda a_\lambda \lambda V_\lambda, \]
so equating coefficients gives
\[ a''_\lambda =  -\lambda a_\lambda. \]
Noting that usually for these sorts of eigenmodes, $-\lambda>0$ and that realistically the potential should be exponentially decaying rather than increasing with $z$, the solution is 
\[ a_\lambda(z) = A_\lambda \e^{-\sqrt{-\lambda}z}, \]
for $z\ge0$ and some constants $A_\lambda$.

\subsection{Initial Conditions}
When passing through $z=0$ we expect an infinite value of $V''$ to occur because of the impulse nature of $\rho$ and the $x$ and $y$ contributions to the Laplacian become negligible.  By symmetry, $V(x,y,-z)=V(x,y,z)$ so $V'(x,y,0^-)=-V'(x,y,0^+)$ and the impulse contribution to $V''=\nabla^2 V$ at $z=0$ is:
\[ 2V'(x,y,0^+)\delta(z) = -\rho/\epsilon_0 = \frac{-\rhot(x,y)\delta(z)}{\epsilon_0}
\qquad \Rightarrow \qquad 
V'(x,y,0^+) = \frac{-\rhot(x,y)}{2\epsilon_0}. \]
Note that $V'=\sum_\lambda a'_\lambda V_\lambda$ and $a'_\lambda=-\sqrt{-\lambda} A_\lambda  \e^{-\sqrt{-\lambda}z}$, so
\[ V'(z=0^+) = \sum_\lambda -\sqrt{-\lambda} A_\lambda V_\lambda = \frac{-\rhot}{2\epsilon_0}. \]
Assuming the eigenpotentials are orthogonal and taking the dot product with a particular $V_\lambda$ gives
\[ \sqrt{-\lambda} A_\lambda V_\lambda\cdot V_\lambda = \frac{\rhot\cdot V_\lambda}{2\epsilon_0}
\qquad \Rightarrow \qquad 
A_\lambda = \frac{\rhot\cdot V_\lambda}{2\epsilon_0\sqrt{-\lambda}|V_\lambda|^2}. \]
Alternatively, if $\rhot=\sum_\lambda \rho_\lambda V_\lambda$ then $A_\lambda = \rho_\lambda/(2\epsilon_0\sqrt{-\lambda})$.

\subsection{Beam Slice Potential}
After obtaining $V_\lambda$ with an eigenmode solver of the pipe and the constants $A_\lambda$ using $\rhot$ and the formula in the previous section, the potential anywhere in the pipe is given by
\[ V(x,y,z) = \sum_\lambda A_\lambda \e^{-\sqrt{-\lambda}|z|}V_\lambda(x,y). \]

\section{Separable Charge Density Case}
\newcommand{\id}{\,\mathrm d}
The result for beam slices generalises to separable charge densities (at rest) \[\rho=\rhot(x,y)f(z), \]
since the problem is linear, translatable in $z$ and $f(z)=\int_{-\infty}^{\infty} f(Z)\delta(z-Z)\id Z$,
\[ V(x,y,z) = \sum_\lambda A_\lambda V_\lambda(x,y) \int_{-\infty}^{\infty} f(Z) \e^{-\sqrt{-\lambda}|z-Z|}\id Z. \]

\subsection{Constant Current Beam}
Suppose that $f(z)$ is constant; without loss of generality, $f(z)=1$ with appropriate rescaling of $\rhot$.  Now
\[ \int_{-\infty}^{\infty} f(Z) \e^{-\sqrt{-\lambda}|z-Z|}\id Z = \int_{-\infty}^{\infty} \e^{-\sqrt{-\lambda}|z|}\id z = 2 \int_0^{\infty} \e^{-\sqrt{-\lambda}z}\id z = \frac{2}{\sqrt{-\lambda}}, \]
so the potential no longer depends on $z$:
\[ V(x,y,z) = \sum_\lambda \frac{2A_\lambda }{\sqrt{-\lambda}}V_\lambda(x,y) = \Vt(x,y) = 
\frac1{\epsilon_0}\sum_\lambda \frac{\rhot\cdot V_\lambda}{-\lambda|V_\lambda|^2}V_\lambda(x,y). \]
If $\rhot=\sum_\lambda \rho_\lambda V_\lambda$ then by orthogonality, $\rhot\cdot V_\lambda = \rho_\lambda|V_\lambda|^2$ and
\[ \Vt = \sum_\lambda \frac{\rho_\lambda}{\epsilon_0(-\lambda)}V_\lambda. \]
This can also be derived quickly from the 2D equation $\lap V = -\rho/\epsilon_0$ that assumes $\dd{}z=0$.

\subsection{Longitudinal Space Charge Force}
For longitudinal space charge the field component $E_z=-V'$ is of interest and
\begin{eqnarray*}
V'(x,y,z) &=& \sum_\lambda A_\lambda V_\lambda(x,y) \sqrt{-\lambda} \left( -\int_{-\infty}^{z} f(Z) \e^{\sqrt{-\lambda}(Z-z)}\id Z + \int_{z}^{\infty} f(Z) \e^{\sqrt{-\lambda}(z-Z)}\id Z \right) \\
&=& \sum_\lambda A_\lambda V_\lambda(x,y) \sqrt{-\lambda} \int_{-\infty}^{\infty} f(Z)\,\mathrm{sgn}(Z-z)\,\e^{-\sqrt{-\lambda}|Z-z|}\id Z.
\end{eqnarray*}

\subsection{Derivative of Line Density Approximation}
\newcommand{\iint}{\int\!\!\int\!}
Since $V'=0$ for $f(z)=1$, by linearity constant contributions to $f$ can be neglected when calculating $E_z$.  Also, $V'$ is only sensitive to the value of $f$ `locally', that is in the region a few times $1/\sqrt{-\lambda}$ (the exponential decay length) from $z$.  This means the next Taylor series term $f(Z) = f'(z)(Z-z)$ can be used as an approximation.  Note that defining the line density $\rhol(z)=\iint\rho(x,y,z)\id x \id y=f(z)\iint\rhot$, we have $f'(z) = \rhol'(z)/\iint\rhot$.  The integral in the last section becomes:
\begin{eqnarray*}
\int_{-\infty}^{\infty} f(Z)\,\mathrm{sgn}(Z-z)\,\e^{-\sqrt{-\lambda}|Z-z|}\id Z
&=& f'(z) \int_{-\infty}^{\infty} (Z-z)\mathrm{sgn}(Z-z)\,\e^{-\sqrt{-\lambda}|Z-z|}\id Z \\
&=& f'(z) \int_{-\infty}^{\infty} |Z-z|\e^{-\sqrt{-\lambda}|Z-z|}\id Z \\
&=& 2f'(z) \int_{0}^{\infty} Z\e^{-\sqrt{-\lambda}Z}\id Z \\
&=& 2f'(z) \left[ \left(-\frac{Z}{\sqrt{-\lambda}}-\frac1{-\lambda}\right)\e^{-\sqrt{-\lambda}Z}\right]_{Z=0}^{\infty} \\
&=& 2f'(z) \frac1{-\lambda}.
\end{eqnarray*}
So the potential is approximately
\[ V'(x,y,z) = 2f'(z) \sum_\lambda \frac{A_\lambda}{\sqrt{-\lambda}} V_\lambda(x,y) \quad=\quad f'(z) \Vt(x,y), \]
that is $f'(z)$ times the formula for $V$ in the constant current case.  Substituting the formula for $A_\lambda$ gives
\[ V'(x,y,z) =  f'(z) \sum_\lambda \frac{\rhot\cdot V_\lambda}{\epsilon_0(-\lambda)|V_\lambda|^2} V_\lambda(x,y) \]
and
\[ \frac{E_z}{\rhol'} = \frac{-V'}{f'(z)\iint\rhot} = \frac1{\epsilon_0\iint\rhot} \sum_\lambda \frac{\rhot\cdot V_\lambda}{\lambda|V_\lambda|^2} V_\lambda(x,y)
\quad=\quad \frac{-\Vt(x,y)}{\iint\rhot}. \]
This still depends on $x$ and $y$.  To further approximate, an averaged value over the whole beam can be obtained by dotting this formula with $\rhot/\iint\rhot$:
\[ \left<\frac{E_z}{\rhol'}\right> = \frac1{\epsilon_0(\iint\rhot)^2} \sum_\lambda \frac{(\rhot\cdot V_\lambda)^2}{\lambda|V_\lambda|^2}
\quad=\quad \frac{-\rhot\cdot \Vt}{(\iint\rhot)^2}. \]

\subsection{Relationship to the $g$ Factor}
A conventional way of writing the longitudinal space charge field (e.g. \cite{WiedemannEz}) is
\[ E_z=E_{z,\mathrm{wall}}-\frac{q}{4\pi\epsilon_0}\frac1{\gamma^2}\left(1+2\ln\frac{r_\mathrm{wall}}{\rb}\right)\pdd{\lambda}{z}, \]
where $q$ is the charge on an individual particle, $\gamma$ is the gamma factor of a moving beam and $\lambda$ is the particle line density.  So far this report has considered perfectly-conducting walls where $E_{z,\mathrm{wall}}=0$ and stationary `beams' for which $\gamma=1$.  The charge line density $\rhol=q\lambda$.  In this case,
\[ E_z=-\frac{1}{4\pi\epsilon_0}\left(1+2\ln\frac{r_\mathrm{wall}}{\rb}\right)\rhol', \]
where the dimensionless term in brackets is called the `$g$ factor'.  This formula assumes a circular pipe and the $g$ factor differs for other shapes.  In general, this may be related to the quantity in the previous section via
\[ E_z=-\frac{1}{4\pi\epsilon_0}g\rhol' \qquad \Rightarrow \qquad g=-4\pi\epsilon_0\left<\frac{E_z}{\rhol'}\right>
\quad=\quad \frac{4\pi}{(\iint\rhot)^2} \sum_\lambda \frac{(\rhot\cdot V_\lambda)^2}{-\lambda|V_\lambda|^2}. \]
In terms of the 2D potential,
\[ g = 4\pi\epsilon_0 \frac{\rhot\cdot \Vt}{(\iint\rhot)^2}
\qquad\mbox{and}\qquad \left<E_z\right>=-\rhol'\frac{\rhot\cdot \Vt}{(\iint\rhot)^2}. \]

\subsection{Example: Circular Pipe}
For a circularly symmetric beam in a circular pipe, $\rhot\cdot V_\lambda=0$ for eigenpotentials that vary with $\theta$, leaving only those expressable as $V(r)$.  The eigenpotentials must satisfy
\[ \lap V(r) = \frac1r\pdd{V}{r}+\pdds{V}{r}=\lambda V \qquad \Rightarrow \qquad r^2\pdds{V}{r}+r\pdd{V}{r}-\lambda r^2 V=0, \]
where $V'(0)=0$ by differentiability at the centre and $V(R)=0$ boundary condition for a pipe of radius $R$.  Changing to a rescaled variable $x$ with $x^2=-\lambda r^2$ turns this into Bessel's equation
\[ x^2\pdds{V}{x}+x\pdd{V}{x}+(x^2-\alpha^2) V=0 \]
for $\alpha=0$, with boundary conditions $V'(0)=0$ and $V(x=\sqrt{-\lambda}R)=0$.  This has the solution $V=J_0(x)=J_0(\sqrt{-\lambda}r)$ but to satisfy the boundary condition, $\sqrt{-\lambda}R=J_{0n}$ for some $n$, where $J_{0n}$ is the $n$\th root of $J_0$ (for which there is no analytic expression).  Thus the eigenvalues are $\lambda_n=-(J_{0n}/R)^2$ and $V_{\lambda_n}$ will be abbreviated to $V_n(r) = J_0((r/R)J_{0n})$.

The established result $1+2\ln(R/\rb)$ for a circular uniform beam is actually calculated for the field at the centre of the pipe and not averaged over the beam, so define the position-dependent $g$ factor as
\[ g(x,y) = -4\pi\epsilon_0\frac{E_z(x,y)}{\rhol'} \quad=\quad \frac{4\pi}{\iint\rhot} \sum_\lambda \frac{\rhot\cdot V_\lambda}{-\lambda|V_\lambda|^2} V_\lambda(x,y)
\quad=\quad 4\pi\epsilon_0\frac{\Vt(x,y)}{\iint\rhot}. \]

\subsubsection{Circular Uniform Beam}
Assume a uniform, normalised beam for which $\rhot=1/(\pi \rb^2)$ for $r\le \rb$ and zero outside.  Noting $\iint\rhot=1$ and $J_0(0)=1$, the value of interest is
\[ g(0,0) = 4\pi\sum_n \frac{\rhot\cdot V_n}{-\lambda_n|V_n|^2} V_n(r=0) 
\quad=\quad 4\pi\sum_n \frac{\rhot\cdot V_n}{(J_{0n}/R)^2|V_n|^2}. \]
Consider the dot product
\[ \rhot\cdot V_n = \frac1{\pi \rb^2}\int_{r=0}^{\rb} J_0((r/R)J_{0n})\,2\pi r \id r
\quad=\quad \frac2{\rb^2}\int_{r=0}^{\rb} J_0((r/R)J_{0n})\,r \id r, \]
rescale the integral with $x=(r/R)J_{0n}$, giving $r\id r = (R/J_{0n})^2\,x \id x$:
\[ \rhot\cdot V_n = \frac2{((\rb/R)J_{0n})^2}\int_{x=0}^{\frac{\rb}{R}J_{0n}} xJ_0(x) \id x. \]
Properties of the Bessel functions give that $xJ_0(x)$ is the derivative of $xJ_1(x)$, so
\[ \rhot\cdot V_n = \frac2{(\rb/R)J_{0n}} J_1((\rb/R)J_{0n})
\quad=\quad \frac2{fJ_{0n}} J_1(fJ_{0n}), \]
where the beam radius fraction $f=\rb/R$ has been defined (so the textbook result $1+2\ln(R/\rb)$ becomes $1-2\ln f$).
The other integral to calculate is
\[ |V_n|^2 = \int_{r=0}^{R} J_0((r/R)J_{0n})^2\,2\pi r \id r
\quad=\quad 2\pi R^2 \int_{x=0}^{1} xJ_0(xJ_{0n})^2\id x, \]
where the substitution $x=r/R$ and $r\id r = R^2\,x\id x$ has been used.  Compare with the orthogonality relation for Bessel functions
\[ \int_0^1 xJ_\alpha(xJ_{\alpha m})J_\alpha(xJ_{\alpha n})\id x = \frac{\delta_{mn}}{2}J_{\alpha+1}(J_{\alpha m})^2, \]
which for $\alpha=0$ and $m=n$ becomes
\[ \int_0^1 xJ_0(xJ_{0n})^2\id x = \frac{1}{2}J_1(J_{0n})^2, \]
therefore
\[ |V_n|^2 = \pi R^2 J_1(J_{0n})^2. \]
Substituting these into the formula gives
\[ g(0,0) = 4\pi\sum_n \frac{\frac2{fJ_{0n}} J_1(fJ_{0n})}{(J_{0n}/R)^2 \pi R^2 J_1(J_{0n})^2}
\quad=\quad 8\,\sum_n \frac{J_1(fJ_{0n})}{fJ_{0n}^3 J_1(J_{0n})^2}. \]
The beam-averaged $g$ factor in this case just has an additional factor of $\rhot\cdot V_n$:
\[ g = 16\,\sum_n \frac{J_1(fJ_{0n})^2}{f^2J_{0n}^4 J_1(J_{0n})^2}. \]
Numerical evaluation (30000 terms) shows that to within $10^{-7}$, this $g(0,0)=1-2\ln f$ as required and $g=\frac12-2\ln f$.

\subsubsection{Circular Waterbag Beam}
From results in the last section, for any beam normalised by $\iint\rhot=1$ in a circular pipe,
\[ g = 4\,\sum_n \frac{(\rhot\cdot V_n)^2}{J_{0n}^2 J_1(J_{0n})^2} \]
is the averaged $g$ factor, where $V_n = J_0((r/R)J_{0n})$.  For a waterbag beam, $\rhot\propto 1-(r/\rb)^2$ for $r\le \rb$ and zero outside, so when normalised, $\rhot = (2/(\pi \rb^2))(1-(r/\rb)^2)$.  The term that needs calculating is
\begin{eqnarray*}
\rhot\cdot V_n &=& \frac2{\pi \rb^2}\int_{r=0}^{\rb} (1-(r/\rb)^2)\,J_0((r/R)J_{0n})\,2\pi r \id r \\
&=& \frac4{\rb^2}\int_{r=0}^{\rb} (1-(r/\rb)^2)\,J_0((r/R)J_{0n})\, r \id r.
\end{eqnarray*}
Now perform the same rescaling as before with $x=(r/R)J_{0n}$, $r\id r = (R/J_{0n})^2\,x \id x$ and $r/\rb=x/(fJ_{0n})$:
\begin{eqnarray*}
\rhot\cdot V_n &=& \frac{4R^2}{\rb^2J_{0n}^2}\int_{x=0}^{fJ_{0n}} (1-(x/(fJ_{0n}))^2)\,J_0(x)\, x \id x \\
&=& \frac{4}{f^2J_{0n}^2}\int_{x=0}^{fJ_{0n}} xJ_0(x) -\frac1{f^2J_{0n}^2}\,x^3J_0(x) \id x.
\end{eqnarray*}
Since $J_0'=-J_1$ and $(xJ_1)'=xJ_0$, it can be seen that $(x^3-4x)J_1+2x^2J_0$ differentiates to $x^3J_0$.  Therefore
\begin{eqnarray*}
\rhot\cdot V_n &=& \frac{4}{f^2J_{0n}^2}\left[xJ_1(x)-\frac{(x^3-4x)J_1(x)+2x^2J_0(x)}{f^2J_{0n}^2}\right]_{x=0}^{fJ_{0n}} \\
&=& \frac{4}{f^2J_{0n}^2}\left(fJ_{0n}J_1(fJ_{0n})-\frac{((fJ_{0n})^3-4fJ_{0n})J_1(fJ_{0n})+2(fJ_{0n})^2J_0(fJ_{0n})}{f^2J_{0n}^2}\right) \\
&=& \frac{4}{f^2J_{0n}^2}\left(fJ_{0n}J_1(fJ_{0n})-fJ_{0n}J_1(fJ_{0n})-2J_0(fJ_{0n})+\frac{4J_1(fJ_{0n})}{fJ_{0n}}\right) \\
&=& \frac{8}{f^2J_{0n}^2}\left(\frac{2J_1(fJ_{0n})}{fJ_{0n}}-J_0(fJ_{0n})\right) \qquad \left(\mbox{c.f.}\quad\frac2{fJ_{0n}} J_1(fJ_{0n})\quad\mbox{for uniform beam}\right).
\end{eqnarray*}
Substituting back into the formula for $g$ gives
\[ g = 256\,\sum_n \frac{\left(\frac{2J_1(fJ_{0n})}{fJ_{0n}}-J_0(fJ_{0n})\right)^2}{f^4J_{0n}^6 J_1(J_{0n})^2}. \]
Numerical computation (30000 terms) shows that to within $10^{-11}$, $g=\frac{11}{12}-2\ln f$.

\subsubsection{Potential Method for General Circularly-Symmetric Beam}
If the constant current potential $\Vt$ is known, the $g$ factor $4\pi\epsilon_0(\rhot\cdot\Vt)/(\iint\rhot)^2$ can be found by calculating $\rhot\cdot\Vt$ directly.  In the case of a circular pipe with perfectly conducting boundary at $r=R$, $\Vt=V(r)$ is determined by
\[ \lap V(r)=\frac1r V'+V''=-\frac{\rho}{\epsilon_0}, \qquad V(R)=0, \]
where prime denotes differentiation by $r$ not $z$ here.  The integral to be evaluated is
\[ \rhot\cdot\Vt = \int_0^R \rho(r) V(r) 2\pi r \id r \]
and while there is no general simplification of this, Baartman provides solutions for a wide range of cases using this method in \cite{Baartman}.

\subsection{Example: Rectangular Pipe}
A rectangular beam pipe $[0,X]\times[0,Y]$ has eigenpotentials 
\[ V_{nm}=\sin(n\pi x/X)\sin(m\pi y/Y) \qquad \mbox{with} \qquad \lambda_{nm}=-(n\pi/X)^2-(m\pi/Y)^2 \]
for integers $n,m\ge1$.  Calculations will also require the value
\[ |V_{nm}|^2 = \int_0^X\int_0^Y V_{nm}^2 \id x \id y = \int_0^X \sin^2(n\pi x/X) \id x \int_0^Y \sin^2(m\pi y/Y) \id y = \frac{XY}{4}. \]
Substituting these into the $g$ factor formula gives
\[ g = \frac{16}{\pi XY(\iint\rhot)^2} \sum_{n,m\ge1} \frac{(\rhot\cdot V_{nm})^2}{(n/X)^2+(m/Y)^2}. \]

\subsubsection{Uniform Rectangular Beam}
Suppose the beam occupies the rectangle $[x_0,x_1]\times[y_0,y_1]$ with $\rhot=1$, so that $\iint\rhot=(x_1-x_0)(y_1-y_0)=\Delta x\Delta y$.  Calculate
\begin{eqnarray*}
\rhot\cdot V_{nm} &=& \int_{x_0}^{x_1}\int_{y_0}^{y_1} \sin(n\pi x/X)\sin(m\pi y/Y) \id x \id y \\
&=& \left[\frac{-X}{n\pi}\cos(n\pi x/X)\right]_{x_0}^{x_1} \left[\frac{-Y}{m\pi}\cos(m\pi y/Y)\right]_{y_0}^{y_1} \\
&=& \frac{XY}{nm\pi^2}(\cos(n\pi x_1/X)-\cos(n\pi x_0/X))(\cos(m\pi y_1/Y)-\cos(m\pi y_0/Y)).
\end{eqnarray*}
So that 
\[ g = \frac{16XY}{\pi^5\Delta x^2\Delta y^2} \sum_{n,m\ge1} \frac{(\cos(n\pi x_1/X)-\cos(n\pi x_0/X))^2(\cos(m\pi y_1/Y)-\cos(m\pi y_0/Y))^2}{n^2m^2((n/X)^2+(m/Y)^2)}. \]

\subsubsection{Elliptical Gaussian Beam}
\newcommand{\M}{\mathbf M}\newcommand{\Q}{\mathbf Q}
\newcommand{\bv}{\mathbf b}\newcommand{\cv}{\mathbf c}
An elliptical Gaussian beam has $\rhot=\e^{-(ax^2+2bxy+cy^2)}=\e^{-\x^T\M\x}$ where $\M=\small\mxtwo{a}{b}{b}{c}$ and $\x=\small\vectwo{x}{y}$.  It is well known that $\iint\e^{-\x^T\x} = \iint\e^{-(x^2+y^2)} = \pi$, so composing an arbitrary linear function $\x\rightarrow\Q\x$ gives $\iint\e^{-\x^T\Q^T\Q\x}=\pi/\det\Q$, thus if $\Q^T\Q=\M$ then $(\det\Q)^2=\det\M$ and $\iint\e^{-\x^T\M\x}=\pi/\sqrt{\det\M}$, that is $\iint\rhot = \pi/\sqrt{ac-b^2}.$

For a Gaussian centred on $(\bar x,\bar y)$ that is not significantly outside the beam pipe, make the approximation to integrate over the whole $x$--$y$ plane
\begin{eqnarray*}
\rhot\cdot V_{nm} &=& \int_{-\bar x}^{X-\bar x}\int_{-\bar y}^{Y-\bar y} \e^{-(ax^2+2bxy+cy^2)} \sin(n\pi (\bar x+x)/X)\sin(m\pi (\bar y+y)/Y) \id x \id y \\
&\simeq& \iint \e^{-(ax^2+2bxy+cy^2)} \sin(n\pi (\bar x+x)/X)\sin(m\pi (\bar y+y)/Y).
\end{eqnarray*}
Note that $\sin x = \frac{1}{2\ii}(\e^{\ii x}-\e^{-\ii x})$ so that
\begin{eqnarray*}
\rhot\cdot V_{nm} &\simeq& \frac1{-4} \sum_{j,k=-1,1}\iint \e^{-(ax^2+2bxy+cy^2)} j \e^{j \ii n\pi (\bar x+x)/X} k \e^{k \ii m\pi (\bar y+y)/Y} \\
&=& -\frac14 \sum_{j,k=-1,1} jk \e^{\ii\pi(jn\bar x/X + km\bar y/Y)} \iint \e^{-(ax^2+2bxy+cy^2-(\ii\pi jn/X)x-(\ii\pi km/Y)y)}.
\end{eqnarray*}
The negative exponent is of the form $\x^T\M\x+\bv^T\x$, where $\bv^T=-\ii\pi[jn/X,km/Y]$.  Noting that $\M=\M^T$, a shift of origin expands as $(\x+\cv)^T\M(\x+\cv) = \x^T\M\x + 2\cv^T\M\x + \cv^T\M\cv$.  This is equal to $\x^T\M\x + \bv^T\x + \cv^T\M\cv$ if one puts $\cv=\frac12\M^{-1}\bv$.  Thus $\x^T\M\x+\bv^T\x = (\x+\cv)^T\M(\x+\cv) - \frac14\bv^T\M^{-1}\bv$.  The shift won't affect the integral while the constant term in the exponent becomes a constant factor:
\begin{eqnarray*}
\rhot\cdot V_{nm} &\simeq& -\frac14 \sum_{j,k=-1,1} jk \e^{\ii\pi(jn\bar x/X + km\bar y/Y)}\e^{\frac14\bv^T\M^{-1}\bv} \iint \e^{-(ax^2+2bxy+cy^2)} \\
&=& -\frac\pi{4\sqrt{ac-b^2}} \sum_{j,k=-1,1} jk \e^{\ii\pi(jn\bar x/X + km\bar y/Y)}\e^{\frac14\bv^T\M^{-1}\bv}.
\end{eqnarray*}
To evaluate the exponent, note that $\M^{-1}=\frac1{ac-b^2}\small\mxtwo{c}{-b}{-b}{a}$ so
\begin{eqnarray*}
\frac14\bv^T\M^{-1}\bv &=& -\frac{\pi^2}{4}[jn/X,km/Y]\frac1{ac-b^2}\mxtwo{c}{-b}{-b}{a}\vectwo{jn/X}{km/Y} \\
&=& -\frac{\pi^2}{4(ac-b^2)}\left(c(jn/X)^2-2b(jn/X)(km/Y)+a(km/Y)^2\right) \\
&=& -\frac{\pi^2}{4(ac-b^2)}\left(c(n/X)^2-2bjk(n/X)(m/Y)+a(m/Y)^2\right),
\end{eqnarray*}
where the last line has used $j^2=k^2=1$.  The summand apart from the $\e^{\ii\pi C}$ term now only depends on $jk$, which shall be renamed to $j$ and collect pairs of terms which differ by both $j$ and $k$ changing sign, which inverts the sign of $C$, making $\e^{\ii\pi C}+\e^{-\ii\pi C}=2 \cos \pi C$.
\begin{eqnarray*}
\rhot\cdot V_{nm} &\simeq& -\frac\pi{2\sqrt{ac-b^2}} \sum_{j=-1,1} j \cos \pi C \exp E \\
&=& -\frac12 \iint\rhot \sum_{j=-1,1} j \cos \pi C \exp E,
\end{eqnarray*}
where $C=n\bar x/X + jm\bar y/Y$ and $E = -\frac{\pi^2}{4(ac-b^2)}\left(c(n/X)^2-2bj(n/X)(m/Y)+a(m/Y)^2\right)$.  Aside from the cancellation of $\iint\rhot$ no further simplification is possible on substituting into the $g$ factor formula.

%\subsubsection{General Separable Beam}

%\subsubsection{Piecewise Polynomial Separable Beam}

\section{Moving Beams (Longitudinal Space Charge)}
Previous sections have only considered a stationary beam, or alternatively a beam in its co-moving frame (assuming negligible relative velocities between particles).  Writing the co-moving quantities with a tilde, the $g$ factor relationship is $\tilde E_z=-\frac{1}{4\pi\epsilon_0}g\tilde\rhol'$.  The velocity of the beam is parallel to the $z$ axis, so the parallel component of the electric field $E_z=\tilde E_z$ is unchanged under frame transformation and produces no other field components.

The line density is less trivial to transform.  Define $Q(z)=\int_{-\infty}^z\rhol(Z)\id Z$ to be the amount of charge at positions $Z\le z$.  $Q$ transforms as a scalar quantity whereas the density $\rhol=\dd{Q}{z}$ is actually a rank 1 covariant tensor and $\rhol'=\dds{Q}{z}$ is rank 2 covariant.  As the beam is larger in its rest frame by $\tilde z=\gamma z$, this means $\rhol'=(\dd{\tilde z}{z})^2\dds{Q}{\tilde z}=\gamma^2\tilde\rhol'$.

Putting this all together gives the $g$ factor relation for a moving beam as 
\[ E_z=-\frac{1}{4\pi\epsilon_0}\frac{g}{\gamma^2}\rhol'. \]

\subsection{Space Charge Energy Gain (or loss) per Turn}
If a turn in the machine has length $L$, then $\Delta E = L \dd{E}{z} = LF_z = LqE_z$, therefore
\[ \Delta E = -\frac{Lq}{4\pi\epsilon_0}\frac{g}{\gamma^2}\rhol'. \]
This is equivalent to a voltage of
\[ V = \frac{L}{4\pi\epsilon_0}\frac{g}{\gamma^2}\rhol'. \]

%\section{General Charge Density Case}% Do beam slices before this?

\begin{thebibliography}{9}   % Use for  1-9  references
%\begin{thebibliography}{99} % Use for 10-99 references
\newcommand{\link}[1]{\texttt{#1}}

\bibitem{WiedemannEz} \emph{Particle Accelerator Physics II}, 2\nd Edition, H. Wiedemann, section 10.2.4 \emph{Longitudinal Space-Charge Field}, equation (10.55).
\bibitem{Baartman} \emph{Form Factor g In Longitudinal Space Charge Impedance}, R. Baartman, TRIUMF design note 1992-TRI-DN-K206, available from \link{http://lin12.triumf.ca/text/design\_notes/}-\link{k206/k206h.pdf}.

\end{thebibliography}


\end{document}