Livermore Triple Gamma Conversion

The class G4BoldyshevTripletModel was developed to simulate the pair production by linearly polarized gamma rays on electrons For the angular distribution of electron recoil we used the cross section by Vinokurov and Kuraev [VK72][VK73] using the Borsellino diagrams in the high energy For energy distribution for the pair, we used Boldyshev [VFBP94] formula that differs only in the normalization from Wheeler-Lamb. The cross sections include a cut off for momentum detections [MLI11].

Method

The first step is sample the probability to have an electron recoil with momentum greater than a threshold define by the user (by default, this value is \(p_{0}=1\) in units of \(mc\)). This probability is

\[\sigma(p\geq p_{0})=\alpha r_{0}^{2}\left(\frac{82}{27}-\frac{14}{9}\ln X_{0}+\frac{4}{15}X_{0}-0.0348X_{0}^{2}+0.008X_{0}^{3}-...\right)\]

\[X_{0}=2\left(\sqrt{p_{0}^{2}+}-1\right).\]

Since that total cross section is \(\sigma=\alpha r_{0}^{2}\left(\frac{28}{4} \ln 2E_{\gamma}-\frac{218}{27}\right)\), if a random number is \(\xi\geq\sigma(p\geq p_{0})/\sigma\) we create the electron recoil, otherwise we deposited the energy in the local point.

Azimuthal Distribution for Electron Recoil

The expression for the differential cross section is composed of two terms which express the azimuthal dependence as follows:

\[d\sigma=d\sigma^{(t)}-Pd\sigma^{(l)}\cos(2\varphi)\]

Where both \(d\sigma_{(t)}\) and \(d\sigma_{(l)}\) are independent of the azimuthal angle, \(\varphi\), referred to an origin chosen in the direction of the polarization vector \(\vec{P}\) of the incoming photons.

Monte Carlo Simulation of the Asymptotic Expression

In this section we present an algorithm for Monte Carlo simulation of the asymptotic expressions calculate by Vinokurov et.al. [VK72][VK73].

We must generate random values of \(\theta\) and \(\varphi\) distributed with probability proportional to the following function \(f(\theta,\varphi\)), for \(\theta\) restricted inside of its allowed interval value [VFBP94] (\(0\), or \(\theta_{max}(p_0)\)):

(14)\[f(\theta, \varphi) = \frac{\sin \theta}{\cos^3 \theta} (F_1(\theta) - {\rm P} \cos (2 \varphi) F_P (\theta) )\]

\[F_1 ( \theta ) = 1 - \frac{1-5\cos^2 \theta}{\cos \theta} \ln ( \cot (\theta / 2) )\]

\[F_P (\theta) = 1 - \frac{\sin^2 \theta}{\cos\theta} \ln ( \cot ( \theta / 2) )\]

As we will see, for \(\theta<\pi/2\), \(F_{1}\) is several times greater than \(F_{P}\), and since both are positive, it follows that \(f\) is positive for any possible value of \(P\) (\(0\leq P\leq1\)).

Since \(F_{1}\) is the dominant term in expression (14), it is more convenient to begin developing the algorithm of this term, belonging to the unpolarized radiation.

Algorithm for Non Polarized Radiation

The algorithm was described in Ref.[DI09]. We must generate random values of \(\theta\) between \(0\) and \(\theta_{max}=\arccos\left( \frac{E_1-mc^2}{p_0} + mc^2\frac{E_1+mc^2}{E_{\gamma}p_0}\right) \), \(E_1 = \sqrt{p_0^2 + (mc^2)^2}\) distributed with probability proportional to the following function \(f_{1}(\theta\)):

\[\begin{split}f_{1}(\theta) &= \frac{\sin(\theta)}{\cos^{3}(\theta)}\left(1-\frac{1-5\cos^{2}(\theta)}{\cos(\theta)}\ln (\cot(\theta/2))\right)\\ &= \frac{\sin(\theta)}{\cos^{3}(\theta)}\times F_{1}\left(\theta\right)\end{split}\]

By substitution \(\cos(\theta/2)=\sqrt{\frac{1+\cos\theta}{2}}\) and \(\sin(\theta/2)=\sqrt{\frac{1-\cos\theta}{2}}\) , we can write:

\[\ln(\cot(\theta/2)) = \frac{1}{2}\ln\left(\frac{1+\cos\theta}{1-\cos\theta}\right)\]

