fl) CrossMark

yll fcliqk for updates

Loading relativistic Maxwell distributions in particle simulations

Seiji Zenitania)

National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

(Received 30 January 2015; accepted 15 April 2015; published online 29 April 2015)

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 «50% for generic cases and 100% for symmetric distributions. They can be combined with arbitrary base algorithms. © 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1063/L4919383]

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 instabilities15 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 rela-tivistic 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 Jutter-Synge distribution function.6,14 In order to load it, perhaps the Sobol12 algorithm is the most popular, at least in Monte-Carlo simulation community. The algorithm was formally proposed by Sobol12 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 Jutter-Synge distribution have not been applied to the relativistic shifted-Maxwellian. Several alternative algorithms have been proposed. Swisdak13 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

a)Electronic mail: seiji.zenitani@nao.ac.jp

1070-664X/2015/22(4)/042116/6

describe two base algorithms to load stationary relativistic Maxwellian, the inverse transform method and the Sobol method. 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 (Juttner-Synge distribution6,14) in the following form:

f (u)d3

4nm2cTK2 (mc2 / T)

ymc2\ 3 exp I--) d u, (1)

where u = yv is the spatial components of the 4-velocity, v is the velocity, y =[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 K2 (x) is the modified Bessel function of the second kind. The normalization constant is set such that the number density is N = jf (u)d3u. 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 (ux, uy, uz) = (u sin hcos u, u sin h sin u, u cos h). Then, Eq. (1) yields

f(u)du =

TK2(1/T)

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 Eq. (2) by either the inverse transform method or the Sobol method. We will describe these methods in Subsections IIA and IIB.

After we obtain u, we generate u on a spherical surface |u| = u in the momentum space. Using uniform random variables X1 (0 < X1 < 1) and X2(0 < X2 < 1), we set u in the following way:

22, 042116-1

© Author(s) 20151

= u (2Xi - 1)

= 2uv/X1(1 - X1) cos(2pX2) = 2uv/X1 (1 - X1 ) sin(2pX2).

Then, we obtain a relativistic Maxwellian which follows Eq. (1).

A. Inverse transform method

We consider the cumulative distribution function F(w) with a practical upper bound wmax,

F(u) =

f (u)du

f (u)du

f(u)du

f (u)du

In the nonrelativistic limit of T ^ 1, umax = 5vh is sufficient, where vth = is the thermal speed. In the relativis-tically hot case of T S1, Eq. (2) behaves like / exp(—u/T)u2 for u ^ 1. This decays slower than the nonrelativistic limit of / exp[—(u/uth)2]v2, and so we increase the upper bound to umax = 20T. We usually prepare a numerical table of F(u) with 2000 or more grid points. Using a uniform random variable X3, we compute

u = F-1(X3), (5)

by referring and interpolating the table.

B. Sobol method

Let us consider the gamma distribution. Its probability function P(x) is given by

P(x; a, b) =

ba Gamma(a)

xa-1e-x=b (x > 0),

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 X;'s (0 < X;- < 1) in the following way:7

Sobol12 noticed that the right hand side of Eq. (2) is sim-

ilar to the third-order Gamma distributions,

P(u; 3, T)= U^jexpf- ^V

For a certain T, we initialize u by using three random variables (X4.X6),

lnX4 - lnX5 - lnX6 = -lnX4X5X6.

Comparing the exponential parts in Eqs. (2) and (8), we obtain a relativistic Maxwellian by the rejection method. By using another random variable X7, we accept the particle if

Then, we obtain u which is distributed by Eq. (2). Using Eq. (9), this criteria can be modified to

\/1 + u2 < (u - TlnX7) = -TlnX4X5X6X7. (11)

This leads to a simple form of the Sobol's criterion,10,11

g2 - u2 > 1; (12)

where g = —TlnX4X^X6X7. Note that g and u share the same variables X4, X5, and X6. Make sure to avoid zero in X4.. X7, because ln0 is undefined. Once Eq. (12) (or Eq. (10)) is met, we continue to the next step of the spherical scattering (Eq. (3)).

Comparing the normalization factors in Eq. (2) with N = 1 and Eq. (8), we obtain the overall efficiency of the rejection method as a function of T,10,11

^ K2 (1 =T).

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 !i.

III. RELATIVISTIC SHIFTED-MAXWELLIAN 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 4-velocity (r, — Tb, 0,0). Without losing generality, we consider the transformation in the direction. In S', we observe the particle distribution, boosted by the 4-velocity (T, +T@, 0,0). We denote the observed properties in S0 by the prime sign (0).

As the total particle number is conserved, we recognize

f (x, u) d3xd3u = f V, u') d3x' d3u'.

Here, d3x = 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 d3x that is moving at

FIG. 1. Acceptance efficiency of the Sobol algorithm as a function of the temperature t. The black squares show numerical results in Sec. IV.

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 A is 1, the 4-volume dtd3x is conserved, i.e., dtd3x = dt'd3/. Since we deal with the u-moving volume, the time element dt is related to the canonical time element ds in the following way: dt = yds. We similarly see dt' = y'ds. Therefore, we obtain

-d3x = -' d3 X.

This also indicates the length contraction for the volume. The transformation is slightly different for d3 u, because u is constrained by uMuM = u2 — y2 = — 1. Without losing generality, one can consider the Lorentz transformation by (r, —rb, 0, 0) in the +x direction:

y = ry(1 + bvx), duX = r(dux + bdy), duy = duy, duZ = duz.

Under the condition of ydy = uxdux, we obtain d3u d3u'

y y' .

From Eqs. (15) and (19), we obtain d3xd3u This ensures

f (x, u)=f'(x', u').

d3x' d3".

We obtain a relativistic shifted-Maxwellian by simply translating Eq. (20),

f (и) = /'(«')

