Non-ideal Particle Distributions from Kinetic Freeze Out Models

In fluid dynamical models the freeze out of particles across a three dimensional space-time hypersurface is discussed. The calculation of final momentum distribution of emitted particles is described for freeze out surfaces, with both space-like and time-like normals, taking into account conservation laws across the freeze out discontinuity.


Introduction
The freeze out of particle distributions is an essential part of continuum or fluid dynamical reaction models. From the point of view of observable consequences this is one of the most essential parts of the model. On the other hand this step is not based on fluid dynamical principles and governed by a large variaty of ad hoc assumptions. The freeze out can be considered as a discontinuity across a hypersurface in space-time.
The general theory of discontinuities in relativistic flow was not worked out for a long time, and the 1948 work of A. Taub 1 discussed discontinuities across propagating hypersurfaces only (which have a space-like normal vector, dσ µ dσ µ = −1). Events happening on a propagating, (2 dimensional) surface belong to this category.
Another type of change in a continuum is an overall sudden change in a finite volume. This is represented by a hypersurface with a time-like normal, dσ µ dσ µ = +1, called confusingly both space-like and time-like surface in the literature. In 1987 Taubs approach was generalized to both types of surfaces, 2 making it possible to take into account conservation laws exactly across any surface of discontinuity in relativistic flow. This approach also eliminates the imaginary particle currents arising from the equation of the Rayleigh line. When the EoS is different on the two sides of the freeze out front these conservation laws yield changing temperature, density, flow velocity across the front.
In fact the freeze out surface is an idealization of a layer of finite thickness where the frozen out particles are formed, and the interactions in the matter become gradually negligible. The dynamics of this layer can be described in different kinetic models or four-volume emission models. 3 The zero thickness limit of such a layer is the idealized freeze out surface.
The invariant number of conserved particles (world lines) crossing a surface element, dσ µ , is dN = N µ dσ µ , and the total number of all the particles crossing the FO hyper-surface, S, is N = S N µ dσ µ . This total number, N , and the total energy and momentum are of course the same at both sides of the freeze out surface. If we insert the kinetic definition of N µ into the integral over the freeze out hypersurface, S, we obtain the Cooper-Frye formula: 4 where f F O (x, p; T, n, u ν ) is the post FO phase space distribution of frozen-out particles which is not known from the fluid dynamical model. Problems usually arise from the bad choice of this distribution. First of all, to evaluate measurables we have to use the correct parameters of the matter after the FO discontinuity! If we know the pre freeze out baryon current and energy-momentum tensor, N µ 0 and T µν 0 , we can calculate locally, across a surface element of normal vector dσ µ the post freeze out quantities, N µ and T µν , from the relations 1,2 : In numerical calculations the local freeze out surface can be determined most accurately via self-consistent iteration. 7,9 This fixes the parameters of our post FO momentum distribution, f F O (x, p; T, n, u ν ).
For example we can illustrate the effect of conservation laws for a situation where the frozen out matter is massless, baryonfree Bose gas. Then, the conservation laws across the freeze-out surface with timelike normal vector dσ µ are [T µν dσ ν ] = 0 , In the most general (three dimensional) case there are four parameters to be determined from the conservation laws: The final, post FO temperature, T , and three components of the velocity, u. The energymomentum tensor on the pre freeze-out side, and the normal to the surface are given. The post freeze-out energymomentum tensor is of the form T µν = (e + p)u µ u ν − pg µν , where the energy density, pressure, and temperature are connected by the EoS: e = σ SB T 4 = 3p, where σ SB is the Stefan-Boltzmann constant. Then T µν = (e+p)u µ u ν −pg µν , can be written as a vector equation: where Taking the normal projection of (2) and the norm of the four-velocity, u µ , the solution for the four quantities we are looking for will bei given by: a. Idealized freeze out across propagating discontinuities. One can go a step further in the study of freeze out process. We usually assume that the pre freeze out momentum distribution as well as the post freeze out distribution are both local thermal equilibrium distributions boosted by the local collective flow velocity on the actual side of the freeze out hypersurface, although, the post freeze out distribution need not be a thermal distribution!
The case of freeze out across a hypersurface with space-like normal shows this clearly because p µ is time-like and dσ µ is space-like, thus p µ dσ µ can be both positive and negative. I.e., p µ may point now both in the post-and pre-FO directions. Thus the integrand in the above integral (1) may change sign in the integration domain, and this indicates that part of the distribution contributes to a current going back, into the front while another part is coming out of the front. On the pre freeze out side p µ is unrestricted and p µ dσ µ may really have both signs, because we may assume that the freeze out front has a certain thickness 8 , and due to internal rescatterings inside this front a current is fed back to the pre freeze out side to maintain the thermal equilibrium there.
On the post freeze out side, however, the distribution, f * F O (x, p; dσ µ ) must vanish for momentum four-vectors, p µ , which point back into the pre FO direction, i.e. do not satisfy the condition, p µ dσ µ > 0. 6,7 Thus, this distribution cannot be a Jüttner-or other ideal gas distribution. * Nevertheless, the above conservation laws, have to be satisfied, even if the post freeze out distribution is not a local thermal distribution. Since, the kinetic definitions of the energy-momentum tensor and conserved current(s) are reliably applicable, the conservation laws across a small element of the freeze out front take the form: * Note, that the contravariant normal when becomes space-like, dσ µ , should point into the pre-FO direction to satisfy the condition, p µ dσµ > 0, while the covariant normal, dσµ, always points into the post-FO direction! Thus, the direction of the contravariant normal dσ µ , in the Cooper-Frye formula goes continuously over from pointing to the pre-FO direction to pointing to the post-FO direction while the covariant normal of the FO surface stays directed always in the post-FO direction when it goes continuously over from time-like to space-like.
Here, the matter is characterized by T µν 0 and N µ 0 on the pre freeze out side of the front. The construction of the post freeze out distribution, f * F O , is a problem in the case of freeze out fronts with space-like normal.
For cut Jüttner distribution the conserved currents were evaluated in ref. 3 . Thus, if we know the 5 parameters of the pre FO flow and the local freeze out surface from kinetic considerations, then assuming that the post FO distribution, f * F O (p, x), is a cut Jüttner distribution, we can completely determine the parameters of the post FO matter from the conservation laws (4,5). Although, this way we would formally satisfy the conservation laws and we would eliminate the particle current pointing back to the pre FO matter, the strange shape of the cut Jüttner distribution makes it difficult to accept it as a physical post FO momentum distribution.

