Loading Relativistic Maxwell Distributions in Particle Simulations

Numerical algorithms to load relativistic Maxwell distributions in particle-in-cell (PIC) and Monte-Carlo simulations are presented. For stationary relativistic Maxwellian, the inverse transform method and the Sobol algorithm are reviewed. To boost particles to obtain relativistic shifted-Maxwellian, two rejection methods are proposed in a physically transparent manner. Their acceptance efficiencies are ${\approx}50\%$ for generic cases and $100\%$ for symmetric distributions. They can be combined with arbitrary base algorithms.


I. INTRODUCTION
Because of an increasing demand in high-energy astrophysics, numerical modeling of relativistic kinetic plasmas has been growing in importance. To date, many simulations on relativistic kinetic processes have been performed, such as the Rankine-Hugoniot problem across a relativistic shock 3 , magnetic reconnection and kinetic instabilities 15 in a relativistically hot current sheet, 4,5 and the kinetic Kelvin-Helmholtz instability in a relativistic flow shear 1 . In these simulations, one has to carefully set up ultrarelativistic bulk flows and/or relativistically hot plasmas in their rest frame. Loading velocity distribution function, i.e., initializing particle velocities by using random variables according to a relativistic distribution function, is essentially important.
In nonrelativistic particle simulations, it is quite natural to begin with a Maxwell-Boltzmann distribution (Maxwellian in short). To load the Maxwellian, the Box-Muller algorithm is widely used. 2 One can easily initialize a distribution with a bulk drift velocity, by applying an offset to the particle velocities.
In relativistic simulations, it is natural to begin with a relativistic Maxwellian, also known as the Jütter-Synge distribution function. 6,14 In order to load it, perhaps the Sobol 12 algorithm is the most popular, at least in Monte-Carlo simulation community. The algorithm was formally proposed by Sobol 12 in a Russian proceeding. Its key results are outlined in Pozdnyakov et al. 10,11 . Meanwhile, it is not so clear how to initialize particles according to the relativistic shifted-Maxwellian or moving population of other distributions. 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 13 applied a rejection method for a log-concave distribution function. Melzani et al. 9 utilized a numerical cumulative distribution function and cylindrical transformation.
In this research note, we describe numerical methods to load relativistic Maxwellians in particle simulations. We first describe two base algorithms to load stationary relativistic Maxwellian, the inverse transform method and the Sobol method. 12 Next we apply the Lorentz transformation to obtain the relativistic shifted-Maxwellian. Simple rejection methods are proposed to deal with the spatial part of the Lorentz transformation. We validate the algorithms by test problems, followed by discussions.

II. STATIONARY RELATIVISTIC MAXWELLIAN
We consider relativistic Maxwell distributions (Jüttner-Synge distribution 6,14 ) in the following form, where u = γv is the spatial components of the 4-velocity, v is the velocity, γ = [1 − (v/c) 2 ] −1/2 is the Lorentz factor, m is the rest mass, c is the light speed, T is the temperature, and K 2 (x) is the modified Bessel function of the second kind. The normalization constant is set such that the number density is N ≡ f (u)d 3 u. Hereafter we set m = 1 and c = 1 for simplicity. We use uppercases for fluid quantities and lowercases for particle properties throughout the paper.
To generate u, we consider the spherical transformation (u x , u y , u z ) = (u sin θ cos ϕ, u sin θ sin ϕ, u cos θ).
Then Equation 1 yields In the special case of N = 1, one can read this equation as a probability function with respect to u.
We generate u whose distribution follows Equation 2 by either the inverse transform method (Sec. II A) or the Sobol method (Sec. II B). We will describe these methods in the next subsections.
After we obtain u, we generate u on a spherical surface |u| = u in the momentum space.
Using uniform random variables X 1 (0 < X 1 ≤ 1) and X 2 (0 < X 2 ≤ 1), we set u in the following way, Then we obtain a relativistic Maxwellian which follows Equation 1.

