\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}