Freeze out distribution from kinetic theory
Following the ideas introduced in ref. 3 we can calculate the kinetic freeze out distribution based on four-volume emission models. The proposed model, on the other hand, requires extended numerical calculation, so here we intend to study some overly simplified models, which might give us some hints about the expected shape of post freeze out distributions.
The freeze-out will turn out to be an exponential process, and after about three mean free pathes the amount of interacting matter reduces to 5 per cent. Thus, the sharp freeze out layer turns out to be an over-idealization of kinetic freeze out in heavy ion reactions, while it is applicable on more macroscopic scales like in astrophysics. † Let us first demonstrate the kinetic model for a drastically oversimplified situation of a plane FO surface. Let us assume an infinitely long tube with its left half (x < 0) filled with nuclear mater and in the right vacuum is maintained. We can remove the dividing wall at t = 0, and then the matter will expand into the vacuum. By continuously removing particles at the right end of the tube and supplying particles on the left end, we can establish a stationary flow in the tube, where the particles will gradually freeze out in an exponential rarefaction wave propagating to the left. We can move with this front, so that we describe it from the reference frame of the front (RFF).
In this frame, we have a stationary supply of equilibrated matter from the left, and a stationary rarefaction front on the right, x > 0. We can describe the freeze out kinetics on the r.h.s. of the tube assuming that we have two components of our momentum distribution, f f ree (x, p) and f int (x, p). However, we assume only that at x = 0 f f ree vanishes exactly and f int is an ideal Jüttner distribution (supplied by the inflow of equilibrated matter), while f int gradually disappears and f f ree gradually builds up as x tends to infinity. We do not assume a priory that f int (x, p) is an ideal Jüttner distribution for all x, so we will have different FO results depending on the assumed FO mechanism. Let us take first the most simple kinetic model describing the evolution of such a system. Starting from a fully equilibrated Jüttner distribution the two components of the momentum distribution develop according to the coupled differential equations: Here the interacting component, f int , will deviate from the Jüttner shape and the solution will take the form This solution is depleted in the forward p-direction, particularly along the x-axis. Inserting it into the second differential equation above, leads to the freeze out solution: At x −→ ∞ this distribution will tend to the cut Jüttner distribution introduced in the previous section. (see Figs.  1, 2, and 3.) The remainder of the original Jüttner distribution survives as f int , even if x −→ ∞. In this model the particle density does not change with x, barely particles moving faster than the freeze out front (i.e. p µ dσ µ > 0) are transferred gradually from component f int to component f f ree . This is a highly unrealistic model, indicating that rescattering and re-thermalization should be taken into account in f int . This would allow particle transfer from the "negative momentum part" (i.e. p µ dσ µ < 0) of f int to f f ree , which is not possible otherwise. . At large distances from the initial point of the freeze out process, x −→ ∞ (C), the distribution becomes a cut Jüttner distribution, but less than half of the distribution is cut off! Note that the boosted Jüttner distribution became Lorentz elongated and asymmetric (see Fig. 2.10 of ref. 5 .) The earlier stages of the freeze out, here also, are characterized by asymmetric distributions, elongated in the freeze out direction, x.