A. Inverse transform method
We consider the cumulative distribution function F (u) with a practical upper bound In the nonrelativistic limit of T 1, u max = 5v th is sufficient, where v th = √ 2T is the thermal speed. In the relativistically hot case of T 1, Equation 2 behaves like ∝ exp(−u/T )u 2 for u 1. This decays slower than the nonrelativistic limit of ∝ exp[−(v/v th ) 2 ]v 2 , and so we increase the upper bound to u max = 20T . We usually prepare a numerical table of F (u) with 2000 or more grid points. Using a uniform random variable X 3 , we compute by referring and interpolating the table.

B. Sobol method
Let us consider the gamma distribution. Its probability function P (x) is given by where a and b are free parameters and Gamma(x) is the Gamma function. The gamma distribution with an integer parameter a can be generated by multiple random variables Sobol 12 noticed that the right hand side of Equation 2 is similar to the third-order Gamma distributions, For a certain T , we initialize u by using three random variables (X 4 . . . X 6 ),

PHYSICS OF PLASMAS
Comparing the exponential parts in Equations 2 and 8, we obtain a relativistic Maxwellian by the rejection method. By using another random variable X 7 , we accept the particle if Then we obtain u which is distributed by Equation 2. Using Equation 9, this criteria can be modified to This leads to a simple form of the Sobol's criterion, 10,11 where η = −T ln X 4 X 5 X 6 X 7 . Note that η and u share the same variables X 4 , X 5 , and X 6 .
Make sure to avoid zero in X 4 . . . X 7 , because ln 0 is undefined. Once Equation 12 (or Eq. 10) is met, we continue to the next step of the spherical scattering (Eq. 3).

Comparing the normalization factors in Equation 2 with N = 1 and Equation 8, we
obtain the overall efficiency of the rejection method as a function of T , 10,11 (13) Figure 1 shows the acceptance efficiency of the Sobol method, as a function of T . The efficiency quickly decreases for T → 0, while it approaches to 1 for T → ∞.

A. Lorentz transformation
Next we discuss general properties of the Lorentz transformation for particle distributions.
We consider the transformation between two frames, S and S . We assume that particles are stationary in the reference frame S, and then we switch to a moving frame S at the 4velocity (Γ, −Γβ, 0, 0). Without losing generality, we consider the transformation in the +x direction. In S , we observe the particle distribution, boosted by the 4-velocity (Γ, +Γβ, 0, 0).
We denote the observed properties in S by the prime sign ( ).
As the total particle number is conserved, we recognize Here, d 3 x = dx dy dz is the spatial volume element. Using the time element dt in the same frame, we consider the 4-dimensional volume element of dt and d 3 x that is moving at the 4-vector of u. The 4-dimensional position (t, x, y, z) follows the Lorentz transformation, and so the 4-volume element vector (dt, dx, dy, dz) also follows the Lorentz transformation. Since the Jacobian of the Lorentz transformation matrix Λ is 1, the 4-volume dt d 3 x is conserved, Since we deal with the u-moving volume, the time element dt is related to the canonical time element dτ in the following way, dt = γdτ . We similarly see dt = γ dτ . Therefore we obtain This also indicates the length contraction for the volume. The transformation is slightly Without losing generality, one can consider the Lorentz transformation by (Γ, −Γβ, 0, 0) in the +x direction: Under the condition of γdγ = u x du x , we obtain From Equations. 15 and 19, we obtain We obtain a relativistic shifted-Maxwellian by simply translating Equation 20, Since we know nice algorithms (Sec. II), we would like to initialize the particle momentum u in the S frame, and then translate it to the S frame by the Lorentz transformation, u → u . This procedure contains the momentum-space transformation (Eq. 19). However, it does not take care of the spatial part of the transformation, d 3 x → d 3 x (Eq. 15). Using the S-frame quantities, the distribution in S appears to the observer in the following way, We recognize a volume transform factor (γ /γ), because the element volume in S is not identical to the element volume in S . This issue is also addressed by Melzani et al. 9 .
One can also interpret that the number density is reciprocal to the volume size ∝ (d 3 x) −1 (See also Eq. 15). Since both spacial and momentum transformation (Equations 15 and 19) depends on u, the factor differs from particle to particle. This may sound tricky, but the above formula describes what the observer looks at. We obtain very different results without the volume transformation, as will be shown in Section IV.

