I was trying to code some Monte Carlo calculations regarding x-ray matter interactions and came to the point where I wanted to sample the energy transfer or the scattering angle. However I had some difficulty with turning the Klein-Nishina differential cross-section into a cumulative distribution function. Namely I ended up finding negative CDF values at a photon energy of 10keV.
How do I best find a normalized CDF for Compton scattering that only depends on pre- and post-scatter energy?
What I tried to do is
- start with the Klein-Nishina differential cross-section $\frac{d\sigma}{d\cos(\theta)d\phi} = A \left(\frac{E'}{E}\right)^2\left[\frac{E'}{E} + \frac{E}{E'} – \sin^2(\theta)\right]$
- substitute the $\theta$ terms with the corresponding energy post scatter $E'$ using $E' = \frac{E}{1+(1-\cos(\theta))E/mc^2}$
- "integrate" both sides with respect to $\phi$
- calculate the indefinite integral with respect to $E'$
What I got in the end was (for simplicity expressed with $\epsilon = E'/E$, $k=\frac{mc^2}{E}$ and $l=\frac{1}{1-2/k}$ and up to a constant factor)
$$F(\epsilon)=\frac{\epsilon^3}{3} – k\left(1+\frac{k}{2}\right)\cdot \epsilon^2 + k\cdot ln\left(\frac{l}{\epsilon}\right) + \left(1-2k-2k^2\right) \cdot \epsilon$$
With the normalized CDF I was looking for being $G(\epsilon)=\frac{F(\epsilon)-F(l)}{F(1)-F(l)}$. $l$ was the same as $\epsilon$ at $\theta = \pi$, so I used it as a lower integration bound.
When I threw that into python to plot it for $E=10keV$ over $E' \in [l\cdot E,E]$ I got this:
It looks mostly sensible, but in the beginning it goes below 0 which doesn't really make sense…
I tried looking around for other people solving this problem and found this old paper. And while the function I'm looking for must be in there somewhere, all of the expressions are so obfuscated and the typeface so hard to read that I didn't really manage to extract it. Maybe some of the substitution names are standard and someone else has an easier time with it ^^
Is that procedure not sound? Did I just a wrong lower integration bound? Or is there just another arithmetic error hidden in my notebook?
I also tried to substitute $E'$ for $\theta$, but I (or better said wolfram alpha) couldn't solve the indefinite integral.
Best Answer
The well-known cross section formula for Compton scattering is $$ \frac{d\sigma}{d\Omega} =\frac{\alpha^2(\hbar c)^2}{2s} \left(\frac{\omega'}{\omega}\right)^2 \left(\frac{\omega}{\omega'}+\frac{\omega'}{\omega}-\sin^2\theta\right) $$ where $$ \omega'=\frac{\omega}{1+\frac{\hbar\omega}{mc^2}(1-\cos\theta)} $$
We can integrate $d\sigma$ to obtain a cumulative distribution function. Let $I(\theta)$ be the following integral of $d\sigma$. (The $\sin\theta$ is due to $d\Omega=\sin\theta\,d\theta\,d\phi$) \begin{equation*} I(\theta)= \int \left(\frac{\omega'}{\omega}\right)^2 \left(\frac{\omega}{\omega'}+\frac{\omega'}{\omega}-\sin^2\theta\right) \sin\theta\,d\theta \end{equation*}
The solution is \begin{multline*} I(\theta)=-\frac{\cos\theta}{R^2} +\log\big(1+R(1-\cos\theta)\big)\left(\frac{1}{R}-\frac{2}{R^2}-\frac{2}{R^3}\right) \\ {}-\frac{1}{2R\big(1+R(1-\cos\theta)\big)^2} +\frac{1}{1+R(1-\cos\theta)}\left(-\frac{2}{R^2}-\frac{1}{R^3}\right) \end{multline*} where \begin{equation*} R=\frac{\hbar\omega}{mc^2} \end{equation*}
The cumulative distribution function is $$ F(\theta)=\frac{I(\theta)-I(0)}{I(\pi)-I(0)},\quad 0\le\theta\le\pi $$
To plot $F(E'/E)$, convert $E'/E$ to $\theta$ using Compton's formula. $$ \theta=\arccos\left(\frac{mc^2}{E}-\frac{mc^2}{E'}+1\right) $$