Sampling a Maxwell-Jüttner distribution

We base our method on that described in the Appendix B of an article by Schnittman and Krolik.

The Maxwell-Jüttner distribution, as a function of the Lorentz factor \(\gamma\), reads

\[f(\gamma) = \gamma^2 \beta \exp\left(- \frac {\gamma}{\theta} \right)\]

where \(\theta\) is the temperature divided by \(mc^2\). It is problematic that the change of variable \(\gamma/\theta\) is impossible, because it requires the cumulative distribution function to be computed for every different temperature.

Instead, the “rejection method” makes it possible to choose another function \(g(\gamma)\) such that \(g(\gamma)>f(\gamma)\) everywhere. It can be chosen so that the cumulative distribution function \(G(\gamma)\) is easy to inverse. First, we take a random number \(U_1\) between 0 and 1, and sample the value \(\gamma_1=G^{-1}(U_1)\). Second, we pick another random number \(U_2\), and if \(U_2<f(\gamma_1)/g(\gamma_1)\), we keep the value \(\gamma_1\). Otherwise, we start over to choose another \(U_1\), and so on until a good value is found.

In this particular case, we choose

\[g(\gamma) = \gamma^2 \exp\left(- \frac {\gamma}{\theta} \right)\]

which verifies \(g(\gamma)>f(\gamma)\) and which has the cumulative distribution function

\[G(\gamma) = \int_1^\gamma g(x) dx = 1 - \exp\left[H(\gamma/\theta)-H(1/\theta)\right]\]

where \(H(u) = -u +\ln(1+u+u^2/2)\).

The rejection methods proceeds as

  1. pick a random \(U_1\)

  2. calculate \(\gamma_1=G^{-1}(U_1)=\theta\; H^{-1}[\ln(1-U_1)+H(1/\theta)]\)

  3. pick a random \(U_2\)

  4. select \(\gamma_1\) if \(U_2<\sqrt{1-\gamma_1^{-2}}\), otherwise restart from point 1

Now, to do this, we need to know \(H^{-1}\), which is not easy. We choose to tabulate it in Smilei. For \(X>-\exp(-26)\), we use the series development \(H^{-1}(X) = (-6X)^{1/3}\). For \(X<-\exp(12)\), we use the fit \(H^{-1}(X) = -X + 11.35(-X)^{0.06}\). For all points in between, the function is linearly interpolated in log-log scale over 1000 tabulated values.

Note that the rejection method requires to pick several random numbers if the functions \(f\) and \(g\) differ significantly. This strongly slows the calculation down when the temperature is non-relativistic. For this reason, we fall back to the Maxwell-Boltzmann distribution when \(\theta<0.1\).