B. Volume transform methods
Here, we describe simple methods to deal with the volume transform factor (Eq. 24).
It is impossible to deal with this by adjusting the cell size in PIC simulation, because the transformation differs from particle to particle. One can also change the weight of particles.
However, we prefer not to do so, because the ratio of the heaviest particle to the lightest particle could be very large.
We propose to adjust the particle number by a rejection method. Using a random variable X 8 (0 < X 8 ≤ 1), we accept the particle if the following condition is met, The left hand side ranges from 0 to 1. If the condition is not met, then we re-initialize the particle momentum. of the acceptance efficiency is 50%, If S is not the fluid rest frame, E[v x ] = 0 and so the efficiency may vary.
We can further improve the efficiency in a special case of a symmetric distribution. When f (u x ) = f (−u x ), we multiply the acceptance factor by 2, The (1/Γ) factor is absorbed in the total particle number. We take advantage of the fact that the second term of the right hand side is odd function of u x (or v x ). When βv x is negative, the acceptance factor ranges between 0 < (1 − |βv x |) ≤ 1. We reject the particles at the probability of |βv x |. On the other hand, when βv x is positive, the factor ranges between 1 ≤ (1 + βv x ) < 2. We accept all particles. In addition, we interpret that we need to add another set of particles at the probability of |βv x |. If f (u x ) = f (−u x ), the number of the rejected particles and the number of the particles to be added are equal. We simply reverse the sign of u x of the rejected particles, and then we add them to the positive-βv x side. This logic is schematically illustrated in the bottom of Figure 2.
We summarize the algorithm in the following way. If the following condition is met for a random variable X 9 , then we change u x → −u x , before computing u x . Here we combine the two conditions of −βv x < 0 and −βv x > X 9 to one condition (Eq. 32). The acceptance efficiency is 100%. We call it the flipping method (Eq. 32) to distinguish it from the rejection method (Eq. 29).

IV. TEST PROBLEMS
In order to validate the algorithms, we carry out several test problems. We initialize 10 6 particles in all cases. The black squares in Figure 1 show the acceptance efficiency of the Sobol method, as a function of T . They are in excellent agreement with the expected curve (Eq. 13).
We then compute relativistic shifted-Maxwellian by using the Sobol method and the flipping method (Eq. 32). We set T = 1 and we boost the particles by the bulk Lorentz factor Γ = (1, 1.1, 10) in the +u x direction. Figure 3 compares numerical results and analytic distributions in the moving frame S , integrated over u y and u z . All distributions are normalized by f d 3 u = ΓN . The following analytic solution is obtained by using a cylindrical transformation (u x , u y , u z ) = (u x , u ⊥ cos φ, u ⊥ sin φ) in Equation 22 .
The numerical results are in excellent agreement with the analytic solutions. The stationary Maxwellian looks OK. As Γ increases, the distribution is stretched in the +u x direction.

From Equation 33
, we see for u x → ∞. Therefore, the slope on the boosted side becomes extremely flat. For Γ = 1.1, the numerical results on the right side (u x ≈ 20) look slightly noisier than on the other side (u x ≈ −8). This is probably a unfair comparison, because the right slope has more gridpoints than the left slope in the low-density range. For Γ = 10, the distribution is highly stretched in +u x . Outside the figure, it still remains f (u x )/ΓN ≈ 4 × 10 −3 at u x = 100.
It will be challenging to initialize such a distribution by a direct rejection method in the S frame, because we have to extend the parameter domain 2Γ times longer in +u x . This gives us another motivation to initialize particles in S and then boost it to the S frame.
Next, we compute several fluid quantities in the moving frame S . After initializing the particles, we compute the flow vector N µ and the stress-energy tensor T µν . Then we evaluate the average velocity N x /N 0 = β and the average energy flux T 0x /N 0 = Γβ(E + P )/N . The former is a direct indicator of the bulk motion, and the latter, the energy flux, plays a decisive role to the system evolution. The results are presented in Table I. We change two key parameters, the bulk Lorentz factor Γ = (1.1, 10, 10 2 ) and the relativistic temperature T = (0.1, 1, 10). In the T = 0.1 case, we use the inverse transform method (Sec. II A), because the efficiency of the Sobol method falls to ≈ 0.001. We also test the T = 10 case without the volume transformation. This incorrect case is denoted by the asterisk sign ( * ). In Table I, the first rows show the computed results. The second rows show the relative error to analytic solutions. As can be seen, the results appear to be accurate, except for the rightmost columns.
Without the volume transformation, we see that the average bulk speed is inaccurate in Table I. This is crucial to initialize a relativistic current sheet 4,5 , in which relativistically hot populations carry the electric current. The energy flux is significantly distorted, too. The average energy flux without the volume transformation is Since (E + P )/E → 4/3 for T 1, we lose 25% of the energy flux, regardless of the bulk speed β. We can similarly evaluate the average energy density without the volume transformation. It deviates from the right value by a factor of [1 + Γ 2 −1 Γ 2 P E ] −1 , and therefore the error approaches 25% for Γ 1 and T 1.

