Generation of random deviates for relativistic quantum-statistical distributions

We provide an algorithm for generation of momenta (or energies) of relativistic particles according to the relativistic Bose-Einstein or Fermi-Dirac distributions. The algorithm uses rejection method with effectively selected comparison function so that the acceptance rate of the generated values is always better than 0.9. It might find its use in Monte-Carlo generators of particles from reactions in high-energy physics.


Motivation
In projects related to multiparticle production in hadronic or nuclear collisions it is often demanded to generate a large number of particles with momenta distributed according to relativistic Bose-Einstein of Fermi-Dirac distribution. Here one has to take into account the total energy (i.e. including the mass) when evaluating the exponent of the distributions where m is the mass of the particles, p = | p|, µ is the chemical potential and T is temperature. Parameter q assumes the value of 1 for fermions and −1 for bosons. For the generation of invariant momentum distributions one also needs this distribution multiplied with the energy As the distribution is spherically symmetric, the angles are trivially integrated and we are left with the distributions for the size of the momentum vector. In order to make it suitable for a general procedure, it is expressed with the help of dimensionless variable Thus we get or for the other distribution In these functions we have suppressed the constant pre-factors which also contain dimensions.
For the Monte Carlo generation we shall proceed with the dimensionless distributions without the dimensionfull pre-factors.
A similar algorithm for the generation of relativistic Maxwellian distribution has been reported in [1,2]. In our work we properly account for quantum statistics and allow for non-zero chemical potential, which can influence the momentum distribution. We describe the procedure for the distribution (1.5) with the energy factor. The procedure for the distribution (1.4) can easily be derived along the same steps as we shall proceed.

The algorithm
We demonstrate in the Apendix that the distribution (1.5) is log-concave for large enough x, i.e. its logarithm is a concave function. For such a distribution there always exists an exponential that is everywhere above the demanded distribution. One can generate random deviates according to the exponential and use the rejection method [2].
In order to achieve smallest possible rejection rate we use piecewise analytic comparison function, as indicated in Fig. 1. The three pieces are determined so that • for x ≤ x − the comparison function is linear; • for x − < x ≤ x + the comparison function is constant and equal to the value of the distribution at the mode; • for x > x + the comparison function is exponential. The joint points x − and x + are chosen so that the comparison function is always continuous.
For the determination of the comparison function we thus need to determine the five points indicated on the horizontal axis.
x m The mode of S(x). This is obtained easily by differentiating and we get Unfortunately this expression cannot be solved analytically and numerical methods must be invoked.
Subsequently, the value of S(x) at the mode can be determined x l The point left from the mode in which the linear comparison function touches the distribution. It is found from the condition for the derivative of the distribution 3) The slope of the liner comparison function is then x − The point in which the linear part of the comparison function and its constant part meet. It can be determined as The knowledge of x − also allows to express x u The point in which the exponential part of the comparison function touches the distribution. Above the mode the distribution is log-concave. Therefore, x u can be chosen anywhere above x m . However, we checked that, the acceptance rate is optimised with x u chosen so that the distribution there drops to 1/e of its maximum value. We get the value by solving the equation for the logarithms of the distribution lnS(x u ) = lnS m −1, which leads to Again, this equation must be solved numerically.
Once x u is determined, we can determine the slope parameter of the exponential comparison function. It is given by the logarithm of the distribution. Thus x + Finally, this is the point in which the constant part of the comparison function and its exponential part join. It is determined from a simple equation Once we have x + , we also know the exponential part of the comparison function which reads S ′ (x) = S m e −λ(x−x + ) . (2.10) Thus we can formulate the comparison function In order to use this comparison function as probability density (after normalisation) for random variate generation we need the values The inverse if the integral of S ′ (x) is For the rejection step we need the probabilities to accept the generated value of x. They are given as P (x) = S(x)/S ′ (x). In the three intervals they read Now we have collected all expressions needed to build up the algorithm. Because of the need to numerically solve a few equations the procedure may be lengthy if one needs to generate just one random value. However, if many values must be generated for the same temperature, chemical potential and particle mass, then all parameters can be calculated first and then used repeatedly. Thus the first part of the algorithm is the calculation of the parameters: 8. For later convenience calculate also the values of y − , y + , and y ∞ from Eqs. (2.12). This is the common part of the preparation. Then, in order to generate a value, follow these steps: 9. Generate uniform random deviate y from the interval [0,y ∞ ].
11. Accept the value of x with the probability given by Eqs. (2.14). If the value is not accepted, return to step 9.

Illustration of results
We have tested this algorithm in wide range of parameters A and M. Particle masses were chosen both smaller than temperature so that large momenta are available and also much larger than temperature so that the momenta are practically non-relativistic. Chemical potentials up to the value of particle mass for bosons, i.e. the point of condensation, were tested, as well ( Figure 2). In all cases the acceptance rate was around 90%. This shows that the comparison function is very well adapted to the present problem.

Conclusions
The presented algorithm has been successfully implemented in an upgrade of the Monte Carlo event generator DRAGON [3], which serves for the generation of hadrons produced in high energy nuclear collisions. It is, however, general and can serve in any other application where relativistic momenta must be generated from quantum-statistical distributions.

A Log-concave distribution
In this Appendix we demonstrate that the distribution S(x) according to Eq. (1.5) is indeed log-concave on the interval above the mode. Therefore, an exponential function which touches S(x) from above in one point will never be smaller than S(x). The calculation is straightforward. We take the second derivative of lnS(x). For fermions (q = 1), this leads to .
Note that s f (x) ≤ 1 for any x. Therefore, (1−s f (x)) ≥ 0 and all terms in the bracket in Eq. (A.1) are non-negative. In summary, we see that for fermions and thus the distribution is log-concave everywhere. The case of bosons is slightly more involved. Again, we take the second derivative Note the change of the sign in (1+s b (x)) and in front of the last term. Due to this, for bosons the second derivative may become positive in some cases. We want to demonstrate that such pathological intervals are always below the mode of S(x). For x→∞ the terms s b (x) go to 0 exponentially (we chose the letter s for "small"), and one can inspect that lim x→∞ d 2 lnS(x) dx 2 = 0 and the value of the limit is being approached from below. Thus the second derivative is either negative everywhere or there is a point x = x c where it crosses the horizontal axis and stays negative for x > x c . It is enough to show that x c < x m . It turns out that the second derivative becomes positive only if the particles are light with A ≤ 3 and M is close to 1, which is quite an extreme case. (Recall that for bosons M must be smaller than 1 and M = 1 corresponds to the condensation point where the distribution does no longer apply.) In such a case, for small values of x one obtains s b ≫1 and the last term in Eq. (A.3) prevails. We have scanned the whole relevant parameter region of A and M and checked that always x c < x m (Figure 3). In all other cases the function S(x) is log-concave everywhere.
We conclude that it is safe to use the exponential comparison function in the interval (x m ,∞).