Freeze out distribution with rescattering
The assumption that the interacting part of the distribution remains the distorted (after some drain) Jüttner distribution, is of course highly unrealistic. Rescattering within this component will lead to re-thermalization and re-equilibration of this component. Thus the re-equilibration and the drain terms are in competition and they mutually determine the evolution of the component, f int .
To include the collision terms explicitly into the transport equations, (6) leads to a combined set of integrodifferential equations. We can, however, take advantage of the relaxation time approximation to simplify the description of the dynamics.
Then the two components of the momentum distribution develop according to the coupled differential equations: Here, the interacting component of the momentum distribution shows the tendency to approach an equilibrated, Jüttner type, distribution with a relaxation length coefficient, λ ′ ≈ λ. Of course due to the energy, momentum and conserved particle drain, this distribution, f eq (x, p) is not the same as the initial Jüttner distribution, but its parameters, n eq (x), T eq (x) and u µ eq (x), change as required by the conservation laws. . At large distances from the initial point of the freeze out process, x −→ ∞ (C), the distribution becomes a cut Jüttner distribution, but more than half of the distribution is cut off! The earlier stages of the freeze out, here also are characterized by asymmetric distributions, but these are not elongated in the freeze out direction, x.

Conservation Laws.
In this case the change of the conserved quantities caused by the particle transfer from component int to component f ree can be obtained in terms of the distribution functions as: If we do not have collision or relaxation terms in our transport equation then the conservation laws are trivially satisfied. If, however, collision or relaxation terms are present these contribute, to the change of T µν and N µ , and this should be considered in the modified distribution function f int (x, p). b. Immediate re-thermalization limit. As a first approximation to the solution of eqs. (9) let us assume that λ ′ −→ 0, i.e. we have immediate re-thermalization after every step dx. Thus the drain is always happening from a component of shape f eq (x, p), with parameters,n(x), T (x) and u µ RF G (x), and we can assume that f int = f eq is of spherical Jüttner form at any x including both positive and negative momentum parts. Above and henceforth the notation is similar to the one in 3 :ñ = 8πT 3 e µ/T (2πh) −3 , a = m T , so thatn(µ, T ) =ña 2 K 2 (a)/2 is the invariant scalar density of the symmetric massless Jüttner i.e. K n (z, z) = K n (z). In this case the change of conserved quantities due to particle drain or transfer can be evaluated for an infinitesimal dx. We assume that the 3-flow is normal to the freeze out surface, and for simplicity we assume v > 0. In this case the change of the conserved particle currents in the RFF is given by K 1 (a, b))+ and for the change of the energy -momentum tensor in the RFF we obtain that and dT zz i = dT yy i . Note that in RFF the flow velocity of the re-thermalized component is The new parameters of distribution f int , after moving to the right by dx can be obtained from dN µ i and dT µν i . The conserved particle density of the re-thermalized spherical Jüttner distribution after a step dx iŝ where the expressions are invariant scalars. After straightforward calculation the differential equation describing the change of the proper particle density is: Although this covariant equation is valid in any frame, we can calculate it in the RFF, where the values of dN µ i s were given above. Note again, that the particle drain from f int (x), described by dN µ i is constrained to the "positive part" in the momentum space, but after re-thermalization we attribute this to the change in the complete spherical Jüttner distribution, f int (x + dx). Thus, in order to conserve momentum, we have to obtain a decreased Eckart flow velocity after the infinitesimal particle drain.
For the re-thermalized interacting component Eckart's flow velocity is the velocity of the RFG, which changes with x, so we can actually denote this frame as RFG(x). For the spherical Jüttner distribution the Landau and Eckart flow velocities are the same, u µ i,E,RF G (x) = u µ i,L,RF G (x) = u µ i,RF G (x). Thus we can evaluate the flow velocity u µ i,RF G (x + dx): which leads to the following covariant expression where ∆ µν , is a projector to the plane orthogonal to u µ i,RF G (x). This equation is valid in any reference frame, nevertheless we know the four-vectors on the r.h.s. in the RFF explicitly. Then the new flow velocity of the matter evaluated according to Eckart's definition is . At large distances from the initial point of the freeze out process, x −→ ∞ (C), the distribution becomes a cut Jüttner distribution, which is not centralized in py and less than half of the distribution is cut off! The earlier stages of the freeze out, here also are characterized by distributions asymmetric both in the directions px and py, and these are also elongated in the direction of the freeze out flow velocity, u µ RF G .
To get the temperature and the change of Landau's flow velocity, we have to analyze the change of the energy momentum tensor. Before the particle drain the energy -momentum tensor at x in the RFG is diagonal, . Adding the drain terms, dT µν i (x), to this arising from the freeze out while we move to the right by dx, yields T µν i (x + dx) which will not be diagonal in the RFG(x) and the pressure part will not be isotropic. We can Lorentz transform this to another frame which diagonalizes T µν i (x + dx). This means to find the Landau flow velocity of the new system, u µ i,L,RF G (x + dx) in the original RFG(x). After a straightforward diagonalization, a somewhat tricky algebra and neglecting second and higher order terms we arrive at the covariant expression ‡ Although, for the spherical Jüttner distribution the Landau and Eckart flow velocities are the same, the change of this flow velocity when calculated from the baryon current and from the energy current are different This is a clear consequence of the asymmetry caused by the freeze out process as we pointed out already at the discussion of the properties of the cut Jüttner distribution. Unfortunately, this also illustrates the weakness of our assumption on the complete re-thermalization to a spherical Jüttner distribution, because we cannot choose the correct velocity change: If we choose du µ E as the new velocity of the (spherical Jüttner distribution) f int (x + dx), then we violate the momentum conservation in our model, on the other hand if we choose du µ L , then we violate the baryon current conservation! Thus a spherical (or even elliptic) distribution cannot be fitted to the freeze out drain, and we would have to use an ansatz, which has (in addition) an asymmetry in the x direction (i.e., an egg shape), for the distribution f int .
Being aware of this weakness of the model, we nevertheless, maintain the assumption of spherical Jüttner shape for f int for the sake of simplicity. We can choose the flow velocity change then according to the physical problem. For example for the freeze out of baryon free plasma this problem does not occur, and we have to choose du µ L . The last item is to determine the change of the temperature parameter of f int . From the relation e ≡ u µ T µν u ν we readily obtain the expression for the change of energy density and from the relation between the energy density and the temperature (see Chapter 3 in ref. 5 ), we can obtain the new temperature at x + dx. Fixing these parameters we fully determined the spherical Jüttner approximation for f int . With this ansatz the pressure asymmetry and pressure balance cannot be realized, thus our model will be only a rather approximate description of the freeze out process. Nevertheless, we can draw some preliminary conclusions about the development of the kinetic distribution during freeze-out.

