# Sample the Maxwell-Jüttner distribution

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.

• I found that this article proposes a solution in Appendix A. If someone has tested it, I'd appreciate some comments, or an answer :) Otherwise, I'll answer myself when I have done the test. – fffred Dec 20 '16 at 17:03
• Seems to me you should be able to use the Newton-Raphson method (I briefly outline it here if you're not familiar with it). – Kyle Kanos Dec 23 '16 at 12:29
• @Kyle, if I understand correctly, you're saying I should use the Newton-Raphson method to find the roots of the CDF. However, I do not have an analytical form of the CDF. – fffred Dec 24 '16 at 11:20

This 2015 paper by Seiji Zenitani1 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

• The Sobol method appears the most interesting because there is no need of tabulated values. This is particularly important when many values of the temperatures are needed. Its efficiency, according to the author, still drops heavily in the non-relativistic case sue to the rejection method. But that can be avoided by using a non-relativistic Maxwellian. – fffred Dec 28 '16 at 12:53

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

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$.