V. DISCUSSION AND SUMMARY
We first reviewed two algorithms to initialize the stationary relativistic Maxwellian. In addition to the simple inverse transform method, we have formally reviewed the Sobol algorithm. In our experience, the inverse transform method is faster than the Sobol method, because it only requires 3 random variables. We don't see any problems, as long as we prepare 10 3 -10 4 grids in the table. The algorithm can deal with any spherically-symmetric distributions. On the other hand, the Sobol method has a strong mathematical background.
It is very simple, and so we can easily avoid a bug. The method is certainly slower than the inverse transform method, because it uses 6 random variables. However, this will not be a big deal, because we use these algorithms for initialization. The only problem is that the Sobol method becomes extremely inefficient for the nonrelativistic limit of T 1. In such a case, we simply switch from the Sobol method to the inverse transform method or the Box-Muller method. Another promising option is the log-concave rejection method, described in Section II and Appendix in Swisdak 13 . The algorithm uses 4 random variables, its acceptance efficiency is ≈90%, and it is nearly insensitive to T .
After initializing the stationary Maxwellian, we apply the volume transformation before boosting the particle momentum. We have proposed the two algorithms, the rejection method (Eq. 29) and the flipping method (Eq. 32). They require one more random variable. The flipping method is our first choice. Since it accepts all particles, the overall efficiency is the same as the base algorithm for the stationary one. As a representative case, the Sobol method with the flipping method are summarized in the pseudocode in Table II. We emphasize that our volume transform methods are generic. The flipping method can be combined with power-law, waterbag, or any other distributions, as long as it is symmetric in u x in the S frame. Even when the distribution is non-symmetric, we can divert to the rejection method (Eq. 29). The acceptance efficiency is typically 50%, but it works in any cases. Swisdak 13 used the log-concave rejection method twice for the shifted Maxwellian.
According to his article, the overall efficiency is ≈ 80% insensitive to T . This is a very good result. However, his algorithm is specialized for Maxwellians or possibly other exponentialtype distributions. In contrast, our simple methods can deal with any kind of distributions.
Using the test problems, we have demonstrated that the combinations of the base methods and the flipping method excellently work. Without the volume transformation, we recognize significant errors up to 25% in the average energy flux. This is because the volume transform repeat generate X 1 , X 2 , X 3 , X 4 , uniform on (0, 1] u ← −T ln X 1 X 2 X 3 η ← −T ln X 1 X 2 X 3 X 4 until η 2 − u 2 > 1.
In summary, we have described numerical algorithms to load relativistic Maxwellians in particle simulations. The inverse transform method and the Sobol method are useful to load the stationary Maxwellian. For shifted Maxwellian, the rejection method (Eq. 29) and the flipping method (Eq. 32) take care of the spatial part of the Lorentz transformation. These methods are simple and physically-transparent. They can be combined with arbitrary base algorithms. We hope that these algorithms are useful in relativistic kinetic simulations in high-energy astrophysics.

ACKNOWLEDGMENTS
The author acknowledges M. Oka and A. Taktakishvili for their assistance to find out