source: trunk/documents/UserDoc/UsersGuides/PhysicsReferenceManual/latex/electromagnetic/muons/pair.tex

\section{Positron - Electron Pair Production by Muons}

Direct electron pair production is one of the most important muon interaction
processes.  At TeV muon energies, the pair production cross section exceeds 
those of other muon interaction processes over a range of energy transfers
between 100 MeV and 0.1$E_{\mu}$.  The average energy loss for pair 
production increases linearly with muon energy, and in the TeV region this 
process contributes more than half the total energy loss rate.

To adequately describe the number of pairs produced, the average energy loss
and the stochastic energy loss distribution, the differential cross section 
behavior over an energy transfer range of 
5~MeV~$\leq \epsilon \leq $~0.1~$\cdot E_{ \mu }$ must be accurately 
reproduced.  This is is because the main contribution to the total cross 
section is given by transferred energies
5 MeV $\leq $ $\epsilon $ $\leq $ 0.01 $\cdot E_{ \mu }$, and because the
contribution to the average muon energy loss is determined mostly in the 
region
$0.001 \cdot E_{ \mu } \leq \epsilon \leq $ 0.1 $\cdot E_{\mu }$ .

For a theoretical description of the cross section, the formulae of 
Ref.~\cite{pair.koko69} are used, along with a correction for finite nuclear 
size \cite{pair.koko71}.  To take into account electron pair production in 
the field of atomic electrons, the inelastic atomic form factor contribution
of Ref. \cite{pair.keln97} is also applied.


\subsection{Differential Cross Section}

\subsubsection{Definitions and Applicability}

In the following discussion, these definitions are used:

\begin{itemize}
\item $m$ and $\mu$ are the electron and muon masses, respectively
\item $E\equiv E_{\mu}$ is the total muon energy, $E = T +\mu $
\item $Z$ and $A$ are the atomic number and weight of the material
\item $\epsilon$ is the total pair energy or, approximately, the muon energy 
  loss $(E-E')$
\item $v = \epsilon/E$
\item $e = 2.718\dots$
\item $A^{\star} = 183$.
\end{itemize}

\noindent
The formula for the differential cross section applies when: 

\begin{itemize}
\item $E_{\mu} \gg  \mu$ ($E \geq$ 2 -- 5 GeV) and 
$E_{\mu} \leq 10^{15}$ -- $10^{17} $ eV.  If muon energies exceed this limit,
the LPM (Landau Pomeranchuk Migdal) effect may become important, depending on 
the material

\item the muon energy transfer $\epsilon$ lies between 
$\epsilon_{\rm min} = 4\,m$ and $\epsilon_{\rm max} = E_{\mu}- \frac{3 \sqrt{e}}{4}\,\mu \,Z^{1/3}$, although the formal lower limit is 
$\epsilon\gg 2\, m$, and the formal upper limit requires $E'_{\mu }\gg\mu$.

\item $Z \leq$ 40 -- 50.  For higher $Z$, the Coulomb correction is important
but has not been sufficiently studied theoretically.
\end{itemize}
                       
\subsubsection{Formulae}

The differential cross section for electron pair production by muons
$\sigma(Z,A,E,\epsilon)$ can be written as :
\begin{equation}
\label{mupair.a} 
\sigma(Z,A,E,\epsilon)=\frac{4}{3\pi}\,\frac{Z(Z+\zeta )}
{A}\, N_{A}\,(\alpha r_{0})^{2}\, \frac{1-v}{\epsilon}
\int_{0}^{\rho_{\rm max}}G(Z,E,v,\rho)\,d\rho ,
\end{equation}

\noindent 
where 
$$G(Z,E,v,\rho ) = \Phi_e + (m/\mu)^2 \Phi_{\mu} , $$
$$\Phi_{e,\mu} = B_{e,\mu} L'_{e,\mu} $$ and
$$\Phi_{e,\mu} = 0 \quad {\rm whenever} \quad \Phi_{e,\mu} < 0 . $$

\noindent 
$B_{e}$ and $B_{\mu}$ do not depend on $Z,A$, and are given by

$$B_{e}=[(2+\rho^{2})(1+\beta)+\xi(3+\rho^{2})]
\ln\left(1+\frac{1}{\xi}\right)+\frac{1-\rho^{2}-\beta}{1+\xi}-(3+\rho^{2}) ; $$
$$B_{e}\approx
\frac{1}{2\xi}\,[(3-\rho^{2})+2\beta(1+\rho^{2})]\quad {\rm 
for}\quad\xi\geq 10^{3};$$

$$B_{\mu}=\left[(1+\rho^{2})
\left(1+\frac{3\beta}{2}
\right)-\frac{1}{\xi}(1+2\beta)
(1-\rho^{2})\right]\ln(1+\xi )
$$
$$
+\frac{\xi(1-\rho^{2}-\beta)}{1+\xi}+
(1+2\beta)(1-\rho^{2});$$
$$B_{\mu}\approx\frac{\xi}{2}\,[(5-\rho^{2})+\beta(3+\rho^{2})]\quad{\rm 
for}\quad\xi\leq10^{-3}
;$$

\noindent 
Also,
% $B_{e},\:B_{\mu}$ do not depend on $\,Z,\:A$.

$$\xi=\frac{\mu^{2} v^{2}}{4m^{2}}\, \frac{(1-\rho^{2})}{(1-v)};\quad
\beta=\frac{v^{2}}{2(1-v)};$$

$$L'_{e}=\ln\frac{A^{*}Z^{-1/3}\sqrt{(1+\xi)(1+Y_{e})}}
{1+\displaystyle\frac{2m\sqrt{e}A^{*}Z^{-1/3}(1+\xi)(1+Y_{e})}{Ev(1-\rho^{2})}}
$$
$$
-\frac{1}{2}\ln\left[1+\left(\frac{3mZ^{1/3}}{2\mu}\right)^{2}(1+\xi)(1+Y_{e})\right];
$$

$$L'_{\mu}=\ln\frac{(\mu/m)A^{*}Z^{-1/3}\sqrt{(1+1/\xi)(1+Y_{\mu})}}
{1+\displaystyle\frac{2m\sqrt{e}A^{*}Z^{-1/3}(1+\xi)(1+Y_{\mu})}{Ev(1-\rho^{2})}}
$$
$$
-\ln\left[\frac{3}{2}\,Z^{1/3}\sqrt{(1+1/\xi)(1+Y_{\mu})}\right].
$$
%
For faster computing, the expressions for $L'_{e,\mu }$ are
further algebraically transformed.  The functions $L'_{e,\mu }$ include the 
nuclear size correction \cite{pair.koko71} in comparison with parameterization
\cite{pair.koko69} :
%
$$Y_{e}=\frac{5-\rho^{2}+4\,\beta\,(1+\rho^{2})}
{2(1+3\beta)\ln(3+{1}/{\xi})-\rho^{2}-2\beta(2-\rho^{2})} ;$$

$$Y_{\mu}=\frac{4+\rho^{2}+3\,\beta\,(1+\rho^{2})}
{(1+\rho^{2})(\frac{3}{2}+2\beta)\ln(3+\xi)+1-\frac{3}{2}\,\rho^{2}};$$
%$Y_{e,\mu}$~--~ approximations (\cite{pair.koko69}); there is
%a possibility that they will be improved , as well as corrections (second
% terms) in $L_{e,\mu}$.
% 
$$\rho_{\rm max}=[1-6\mu^{2}/E^{2}(1-v)]\sqrt{1-4m/Ev}.$$

\subsubsection{Comment on the Calculation of the Integral $\int\!d\rho$ in 
Eq.~\ref{mupair.a}}
 
The integral $\int\limits_{0}^{\rho_{\max}}G(Z,E,v,\rho)~d\rho$ is computed 
with the substitutions:
\begin{eqnarray*}
t&=& \ln(1-\rho),\\
1 - \rho &= &\exp(t),\\
1 + \rho& =& 2 - \exp(t),\\
1 - \rho^{2} &=&e^{t}\,(2-e^{t}). 
\end{eqnarray*}

\noindent 
After that,
\begin{equation}
\label{mupair.b}
\int_{0}^{\rho_{\rm max}}G(Z,E,v,\rho)~d\rho = 
\int_{t_{\rm min}}^{0}G(Z,E,v,\rho )\,e^{t}\,dt , 
\end{equation}

\noindent
where

$$t_{\rm
min}=\ln\frac{\displaystyle\frac{4m}{\epsilon}+\frac{12\mu^{2}}{EE'}\left(1-\frac{4m}
{\epsilon}\right)}
{\displaystyle 1+\left(1-\frac{6\mu^{2}}{EE'}\right)\sqrt{1-\frac{4m}{\epsilon}}}.
$$

To compute the integral of Eq.~\ref{mupair.b} with an accuracy better than 
0.5\%, Gaussian quadrature with $N=8$ points is sufficient.

The function $\zeta(E,Z)$ in Eq.~\ref{mupair.a} serves to take into account 
the process on atomic electrons (inelastic atomic form factor contribution).  
To treat the energy loss balance correctly, the following approximation, 
which is an algebraic transformation of the expression in 
Ref.~\cite{pair.keln97}, is used:

$$\zeta(E,Z)=\frac{\displaystyle 0.073\ln\frac{E/\mu}{ 
1+\gamma_{1}Z^{2/3}E/\mu}
-0.26}
{\displaystyle 0.058\ln\frac{E/\mu}{1+\gamma_{2}Z^{1/3}E/\mu}
-0.14}; $$
$$\zeta(E,Z)=0\quad\mbox{if the numerator is negative.}$$

\noindent
For E $\leq 35\,\mu, ~ \zeta(E,Z)=0$.  Also
$\gamma_{1}= 1.95\cdot 10^{-5} $ and $\gamma_{2}= 5.30\cdot 10^{-5} $.

The above formulae make use of the Thomas-Fermi model which is not good
enough for light elements.  For hydrogen ($Z = 1$) the following parameters
must be changed: \\
$A^{*}=183 ~\Rightarrow ~202.4;$\\
$\gamma_{1}= 1.95\cdot 10^{-5}  ~\Rightarrow ~ 4.4\cdot 10^{-5};$\\
$\gamma_{2}= 5.30\cdot 10^{-5}  ~\Rightarrow  ~4.8\cdot 10^{-5}.$\\

\subsection{Total Cross Section and Restricted Energy Loss}

If the user's cut for the energy transfer $\epsilon_{\rm cut}$ is
greater than $\epsilon_{\min}$, the process is represented by
continuous restricted energy loss for interactions with
$\epsilon\le\epsilon_{\rm cut}$, and discrete collisions with
$\epsilon>\epsilon_{\rm cut}$. Respective values of the total cross
section and restricted energy loss rate are defined as:
$$
\sigma_{\rm tot}=\int_{\epsilon_{\rm cut}}^{\epsilon_{\max}}
\sigma(E,\epsilon)\,d\epsilon;\quad
(dE/dx)_{\rm restr}=\int_{\epsilon_{\min}}^{\epsilon_{\rm cut}}
\epsilon\,\sigma(E,\epsilon)\,d\epsilon.
$$ 
For faster computing, $\ln\epsilon$ substitution and Gaussian quadratures are 
used.

\subsection{Sampling of Positron - Electron Pair Production}

The e$^+$e$^-$ pair energy $\epsilon_P$, is found numerically by solving the
equation
\begin{equation}
\label{mupair.c}
 P = \int_{\epsilon_P}^{\epsilon_{\rm max}} \sigma (Z,A,T,\epsilon) d\epsilon
 \quad /  \int_{cut}^{\epsilon_{\rm max}} \sigma (Z,A,T,\epsilon) d\epsilon
\end{equation}
or
\begin{equation}
\label{mupair.d}
  1 - P = \int_{cut}^{\epsilon_P} \sigma (Z,A,T,\epsilon) d\epsilon 
 \quad /  \int_{cut}^{\epsilon_{\rm max}} \sigma (Z,A,T,\epsilon) d\epsilon
\end{equation}

To reach high sampling speed, solutions of Eqs.~\ref{mupair.c}, \ref{mupair.d}
are tabulated at initialization time.  Two 3-dimensional tables (referred to 
here as A and B) of $\epsilon_{P}(P,T,Z)$ are created, and then interpolation
is used to sample $\epsilon_P$.

\noindent
The number and spacing of entries in the table are chosen as follows:
\begin{itemize}
\item a constant increment in $\ln T$ is chosen such that there are four 
points per decade in the range $T_{\rm min}- T_{\rm max}$.  The default range 
of muon kinetic energies in Geant4 is $T=1\:{\rm GeV}-1000\:{\rm PeV}$.

\item a constant increment in $\ln Z$ is chosen.  The shape of the sampling 
distribution does depend on $Z$, but very weakly, so that eight points in the 
range $1\leq Z\leq 128$ are sufficient.  There is practically no dependence 
on the atomic weight $A$.

\item for probabilities $P \leq 0.5$, Eq.~\ref{mupair.c} is used and Table~A 
is computed with a constant increment in $\ln P$ in the range 
$10^{-7}\leq P \leq 0.5$.  The number of points in $\ln P$ for Table~A is
about 100.

\item for $P \geq 0.5$, Eq.~\ref{mupair.d} is used and Table~B is computed
with a constant increment in $\ln(1-P)$ in the range 
$10^{-5} \leq (1-P) \leq 0.5$.  In this case 50 points are sufficient.
\end{itemize}

\noindent 
The values of $\ln (\epsilon_{P}-cut$) are stored in both Table~A and Table~B. 


To create the ``probability tables" for each $(T, Z)$ pair, the following 
procedure is used:

\begin{itemize}
\item a temporary table of $\sim$ 2000 values of 
$\epsilon \cdot \sigma(Z,A,T,\epsilon)$ is constructed with a constant
increment ($\sim$ 0.02) in $\ln \epsilon$ in the range 
$(cut, \epsilon_{\max})$.  $\epsilon$ is taken in the middle of the
corresponding bin in $\ln \epsilon$.

\item the accumulated cross sections
$$\sigma_{1}=\int_{\ln\epsilon}^{\ln\epsilon_{\max}} 
 \epsilon \, \sigma (Z,A,T,\epsilon)\, d (\ln\epsilon)  $$
and
 $$\sigma_{2}=\int_{\ln(cut)}^{\ln\epsilon} 
 \epsilon \, \sigma (Z,A,T,\epsilon)\, d (\ln\epsilon )  $$

are calculated by summing the temporary table over the values above 
$\ln \epsilon$ (for $\sigma_1$) and below $\ln \epsilon$ (for $\sigma_2$)
and then normalizing to obtain the accumulated probability functions.

\item finally, values of $\ln(\epsilon_{P} - cut $) for corresponding values
of $\ln P$ and $\ln (1-P)$ are calculated by linear interpolation of the 
above accumulated probabilities to form Tables A and B.  The monotonic 
behavior of the accumulated cross sections is very useful in speeding up
the interpolation procedure. 

\end{itemize}

The random transferred energy corresponding to a probability $P$, is then 
found by linear interpolation in $\ln Z$ and $\ln T$, and a cubic 
interpolation in $\ln P$ for Table A or in $\ln (1-P)$ for Table B.
For $P \leq 10^{-7}$ and $(1-P) \leq 10^{-5}$, linear extrapolation using
the entries at the edges of the tables may be safely used.  Electron pair 
energy is related to the auxiliary variable
$x = \ln (\epsilon_{P} - cut)$ found by the trivial interpolation
$\epsilon_{P} = e^{x} + cut$.

Similar to muon bremsstrahlung (section \ref{secmubrem}), this sampling 
algorithm does not re-initialize the tables for user cuts greater than 
$cut_{min}$.  Instead, the probability variable is redefined as
$$ P' = P \sigma_{\rm tot}(cut_{user}) / \sigma_{\rm tot}(cut_{min}),$$
and $P'$ is used for sampling.

In the simulation of the final state, the muon deflection angle (which is
of the order of $m/E$) is neglected. The procedure for sampling the energy 
partition between $e^+$ and $e^-$ and their emission angles is similar to 
that used for the $\gamma \to e^+\,e^-$ conversion.

\subsection{Status of this document}

12.10.98 created by R.Kokoulin and A.Rybin\\
18.05.00 edited by S.Kelner, R.Kokoulin, and A.Rybin\\
27.01.03 re-written by D.H. Wright

\begin{latexonly}

\begin{thebibliography}{99}
 
\bibitem{pair.koko69} R.P.Kokoulin and A.A.Petrukhin,
               Proc. 11th Intern. Conf. on Cosmic Rays, Budapest, 1969
               [Acta Phys. Acad. Sci. Hung.,{\bf  29, Suppl.4},
               p.277, 1970].

\bibitem{pair.koko71} R.P.Kokoulin and A.A.Petrukhin,
               Proc. 12th Int. Conf. on Cosmic Rays, Hobart, 1971,
               {\bf vol.6}, p.2436.

\bibitem{pair.keln97} S.R.Kelner, Phys. Atomic Nuclei, 
  {\bf 61} (1998) 448.

\end{thebibliography}

\end{latexonly}

\begin{htmlonly}

\subsection{Bibliography}
 
\begin{enumerate}
\item R.P.Kokoulin and A.A.Petrukhin,
      Proc. 11th Intern. Conf. on Cosmic Rays, Budapest, 1969
      [Acta Phys. Acad. Sci. Hung.,{\bf  29, Suppl.4}, p.277, 1970].

\item R.P.Kokoulin and A.A.Petrukhin,
      Proc. 12th Int. Conf. on Cosmic Rays, Hobart, 1971, {\bf vol.6}, p.2436.

\item S.R.Kelner, Phys. Atomic Nuclei, {\bf 61} (1998) 448.
\end{enumerate}

\end{htmlonly}