# Sampling Maxwell-Jüttner distribution for non-unity mass and speed of light

I am trying to sample Maxwell-Jüttner distribution using the Sobol method as described in Zenitani Loading relativistic Maxwell distributions in particle simulations (2015). Equation (2) in the paper assumes $$m=c=1$$ and then goes on to describe how to sample the four-velocity. I have a few questions regarding the Sobol algorithm:

1. If $$m$$ and $$c$$, are not taken to be unity, how will that affect the value of generated $$u$$? I can see how this will affect $$f(u)$$, but does $$u$$ also need to be scaled?
2. In order to get $$v$$, which the normal three velocity vector, do we just need to use the relation $$u = \gamma v$$, and express $$v$$ in terms of $$u$$?

I am attaching my distribution plot for reference. • – Kyle Kanos Aug 27 '19 at 19:32
• @KyleKanos Yeah I saw that post, and I don't have any problem as such in sampling. I am just confused about proper normalization, as to how should I scale my u(and therefore v), and f(u) when the mass and speed of light is not taken to be unity. – Prav001 Aug 27 '19 at 22:24
• I linked a related query since it's about the same paper in both. I wasn't trying to suggest it answers your own question. – Kyle Kanos Aug 27 '19 at 22:45

OP's equations are \begin{align} f(\mathbf{u})\,\mathrm{d}^3u&=\frac{N}{4\pi m^2cTK_2(mc^2/T)}\exp\left(-\frac{\gamma mc^2}{T}\right)\,\mathrm{d}^3u\tag{1} \\ \Rightarrow f(u)\,\mathrm{d}u&=\frac{N}{TK_2(1/T)}\exp\left(-\frac{\sqrt{1+u^2}}{T}\right)u^2\,\mathrm{d}u\tag{2} \end{align} where $$\mathbf{u}=\gamma\mathbf{v}$$ is the spatial components of the 4-velocity, $$\gamma$$ the Lorentz factor, $$N=\int f(\mathbf{u})\,\mathrm{d}u$$ is the total number density, $$T$$ the temperature and $$K_2(x)$$ the modified Bessel function. The $$4\pi$$ factor is dropped in (2) due to use of spherical coordinates and it is plain to see that, $$u=\gamma v\implies \gamma=\sqrt{1+u^2/c^2}$$
Since $$u$$ is the velocity, then $$m\neq1$$ should be straight-forward re-insertion. Additionally, since $$u=\gamma v$$, then there should be no changes to $$u$$, only to $$f(u)$$. Lastly, indeed you should invert the relation between $$u$$ and $$v$$ to get the 3-velocity, though the work is done for you in Zenitani's paper (and above): $$v_i=\frac{u_i}{\gamma}=u_i\cdot\left(1+u_i u^i\right)^{-1/2}$$
In my opinion, it may not be the best of ideas to not use $$c=1$$, since using such scaling is intended for convenience (in both notation and in coding) in relativistic settings. Since the Maxwell-Jüttner distribution describes the relativistic distribution of velocities, it probably would be better to stick with $$c=1$$.
• Thanks for the answer. The reason why I am reluctant to use c= 1 is that because this sampling is part of a bigger code, and in other places, we have used c= 3e10. When I tried rescaling f(u) by $\frac{1}{(mc)^3}$ and replacing u by u/c in the definition of gamma, I am not scaled version of my original plot. New plot: imgur.com/WYY6cYA Something tells me I need to rescale u too. – Prav001 Aug 29 '19 at 19:00
• I've always been taught that you should keep units out of codes, rescaling for visualizing only (maybe not even then). $u$ will have same units as $v$, of so I don't think additional factors of $c$ is needed. I don't know why you're scaling $f(u)$ by $1/(mc)^3$, though, that's not what's called for... – Kyle Kanos Aug 29 '19 at 19:22
• The $\frac{1}{(mc)^3}$ factor comes from equation (1). Eqn. 2 is obtained by moving to spherical coordinates and setting m=c=1. So I just rescaled f(u) in Eqn (2) by replacing N by $N * \frac{1}{(mc)^3}$ and $\sqrt {1+u^2}$ as $\sqrt {1+\frac{u^2}{c^2}}$. I think I will stick with m=c=1 to avoid any confusion. – Prav001 Aug 29 '19 at 21:16
• I'm still confused. There's a clear $m^2c$ in the denominator in (1) & nothing else, so you've got a few too many factors floating there. There's also a factor of $mc^2$ in the Bessel fcn & another in the exponent, have you accounted for those? – Kyle Kanos Aug 29 '19 at 21:54
• Sorry for the confusion, but I am not using T, rather I have defined a dimensionless temperature $\theta = \frac{T}{mc^2}$. So $T = mc^2 \theta$ and I can effectively replace T in exponential and the Bessel function in Eqn 1 by $\theta$. So equation 1 becomes $f(u) = \frac{N}{4 \pi (mc)^3 \theta K_{2}(\frac{1}{\theta})}* exp(- \frac{\gamma}{\theta})$ I am using $\theta$ = 0.4 – Prav001 Aug 29 '19 at 22:14