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

\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{\B}{\mathbf B}
\newcommand{\F}{\mathbf F}
\newcommand{\x}{\mathbf x}
\newcommand{\p}{\mathbf p}
\newcommand{\vv}{\mathbf v}
\newcommand{\uv}{\mathbf u}

\begin{document}

\title{Beam Distribution with $e^{-x}$-like Tails}
\author{Stephen Brooks}
\maketitle

\section{Hyperspherically Symmetric 4D Distributions}
The distributions in this document will be generated in normalised 4D space $\x=(x_1,x_2,x_3,x_4)$, requiring a linear transformation to be scaled to physical phase space e.g. $(x,x',y,y')$.  Their density will be a function of radius $r=|\x|$, which is proportional to single particle `amplitude' or Courant--Snyder invariant in real space.

The hypervolume of a 4D ball of radius $r$ is $\frac12\pi^2r^4$ and the hyperarea (actually volume) of the 4D sphere constituting its surface is $2\pi^2r^3$.  Thus the probability density function $f(r)$ of the 4D distribution must satisfy
\[ \int_0^\infty f(r) 2\pi^2r^3 \mathrm d r = 1. \]
One may also define the cumulative probability function
\[ F(R) = \int_0^R f(r) 2\pi^2r^3 \mathrm d r, \]
which is the proportion of the distribution for which $r\le R$.

\subsection{Example: Gaussian}
A Gaussian distribution has $f(r) = Ae^{-r^2/2\sigma^2} = A\prod_{n=1}^4 e^{-x_n^2/2\sigma^2}$, the latter property meaning its projections (distributions of $x_1$ etc.) are also Gaussian with standard deviation $\sigma$.  Normalisation requires
\begin{eqnarray*}
1 &=& \int_0^\infty Ae^{-r^2/2\sigma^2} 2\pi^2r^3 \mathrm d r \\
&=& 2\pi^2 A \int_0^\infty r^3 e^{-r^2/2\sigma^2} \mathrm d r \\
&=& 2\pi^2 A \left[ (-2\sigma^4-\sigma^2r^2)e^{-r^2/2\sigma^2} \right]_{r=0}^\infty \\
&=& 2\pi^2 A (2\sigma^4) \\
\Rightarrow \quad A &=& \frac1{4\pi^2\sigma^4}.
\end{eqnarray*}
Thus for a 4D Gaussian, $ f(r) = \frac{1}{4\pi^2\sigma^4}e^{-r^2/2\sigma^2}$ and
\begin{eqnarray*}
F(R) &=& 2\pi^2 A \left[ (-2\sigma^4-\sigma^2r^2)e^{-r^2/2\sigma^2} \right]_{r=0}^R \\
&=& \frac1{2\sigma^4} \left( (-2\sigma^4-\sigma^2R^2)e^{-R^2/2\sigma^2} + 2\sigma^4 \right) \\
&=& 1-\left(1+\frac{R^2}{2\sigma^2}\right)e^{-R^2/2\sigma^2}.
\end{eqnarray*}
As a check, for $R\ll\sigma$, terms $(R/\sigma)^6\simeq0$ and smaller can be neglected.  Thus $e^{-R^2/2\sigma^2}\simeq1-\frac{R^2}{2\sigma^2}+\frac{R^4}{8\sigma^4}$, so 
$ F(R)\simeq1-(1+\frac{R^2}{2\sigma^2})(1-\frac{R^2}{2\sigma^2}+\frac{R^4}{8\sigma^4})\simeq1-(1-\frac{R^4}{8\sigma^4})=\frac{R^4}{8\sigma^4}. $
The ball of radius $R$ has volume $\frac12\pi^2R^4$, so this suggests a density of approximately $\frac{1}{4\pi^2\sigma^4}$, which is equal to $f(0)=A$ as required.

\subsection{Example: Unit Waterbag}
A unit waterbag distribution has constant 4D density $f(r)=\frac{2}{\pi^2}$ for $r\le1$ and $f(r)=0$ otherwise.
\[ F(R) = \int_0^R \frac{2}{\pi^2} 2\pi^2r^3 \mathrm d r = \int_0^R 4r^3 \mathrm d r = R^4, \]
for $R\le1$ and $F(R)=1$ otherwise.

\section{Exponential Tails}
Beam-beam simulations for eRHIC produce disrupted distributions with tails that look like $e^{-r}$ in a log plot of the 1D projection.  No standard beam distribution that the author knows of has such tails, so this document tries to provide some examples.  The obvious choice is $f(r)=Ae^{-r}$, giving the normalisation condition
\begin{eqnarray*}
1 &=& \int_0^\infty Ae^{-r} 2\pi^2r^3 \mathrm d r \\
&=& 2\pi^2 A \int_0^\infty r^3 e^{-r}\mathrm d r \\
&=& 2\pi^2 A \left[ -(r^3+3r^2+6r+6)e^{-r} \right]_{r=0}^\infty \\
&=& (2\pi^2 A) 6 \\
\Rightarrow \quad A &=& \frac1{12\pi^2}.
\end{eqnarray*}
So the 4D density function is $f(r) = \frac{1}{12\pi^2}e^{-r}$ and
\begin{eqnarray*}
F(R) &=&  2\pi^2 A \left[ -(r^3+3r^2+6r+6)e^{-r} \right]_{r=0}^R \\
&=& \frac16 \left( -(R^3+3R^2+6R+6)e^{-R} + 6 \right) \\
&=& 1-\left(\frac16R^3+\frac12R^2+R+1\right)e^{-R}.
\end{eqnarray*}

\subsection{Computational Generation of Samples}
One way to generate points from these distributions is to select $X\in[0,1]$ uniformly at random and then find $R$ such that $F(R)=X$.  Then a random 4-vector $\mathbf n$ on the sphere with $|\mathbf n|=1$ is required, giving the distribution point as $R\mathbf n$.  This requires $R=F^{-1}(X)$ to be calculated, meaning than the inverse of $F$ should be easily computable.  The inverse of the $F$ in the last section, with the product of a cubic polynomial and an exponential, is probably not even analytically expressible, so other choices may be preferable.

\subsection{Modified Distribution for Easier Computation}
The important feature of $F$ for these distributions is that it decays exponentially at large radii, say $F(R)\simeq1-Be^{-CR}$ and corresponds to a finite constant density at radii near zero (so there isn't a singularity or hole), in other words $F(R)\simeq AR^4$.  One of the simplest choices satisfying both of these conditions appears to be $F(R)=(1-e^{-R})^4$.  Its inverse is easily computable as $F^{-1}(X) = -\ln(1-X^{1/4})$.

To find the phase space density function, note the derivative of $F$ satisfies
\[ F'(r) = f(r) 2\pi^2r^3, \]
so that in this case
\[ f(r) = \frac{F'(r)}{2\pi^2r^3} = \frac{4(1-e^{-r})^3e^{-r}}{2\pi^2r^3}
 = \frac2{\pi^2}\frac{(1-e^{-r})^3e^{-r}}{r^3}. \]
The density near $r=0$ is thus $f(0) = \frac2{\pi^2}$, remembering that $1-e^{-r}\simeq r$ for small $r$ (more precisely, $(1-e^{-r})/r \rightarrow 1$ as $r\rightarrow0$).  This is the same density as the unit waterbag.  In the tails, $f(r)$ is asymptotically proportional to $e^{-r}/r^3$.

Note that both the above and the $f(r)\propto e^{-r}$ distribution are slightly unphysical because $f'(0)\ne0$, so the 4D density is not differentiable at $r=0$.  This does not produce anomalies in most uses, particularly as it is smoothed out by projecting the distribution into 2D or 1D.

\section{Weighted Sampling of Tails}
\newcommand{\rmax}{r_\mathrm{max}}
\newcommand{\rcore}{r_\mathrm{core}}
\newcommand{\alts}[1]{\left\{\begin{array}{ll}#1\end{array}\right.}
Often only the outer tails of a beam distribution will exhibit interesting behaviour, making it wasteful to include the full number of particles in the beam core.  In these cases, better statistics may be obtained by assigning each particle a weight
\[ w(r) = \frac{F'(r)}{G'(r)} = \frac{f(r)}{g(r)}, \]
where $f(r)$ is the desired density and $g(r)$ is that actually generated.  An obvious choice for $g$ would come from $G(r)=r/\rmax$ for $r\le\rmax$ but this samples an infinite phase space density at the origin since $g(r)=\frac{G'(r)}{2\pi^2r^3}=\frac{1}{2\pi^2\rmax r^3}$.  This can be fixed by making $g(r)$ constant (a waterbag) for $r\le\rcore$.  Gluing the two parts of the function together and normalising gives
\[ G'(r) = A\alts{(r/\rcore)^3&r\le\rcore\\
1&r\in[\rcore,\rmax]\\
0&r\ge\rmax}
\quad\Rightarrow\quad
G(r) = A\alts{\frac{\rcore}{4}(r/\rcore)^4&\\
\frac{\rcore}{4}+r-\rcore&\\
\rmax-\frac34\rcore&}, \]
so $A=1/(\rmax-\frac34\rcore)$ and
\[ g(r) = \frac{1}{2\pi^2(\rmax-\frac34\rcore)}\alts{1/\rcore^3&r\le\rcore\\
1/r^3&r\in[\rcore,\rmax]\\
0&r\ge\rmax}. \]
To generate the distribution computationally requires the inverse for $X\in[0,1]$
\[ G^{-1}(X) = \alts{
\rcore\sqrt[4]{(4\frac{\rmax}{\rcore}-3)X} & X\le1/(4\frac{\rmax}{\rcore}-3) \\
(\rmax-\frac34\rcore)X + \frac34\rcore & \mbox{elsewhere.}} \]

\subsection{Asymptotically Exponential Sampling}
A weakness of the above method is that samples for $r>\rmax$ are not generated at all, ignoring $1-F(\rmax)$ of the distribution.  In principle, any infinite tail (with finite integral) could be added to $g(r)$ but $e^{-kr}$ will be used here, since with $k=1$ it matches the exponential tails distributions defined earlier.  Analogously to before,
\[ G'(r) = A\alts{(r/\rcore)^3&r\le\rcore\\
1&r\in[\rcore,\rmax]\\
e^{-k(r-\rmax)}&r\ge\rmax}
\quad\Rightarrow\quad
G(r) = A\alts{\frac{\rcore}{4}(r/\rcore)^4&\\
\frac{\rcore}{4}+r-\rcore&\\
\rmax-\frac34\rcore+\frac1k-\frac1ke^{-k(r-\rmax)}&} \]
and choose $A$ so that $G(\infty)=A(\rmax-\frac34\rcore+\frac1k)=1$.  This gives
\[ g(r) = \frac{1}{2\pi^2(\rmax-\frac34\rcore+\frac1k)}\alts{1/\rcore^3&r\le\rcore\\
1/r^3&r\in[\rcore,\rmax]\\
e^{-k(r-\rmax)}/r^3&r\ge\rmax}. \]
The inverse cumulative distribution is now
\[ G^{-1}(X) = \alts{
\rcore\sqrt[4]{\frac4{\rcore}(X/A)} & X/A\le\frac{\rcore}4 \\
(X/A) + \frac34\rcore & X/A\in[\frac{\rcore}4,\rmax-\frac34\rcore] \\
\rmax - \frac1k\ln(k(1-X)/A) & \mbox{elsewhere,}} \]
where $X/A=(\rmax-\frac34\rcore+\frac1k)X$ and $X\in[0,1]$.

\subsection{Weighting for Exponential Tails Distribution}
\newcommand{\wtail}{w_\mathrm{tail}}
\newcommand{\wmax}{w_\mathrm{max}}
For the distribution $F(R)=(1-e^{-R})^4$, the sample weights are given by
\[ w(r) = \frac{F'(r)}{G'(r)} = \frac{4(1-e^{-r})^3e^{-r}}{A\{(r/\rcore)^3,1,e^{-k(r-\rmax)}\}}. \]
Setting $k=1$ and evaluating $A$ gives
\[ w(r) = \frac{4(1-e^{-r})^3e^{-r}(\rmax-\frac34\rcore+1)}{\{(r/\rcore)^3,1,e^{-(r-\rmax)}\}}. \]
For large $r$ this tends to $\wtail=4e^{-\rmax}(\rmax-\frac34\rcore+1)$ and for small $r$, where $1-e^{-r}\simeq r$, it tends to $\wmax=4\rcore^3(\rmax-\frac34\rcore+1)$.

Since $F(2)\simeq0.559$ and $F(20)\simeq1-8.24\times10^{-9}$, one choice of sensible values would be $\rcore=2$ and $\rmax=20$.  This would give $\wmax=624$ and $\wtail\simeq1.67\times10^{-7}$.  The proportion of samples in the core would be $G(2)\simeq0.0256$ and the proportion in the far tails would be $1-G(20)\simeq0.0513$.

\end{document}