Conclusions
We turned to the problem of estimating the freeze out distribution. Obviously the real freeze out distribution depends strongly on the details of the freeze out (and hadronization) dynamics. In heavy ion reactions, the curvature of the freeze out surface and the conditions varying in time do affect the freeze out distribution, nevertheless, as a first step, we assumed that the process is stationary and the curvature of the front is negligible. These approximations are extreme, but still enable us to draw some preliminary conclusions. ‡ Let the energy-momentum tensor of a system be T µν . The energy and momentum flow is characterized by the Landau flow velocity, a unit four vector, uµ. We are looking for a relationship between the infinitesimal change of the flow velocity duµ and the corresponding shift in the energy-momentum tensor dT µν . We introduce the projector, ∆ µν = g µν − u µ u ν , with the properties 5 ∆ µν uν = 0 and duµ = ∆ ν µ duν since uµdu µ = 0 . The Landau flow velocity is parallel to the flow of the momentum. Thus uµ = Const. × T ν µ uν , therefore ∆ρµT µν uν = 0 . We differentiate the above equation and take into consideration the identities e ≡ uµT µν uν and ∆ µ ρ T ρσ ∆ ν σ = −P ∆ µν , where e and P are the energy density and pressure of the dissipationless, fully equilibrated fluid. Then using the properties of ∆ µν we get duρ(e + P ) + uρduµT µν uν = ∆ρµdT µν uν . Since the flow velocity and the momentum flow are parallel the second term on the l.h.s vanishes. Thus the equation describing the change of Landau's flow velocity becomes duρ = ∆ρµdT µν uν/(e + P ) .
Following the lines and ideas presented in ref. 3 , the first simple kinetic freeze out model reproduces the cut Jüttner distribution as the limiting distribution, f f ree , after complete freeze out at large distances. However, the model at the same time leads to unrealistic consequences, namely that the interacting part of the distribution, f int also survives fully, as the other part of the Jüttner distribution. Thus having both components at the end in this model, the physical freeze out is actually not realized. This turns out to be a consequence of the fact that the effect of rescattering and thermalization in the interacting part of the distribution was ignored.
In an improved but still rather approximate kinetic freeze out model which takes rescatterings into account, the interacting component is assumed to be instantly re-thermalized taking a spherical Jüttner shape at each time step with changing parameters. The model leads to a set of coupled differential equations (10,11,12,13). Equations (11) and (12) can be used in some combined form, or one of them can be selected which fits the physical situation the best. Then the three parameters of the interacting component, f int , can be obtained in each time step analytically (considering K n (x, y) an analytic function). Now the density of the interacting component will gradually decrease and disappear according to eq. (10), the flow velocity will also decrease in both cases, (11) or (12), because only forward going particles freeze out, and the energy density will decrease also according to eq. (13). Thus, the initial contribution to f f ree at small x will resemble the distribution shown in Fig. 2A, then as x increases and the velocity decreases it will become to similar to Fig. 1B, while at the final stages it will approach Fig. 3C. As a consequence the integrated distribution will not resemble a cut Jüttner distribution.
Thus the arising post freeze out distribution, f f ree will be a superposition of cut Jüttner type of components, from a series of gradually slowing down Jüttner distributions. This will lead to a comet shaped final momentum distribution, with a more dominant leading head and a tail. In these rough models a large fraction (∼ 95%) of the matter is frozen out by x = 3λ, thus the distribution f f ree at this distance can be considered as a first estimate of the post freeze out distribution. One should also keep in mind that the models presented here do not have realistic behavior in the limit x −→ ∞, due to their one dimensional character. Nevertheless, this improved model with rescattering enables complete freeze out (unlike the simpler model in sect. 2 where only the originally forward moving particles freeze out even at large distances).
In case of rapid hadronization of QGP and simultaneous freeze out, the idealization of a freeze out hypersurface may be justified, however, an accurate determination of the post freeze out hadron momentum distribution would require a nontrivial dynamical calculation.