4pTK2(1/T) eX^ T

4pTK2 (1/T)

Г(-- ß<

(21) (22)

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, d3x ! d3x' (Eq. (15)). Using the S-frame quantities, the distribution in S' appears to the observer in the following way:

/'(и')d3U = f(и) - d3u

We recognize a volume transform factor (y'/y), 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

Г(1 + ßvx).

One can also interpret that the number density is reciprocal to the volume size / (d3x)_1 (see also Eq. (15)). Since both spacial and momentum transformation (Eqs. (15) and (19)) depend on u, the factor differs from part/cle to part/cle. 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 Sec. IV.

In this line, we briefly outline relativistic fluid properties. We assume isotropic Maxwellian distribution. From Eq. (19), the number flux 4-vector N1 yields

ft \ i d3" f (и)и1 —.

We see V'1 = (V,V'V') = V(r, Tb). Equation (25) ensures that V1 follows the Lorentz transformation, i.e., V'1 = A^V", where A is the Lorentz tensor.

Similarly, the stress-energy tensor T1^ yields,

f (u)wV

Clearly, it follows the Lorentz transformation, i.e., T'^ = A^Ab^"b. In this case,

•/00

■'0i

Г2(Е + P)-P, = Г2(£ + P)ß',

where E = J/(u)yd3w = V{[K3(1/T)/K2(1/T)] - T} is the internal energy density and P = J/(u)wx(wx/y)d3w = NT is the pressure in the rest frame.

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 particle-in-cell (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 (0 < < 1), we accept the particle if the following condition is met:

= ! (1 + ßVx) > X8.

The left hand side ranges from 0 to 1. If the condition is not met, then we re-initialize the particle momentum. The factor (1 /2T) can be absorbed in the normalization constant, because we usually know the value of 2TV before loading particles. The expected value E[x] of the acceptance efficiency is 50%,

2Г V С

-(1 + ßE[v*]) = 0.5.

If S is not the fluid rest frame, E[vx] = 0 and so the efficiency may vary.

Fluid rest frame (S)

Observer frame (S')

Acceptance factor

Recycle

Momentum part of the transformation

Y«-r(7 + /3ux) u'x^T(ux+P>y)

Spatial part of the transformation

FIG. 2. Lorentz transformation of a relativistic hot plasma distribution. The bottom panel illustrates the flipping method, which is responsible for the spatial part of the Lorentz transformation.

We can further improve the efficiency in a special case of a symmetric distribution. When / («x) = /(—«x), we multiply the acceptance factor by 2,

= (1 + ßvx).

The (1/r) 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 «x (or vx). When @vx is negative, the acceptance factor ranges between 0 < (1 — |bvx|) < 1. We reject the particles at the probability of |bvx|. On the other hand, when bvx is positive, the factor ranges between 1 < (1 + bvx) < 2. We accept all particles. In addition, we interpret that we need to add another set of particles at the probability of |bvx|. If /(«x) =/(—«x), the number of the rejected particles and the number of the particles to be added are equal. We simply reverse the sign of «x of the rejected particles, and then we add them to the positive-bvx 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 X9,

-ßüx > X9,

1e-01 1e-02 1e-03 1e-04 1e-05

1 \ ■ Г=1

\ \ 0 Г=1.1

T \ X Г=ю

', , ,1 1 Ф, , , », , . 1 . .r1..............

FIG. 3. Distribution functions /'(«^) of Lorentz-boosted Maxwellians as a function of «x. Numerical results are overplotted on the analytic curves (Eq. (33)). We set t = 1 for all cases.

then we change before computing Here, we

combine the two conditions of — bvx < 0 and — @vx > X9 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 106 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 r =(1,1.1,10) in the +wx direction. Figure 3 compares numerical results and analytic distributions in the moving frame S', integrated over and m^. All distributions are normalized by J/'d3w' = CN. The following analytic solution is obtained by using a cylindrical transformation (mX, My, mZ) = (X, m? cos /, m? sin /) in Eq. (22).

f (U )u'Ld/d.u'L

N Г у/1 + u'2 + T 2Г2К2(1/Г)

+ u'2 - ßux)

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'r direction.

From Eq.

-ßuX)/T)]