In order to simulate the \(f_{1}\) function, it may be decomposed in two factors: the first, \(\sin(\theta)/\cos^{3}(\theta)\), easy to integrate, and the other, \(F_{1}(\theta)\), which may constitute a reject function, on despite of its \(\theta=0\) divergence. This is possible because they have very low probability. On other hand, \(\theta\) values near to zero are not useful to measure polarization because for those angles it is very difficult to determine the azimuthal distribution (due to multiple scattering).

Then, it is possible to choose some value of \(\theta_{0}\), small enough that it is not important that the sample is fitted rigorously for \(\theta<\theta_{0}\), and at the same time \(F_{1}(\theta_{0})\) is not too big.

Modifying \(F_{1}\) so that it is constant for \(\theta\leq\theta_{0}\), we may obtain an adequate reject function. Doing this, we introduce only a very few missed points, all of which lie totally outside of the interesting region.

Expanding \(F_{1}\) for great values of \(\theta\), we see it is proportional to \(cos^{2}\theta\):

\[F_1 (\theta) \to \frac{14}{3} \cos^2 \theta \left( 1+\frac{33}{35} \cos^2 \theta + \ldots \right) , \quad{\rm if} \; \theta \to \pi / 2\]

Thus, it is evident that \(F_{1}\) divided by \(\cos^{2}(\theta\)) will be a better reject function, because it tends softly to a some constant value (\(14/3=4,6666...\)) for large \(\theta\), whereas its behavior is not affected in the region of small \(\theta\), where \(\cos(\theta)\rightarrow1\).

It seems adequate to choose \(\theta_{0}\) near \(5^{\circ}\), and, after some manipulation looking for round numbers we obtain:

\[\frac{F_{1}\left(4.47^{\circ}\right)}{\cos^{2}\left(4.47^{\circ}\right)}\cong14.00\]

Finally we define a reject function:

\[\begin{split}\begin{array}{lllll} r(\theta) & = & \frac{1}{14}\frac{F_{1}(\theta)}{\cos^{2}(\theta)}=\frac{1}{14\cos^{2}(\theta)}\\ & & \left(1-\frac{1-5\cos^{2}(\theta)}{2\cos(\theta)}\ln\left(\frac{1+\cos\theta}{1-\cos\theta}\right)\right) & ; & {\rm for}\;\theta\ge4.47^{0}\\ r\left(\theta\right) & = & 1 & ; & {\rm for} \; \theta\le4.47^{0} \end{array}\end{split}\]

Now we have a probability distribution function (PDF) for \(\theta\), \(p(\theta)=Cf_{1}(\theta\)), expressed as a product of another PDF, \(\pi(\theta\)), by the reject function:

\[p\left(\theta\right)=Cf_{1}\left(\theta\right)\cong C^{'}\pi\left(\theta\right)r\left(\theta\right)\]

where \(C\) is the normalization constant belonging to the function \(p(\theta\)).

One must note that the equality between \(C\sim f_{1}(\theta\)) and \(C^{'}\pi(\theta)r(\theta\)) is not exact for small values of \(\theta\), where we have truncated the infinity of \(F_{1}(\theta\)); but this can not affect appreciably the distribution because \(f_{1}\rightarrow0\) there. Now the PDF \(\pi(\theta\)) is:

\[\pi(\theta)=C_{\pi}\;\frac{14\sin(\theta)}{\cos(\theta)}\]

From the normalization, the constant \(C_{\pi}\) results:

\[C_{\pi}=\frac{1}{14\int_{0}^{\theta_{max}}\frac{\sin(\theta)}{\cos(\theta)}d\theta}=\frac{-1}{14\ln\left(\cos(\theta_{max})\right)}=\frac{1}{7}\ln\left(\frac{\omega}{4m}\right)\]

And the relation with \(C\) is given by:

\[C=\frac{1}{\int_{0}^{\theta_{max}}f_{1}(\theta)d\theta}\cong C'C_{\pi}\]

Then we obtain the cumulative probability by integrating the PDF \(\pi(\theta\)):

\[P_{\pi}=\int_{0}^{\theta}\pi(\theta')d\theta'=\frac{-14\ln(\cos(\theta))}{7\ln\left(\frac{\omega}{4m}\right)}=\frac{2\ln(\cos(\theta))}{\ln\left({4m/\omega}\right)}\]

Finally for the Monte Carlo method we sample a random number \(\xi_{1}\) (between 0 and 1), which is defined as equal to \(P_{\pi}\) , and obtain the corresponding \(\theta\) value:

\[\xi_1 = \frac{2 \ln(\cos \theta) }{\ln (4m / \omega)} = \frac{\ln (\cos \theta)}{\ln(\cos(\theta_{\rm max}))}\]

Then,

\[\theta=\arccos\left(\left(\frac{4{\kern1pt }m}{\omega}\right)^{\frac{\xi_{1}}{2}}\right)\]

Another random number \(\xi_{2}\) is sampled for the reject process: the \(\theta\) value is accepted if \(\xi_{2}\leq r(\theta\)), and reject in the contrary.

For \(\theta\leq4.47^{\circ}\) all values are accepted. It happens automatically without any modification in the algorithm previously defined (it is not necessary to define the truncated reject function for \(\theta<\theta_{0})\).

Algorithm for Polarized Radiation

The algorithm was also described in Ref.[DI09]. As we have seen, the azimuthal dependence of the differential cross section is given by the expressions and:

\[f(\theta,\varphi) = \frac{\sin \theta}{\cos^3 \theta} (F_1(\theta) - {\rm P} \cos (2\varphi) F_P (\theta))\]

\[F_P (\theta) = 1 - \frac{\sin^2 \theta}{\cos \theta} \ln (\cot (\theta / 2))\]

We see that \(F_{P}\) tends to \(1\) at \(\theta=0\), decreases monotonically to \(0\) as \(\theta\) goes to \(\pi/2\).

Furthermore, the expansion of \(F_{P}\) for \(\theta\) near \(\pi/2\) shows that it is proportional to \(\cos^{2}(\theta)\), in virtue of which \(F_{P}/\cos^{2}(\theta)\) tends to a non null value, \(2/3\). This value is exactly 7 times the value of \(F_{1}/\cos^{2}(\theta)\).

This suggests applying the combination method, rearranging the whole function as follows:

\[f(\theta,\varphi)=\tan(\theta)\frac{{F}_{1}(\theta)}{\cos^{2}(\theta)}\left(1-\cos(2\varphi)P\frac{F_{P}(\theta)}{F_{1}(\theta)}\right)\]

and the normalized PDF \(p(\theta,\varphi)\):

\[p(\theta,\varphi)=Cf(\theta,\varphi)\]

where \(C\) is the normalization constant

\[\frac{1}{C}=\int_{0}^{\theta_{\max}}\int_{0}^{2\pi}f(\theta,\varphi)\, d\varphi d\theta\]

Taking account that \(\int_{0}^{2{\kern1pt }\pi}\cos(2\varphi)\; d\varphi=0\), then:

\[\frac{1}{C}=2\pi\int_{0}^{\theta_{\max}}\tan(\theta)\frac{F_{1}(\theta)}{\cos^{2}(\theta)}d\theta\]

On the other hand the integration over the azimuthal angle is straightforward and gives:

\[q(\theta)=\int_{0}^{2\pi}p(\theta,\varphi)d\varphi=2\pi C\tan(\theta)\frac{F_{1}(\theta)}{\cos^{2}(\theta)}\]

and \(p(\varphi/\theta)\) is the conditional probability of \(\varphi\) given \(\theta\):

\[\begin{split}p(\varphi/\theta) &= \frac{p(\theta,\varphi)}{q(\theta)}=\frac{1}{2\pi C\tan(\theta)\frac{F_{1}(\theta)}{\cos^{2}(\theta)}}C\frac{\sin(\theta)}{\cos^{3}(\theta)}F_{1}(\theta)\left(1-\cos(2\varphi)P\frac{F_{P}(\theta)}{F_{1}(\theta)}\right)\\ &= \frac{1}{2\pi}\left(1-\cos(2\varphi)P\frac{F_{P}(\theta)}{F_{1}(\theta)}\right)\end{split}\]

Now the procedure consists of sampling \(\theta\) according the PDF \(q(\theta\)); then, for each value of \(\theta\) we must sample \(\varphi\) according to the conditional PDF \(p(\varphi/\theta\)).

Knowing that \(F_{1}\) is several times greater than \(F_{P}\), we can see that P \(F_{1}/F_{P} \ll 1\), and thus \(p(\varphi/\theta\)) maintains a nearly constant value slightly diminished in some regions of \(\varphi\). Consequently the \(\varphi\) sample can be done directly by the rejecting method with high efficiency.

On the other hand, \(q(\theta\)) is the same function \(p(\theta\)) given by , that is the PDF for unpolarized radiation, \(q(\theta)\cong C'\pi(\theta)r(\theta)\), so we can sample \(\theta\) with exactly the same procedure, specified as follows:

1. We begin sampling a random number \(\xi_{1}\) and obtain \(\theta\) from :

   \[\theta=\arccos\left(\left(\frac{4m}{\omega}\right)^{\frac{\xi_{1}}{2}}\right)\]

2. Then we sample a second random number \(\xi_2\) and accept the values of \(\theta\) if \(\xi_{2}\leq r(\theta)\), where \(r(\theta\)) is the same expression defined before:

   \[r(\theta) = \frac{1}{14\, \cos^2 \theta}\,\left[1-\frac{1- 5\cos^2 \theta}{2 \cos \theta} \ln \left(\frac{1+\cos \theta}{1- \cos \theta}\right)\right]\]

   For \(\theta\geq4.47^{\circ}\) and for \(\theta\leq4.47^{\circ}\) all values are accepted.

3. Now we sample \(\varphi\). According to the reject method, we sample a third random number \(\xi_{3}\) (which is defined as \(\varphi/2\pi\)) and evaluate the reject function (which is essentially):

   \[r_{\theta}(\xi_{3})=\frac{1}{2\pi}\left(1-\cos\left(4\pi\xi_{3}\right)P\frac{F_{P}\left(\theta\right)}{F_{1}\left(\theta\right)}\right)\]
   \[=\,\frac{1}{2\pi}\left(1-\cos(4\pi\xi_{3})P\frac{\cos\theta-\sin^{2}\theta\ln\left(\cot\left(\frac{\theta}{2}\right)\right)}{\cos\theta-\left(1-5\cos^{2}\theta\right)\ln\left(\cot\left(\frac{\theta}{2}\right)\right)}\right)\]

4. Finally, with a fourth random number \(\xi_{4}\) , we accept the values of \(\varphi=2\pi\xi_{4}\) if \(\xi_{4}\leq r_{\theta}(\xi_{3})\).

Sampling of Energy

For the electron recoil we calculate the energy from the maximum momentum that can take according with the \(\theta\) angle

\[E_r = mc^2 \frac{\left( S + (mc^2)^2\right)}{D2}\]

where

\[\begin{split}\begin{array}{lll} S &=& mc^2 \left( 2 E_{gamma} + mc^2\right) \\ D2 &=& 4 S mc^2 + \left( S - (mc^2)^2\right)^2 \sin^2(\theta) \end{array}\end{split}\]

The remnant energy is distributed to the pair according to the Boldyshev formula [VFBP94] (\(x\) is the fraction of the positron energy):

\[2\pi \frac{d^2\sigma}{dx d\phi}= 2 \alpha r_0^2 \left\lbrace \left[ 1-2x\left( 1-x\right) \right] J_1(p_0) + 2x\left(1-x \right)\left[ 1 - P \cos(\phi)\right] J_2(p_0)\right\rbrace\]

\[J_1(p_0)= 2 \left( t\frac{\cosh(t)}{\sinh(t)}-ln(2 \sinh(t))\right)\]

\[J_2(p_0)= -\frac{2}{3} \ln (2 \sinh(t)) + t \frac{\cosh(t)}{\sinh(t)} + \frac{\sinh(t) - t \cosh^3(t)}{3 \sinh^3(t)} , \hspace{1cm} \sinh(2t) = p_0\]

This distribution can by written like a PDF for \(x\):

\[P(x)=N\left( 1-Jx(1-x)\right)\]

where \(N\) is a normalization constant and \(J=(J_1 - J_2)/J_1\). Solving for \(x\) (\(\xi\) is a random number):

\[x = \frac{c_1^{1/3}}{2 J}+\frac{J - 4}{2 c_1^{1/3}} + \frac{1}{2}\]

\[\begin{split}\begin{array}{llll} c_1 &=& \left( -6 + 12 r_n + J + 2 a\right) J^2 \\ a &=& \left( \frac{16 -3 J - 36 r_n + 36 J r_n^2 + 6 r_n J^2}{J}\right) \\ r_n &=& \xi\left(1 - \frac{J}{6}\right) \end{array}\end{split}\]

Bibliography

DI09(1,2)

G.O. Depaola and M.L. Iparraguirre. Nucl. Instr. Meth. A, 611():84, 2009.

MLI11

G.O. Depaola M.L. Iparraguirre. The European Physical Journal C, 71():1778, 2011.

VFBP94(1,2,3)

N.P. Merenkov V.F. Boldyshev, E.A. Vinokurov and Yu.P. Peresunko. Phys. Part. Nucl., 25():292, 1994.

VK72(1,2)

E.A. Vinokurov and E.A. Kuraev. Zh. Eksp. Teor. Fiz., 63():1142, 1972. in Russian.

VK73(1,2)

E.A. Vinokurov and E.A. Kuraev. Sov. Phys. JETP, 36:602, 1973. in Russian.