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