(33), we see /'K)z «x exp[—(r^ 1 + «'2 i «x exp[—(«x/2rT)] for «x !i. Therefore, the slope on the boosted side becomes extremely flat. For r = 1.1, the numerical results on the right side («x « 20) look slightly noisier than on the other side («x « —8). This is probably an unfair comparison, because the right slope has more gridpoints than the left slope in the low-density range.

TABLE I. Computed fluid quantities and relative errors.

r = 1.1 r = 10

r = 100

r = 1.1 r = 10

r = 100

T = 0.1 T = 1.0 T = 10 T = 10*

0.416532 0.416708 0.416566 0.288896

1.6 x 10—4 2.6 x 10—4 7.7 x 10—5 0.307

0.994989 0.994996 0.994958 0.975918

1.6 x 10—6 8.9 x 10—6 2.9 x 10—5 0.0192

0.999950 0.999950 0.999950 0.999658

9.3 x 10—10 4.0 x 10—8 1.3 x 10—7 2.9 x 10—4

T = 0.1 T = 1.0 T = 10 T = 10*

0.580438 2.00167 18.3798 13.7357

2.9 x 10—4 5.5 x 10—4 1.4 x 10—3 0.252

12.6063 43.4805 398.147 298.838

2.6 x 10—6 1.1 x 10—4 8.5 x 10—4 0.250

126.674 437.062 4001.75 3003.43

1.5 x 10—4 9.2 x 10—5 7.4 x 10—4 0.250

