Sample the Maxwell-Jüttner distribution

Question

The Maxwell-Bolztmann distribution is easy to sample numerically, as it can be expressed directly with a unique variable $mv^2/k_B T$. One simply has to tabulate the cumulative distribution function once, and pick randomly the velocities from its inverse. If the temperature varies, one can simply multiply the final velocities by the correct coefficient.

Now, the Maxwell-Jüttner distribution, which is its relativistic counterpart, seemingly cannot be expressed with a unique variable such as $\gamma/T$ (where $\gamma$ is the Lorentz factor). There is always a separate dependence in momentum and temperature. As a consequence, for each temperature, the cumulative distribution function must be re-calculated.

Is there any known technique to sample this distribution that avoids re-tabulating everything for each temperature?
For example, is it possible to express the cumulative distribution function with a unique variable? Or maybe a technique that involves physical arguments to simplify the sampling?

Note: In this particular case, I am not convinced that this question belongs to scicomp.SE.

Answer 1

I finally found a way that turns out more mathematical than physical.

The idea is from Appendix B of this article.

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. I choose to tabulate it. For $X>-\exp(-26)$, I use the series development $H^{-1}(X) = (-6X)^{1/3}$. For $X<-\exp(12)$, I 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 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, I fall back to the Maxwell-Boltzmann distribution when $\theta<0.1$.

Answer 2

This 2015 paper by Seiji Zenitani¹ suggests that the MJ distribution can be written in spherical coordinates as, $$ f(u)\,\mathrm du\sim\frac{1}{T\,K_2(1/T)}\exp\left[-\frac{\sqrt{1+u^2} }T\right]\,\mathrm du $$ and one can sample using one of the following methods:

1. Using the acceptance criterion $\sqrt{1+u^2}<=T\ln X_1X_2X_3X_4$ for (iid) random numbers $X_i$ (among other things, see the paper)
2. Using a maximum value (e.g., $u_{max}=20u_{therm}$) in the integral for finding the CDF (which likely will still be solved numerically); after that a root-finding algorithm can be used.

The author states on sampling relativistic particles,

In relativistic simulations, it is natural to begin with a relativistic Maxwellian, also known as the Jütter-Synge distribution function. In order to load it, perhaps the Sobol algorithm is the most popular, at least in Monte–Carlo simulation community...To the best of our knowledge, the algorithms for the Jütter-Synge distribution have not been applied to the relativistic shifted-Maxwellian. Several alternative algorithms have been proposed. Swisdak applied a rejection method for a log-concave distribution function. Melzani et al. utilized a numerical cumulative distribution function and cylindrical transformation.

where the "Sobol algorithm" is 1. above (and, again, is detailed in the paper); perhaps these latter references would be of some help as well. The author does suggest thus method is more efficient than acceptance-rejection method.

^{1. http://dx.doi.org/10.1063/1.4919383}