For r = 10, the distribution is highly stretched in +^1. Outside the figure, it still remains /'(mD/CN « 4 x 10—3 at = 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 2r times longer in +m1. 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'1 and the stress-energy tensor T'1^. Then we evaluate the average velocity N'x/N'0 = b and the average energy flux r°7N'° = rb(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 r = (1.1,10,102) 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 relativisti-cally hot populations carry the electric current. The energy flux is significantly distorted, too. The average energy flux without the volume transformation is

J / (w)w'0(w'7w'0)d3w_ 1 " = N

J/ (u)d3M

/ (w)r(bM0 + M*)d3M

Since (E+ P)/E ! 4/3 for T » 1, we lose 25% of the energy flux, regardless of the bulk speed b. We can similarly evaluate the average energy density without the volume transformation. It deviates from the right value by a factor of [1 + (p) ]_1, and therefore the error approaches 25% for r » 1 andT » 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 do not see any problems, as long as we prepare 103-104 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 Sec. 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 is 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 Mx 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. Swisdak13 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 exponential-type distributions. In contrast, our simple methods can deal with any kind of distributions.

TABLE II. Sobol algorithm with the flipping method. repeat

generate xj, x2, x3, x4, uniform on (0, 1]

m <--t ln xix2x3

g <--t ln xjx2x3x4

until g2 — m2 > 1.

generate x5, x6, x7, uniform on [0, 1] mx ^ m (2x5 - 1) my ^ 2^^x5(1 — x5)cos(2px6) mz ^ 2^^x5(1 — x5) sin(2px6) if ( — bvx > x7), mx <--mx

mx ^ r(mx + bp1 + m2) return mx, my,mz

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 factor (Eq. (24)) is no longer constant for C > 1 and T ^ 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 Sobol's original article and T. N. Kato for his insightful comments on the manuscript. This work was supported by Grant-in-Aid for Young Scientists (B) (Grant No. 25871054).

1E. P. Alves, T. Grismayer, S. F. Martins, F. Fiuza, R. A. Fonseca, and L. O. Silva, "Large-scale magnetic field generation via the kinetic Kelvin-Helmholtz instability in unmagnetized scenarios," Astrophys. J. 746, L14 (2012).

2G. E. P. Box and M. E. Muller, "A note on the generation of random normal deviates," Ann. Math. Stat. 29, 610 (1958).

3Y. A. Gallant, M. Hoshino, A. B. Langdon, J. Arons, and C. E. Claire, "Relativistic, perpendicular shocks in electron-positron plasmas," Astrophys. J. 391, 73 (1992).

4E. G. Harris, "On a plasma sheath separating regions of oppositely directed magnetic field," Nuovo Cimento 23, 115 (1962).

5F. C. Hoh, "Stability of sheet pinch," Phys. Fluids 9, 277 (1966).

6F. Juttner, "Das Maxwellsche Gesetz der Geschwindigkeitsverteilung in der Relativtheorie," Ann. Phys. 339, 856 (1911).

7W. J. Kennedy, Jr. and J. E. Gentle, Statistical Computing (Marcel Dekker Inc., 1980).

8L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields, 3rd ed. (Pergamon press, Oxford, 1971), Secs. 6 and 10.

9M. Melzani, C. Winisdoerffer, R. Walder, D. Folini, J. M. Favre, S. Krastanov, and P. Messmer, "Apar-T: Code, validation, and physical interpretation of particle-in-cell results," Astron. Astrophys. 558, A133 (2013).

10L. A. Pozdnyakov, I. M. Sobol, and R. A. Sunyaev, "Effect of multiple Compton scatterings on an X-ray emission spectrum—Calculations by the Monte Carlo method," Sov. Astron. 21, 708 (1977).

11L. A. Pozdnyakov, I. M. Sobol, and R. A. Sunyaev, "Comptonization and the shaping of X-ray source spectra—Monte Carlo calculations," Astrophys. Space Phys. Rev. 2, 189 (1983).

12I. M. Sobol, "On modeling certain distributions similar to gamma distribution," in Monte Carlo Methods in Computational Mathematics and Mathematical Physics (Novosibirsk, 1976), pp. 24-29 (in Russian).

13M. Swisdak, "The generation of random variates from a relativistic Maxwellian distribution," Phys. Plasmas 20, 062110 (2013).

14J. L. Synge, The Relativistic Gas (Interscience, New York, 1957).

15S. Zenitani and M. Hoshino, "Particle acceleration and magnetic dissipation in relativistic current sheet of pair plasmas," Astrophys. J. 670, 702 (2007).

Physics of Plasmas is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.