Freeze out in hydrodynamical models

We study the effects of strict conservation laws and the problem of negative contributions to final momentum distribution during the freeze out through 3-dimensional hypersurfaces with space-like normal. We study some suggested solutions for this problem, and demonstrate it on one example.


Introduction
Fluid dynamical models, especially their simpler versions are very popular in heavy ion physics, because they connect directly collective macroscopic matter properties, like the Equation of State (EoS) or transport properties, to measurables.
Particles which leave the system and reach the detectors, can be taken into account via source (drain) terms in the 4-dimensional space-time based on kinetic considerations, or in a more simplified way via freeze out (FO) or final break-up schemes, where the frozen out particles are formed on a 3-dimensional hypersurface in space-time. This information is then used as input to compute measurables such as two-particle correlation, transverse-, longitudinal-, radial-, and cylindrical-flow, transverse momentum and transverse mass spectra, etc.
In this paper we concentrate on freeze out. A basic standard assumption in this case is that freeze out happens across a hypersurface as already mentioned, so it can be pictured as a discontinuity where the kinetic properties of the matter, such as energy density and momentum distribution change suddenly. The hypersurface is an idealization of a layer of finite thickness (of the order of a mean free path or collision time) where the frozen-out particles are formed and the interactions in the matter become negligible. The dynamics of this layer is described in different kinetic models such as Monte Carlo models 1,2 or four-volume emission models [3][4][5][6][7] . In fact, the zero thickness limit of such a layer is an over-idealization of kinetic freeze out in heavy ion reactions, while it is applicable on more macroscopic scales like in astrophysics * . * On the other hand, if kinetic freeze out coincides with a rapid phase transition, like in the case of rapid deconfinement transition of supercooled quark-gluon plasma, the sharp freeze out hypersurface idealization may still be applicable even for heavy ion reactions. It is, however, beyond the scope of this work to study the freeze out dynamics and kinetics in this latter case.
Two types of hypersurfaces are distinguished: those with a space-like normal vector, dσ µ dσ µ = −(dσ) 2 (e.g. events happening on a propagating 2-dimensional surface) and those with a time-like normal vector dσ µ dσ µ = (dσ) 2 (a common example of which is an overall sudden change in a finite volume).
Once the freeze out surface is determined, one can compute measurables. Landau when drafting his hydrodynamical model 8 , just evaluated the flow velocity distribution at freeze out, and this distribution served as a basis for all observables. This approach was used in early fluid dynamical simulations of heavy ion collisions also [9][10][11] . This procedure was improved to add thermal velocities to the flow velocities at freeze out, by Milekhin 12,13 and later by Cooper and Frye 14 . This method is widely used, however it rises at least three problems 15 .
First, in some cases before the 90's, the possible existence of discontinuities across hypersurfaces with time-like normal vectors was not taken into account or considered unphysical †19-22 . This point was studied recently 23 so we do not discuss it further.
Second, since the kinetic properties of the matter are different on the two sides of the front, the explicit evaluation of conservation laws across the freeze out surface should be taken into account which is not always easy to implement. In some (simple) cases [24][25][26] , these conservation laws are enforced and discussed. For example in 24 , it was pointed out that the freeze out momentum distribution for hypersurfaces with time-like normal may become locally anisotropic. In 27 a solution for post FO massless, baryonfree Bose gas was presented. Here we remind the procedure that should be followed in section 2.
The third problem is a conceptual problem arising in the Cooper-Frye freeze out description when we apply it to a hypersurface with space-like normal: it is the problem of negative contributions (see section 2). This is the main subject of this paper. This problem appears in all freeze out calculations up to now we are aware of, and to our knowledge it was not satisfactorily discussed yet in the literature. It was recognized by some of those who applied the Cooper-Frye freeze out description before 25,26,28 . A possible partial solution was presented in part 2 of ref. 26 for noninteracting massless particles, in 1+1 dimension using the post FO cut Jüttner ansatz. In 27 it was shown that in an oversimplified kinetic freeze out model one can obtain the cut Jüttner distribution as post FO distribution. In section 4 we complete and generalize the results of 26 and present an example for the solution of the freeze out problem. In section 5 we suggest improvements that go beyond the cut Jüttner ansatz.

Conservation laws across idealized freeze out discontinuities
In the zero width limit of the freeze out domain (freeze out surface), the energy -momentum tensor changes discontinuously across this surface. Consequently, the four-vector of the flow velocity may also change 17,29,30 . These changes should be discussed in terms of the conservation laws.
The invariant number of conserved particles (world lines) crossing a surface element, dσ µ , is and the total number of all the particles crossing the FO hyper-surface, S, is If we insert the kinetic definition of N µ † Taub 16 discussed discontinuities across propagating hypersurfaces, which have a space-like normal vector. If one applies Taub's formalism from 1948, to freeze out surfaces with time-like normal vectors, one gets a usual Taub adiabat but the equation of the Rayleigh line yields imaginary values for the particle current across the front. Thus, these hypersurfaces were thought unphysical. However more recently, Taub's approach has been generalized to these hypersurfaces 17 (see also 18 ) while eliminating the imaginary particle currents arising from the equation of the Rayleigh line. Thus, it is possible to take into account conservation laws exactly across any surface of discontinuity with relativistic flow. into eq. (1) we obtain the Cooper-Frye formula 14 : where f F O (x, p) is the post FO phase space distribution of frozen-out particles which is not known from the fluid dynamical model. The problem is to choose its form correctly. Usually one assumes that the pre-FO momentum distribution as well as the post FO distribution are both local thermal equilibrium distributions, with the same temperature, boosted by the local collective flow velocity on the actual side of the freeze out surface, although the post FO distribution need not be a thermal distribution. Parametrizing the post FO distribution as thermal, f F O (x, p; T, n, u µ ), where T does not necessarily coincide with the pre FO temperature, and knowing the pre-FO baryon current and energy-momentum tensor, N µ 0 and T µν 0 , we can calculate the post freeze out quantities N µ and T µν from the relations 16,17 [N µ dσ µ ] = 0 and [T µν dσ µ ] = 0, across a surface element ‡ of normal vector dσ µ . Here [A] ≡ A − A 0 . The post FO distribution is not a thermal equilibrium distribution, so temperature does not exist, nevertheless, the conservation laws fix the parameters, e.g. T , n, u µ , of our momentum distribution, f F O (x, p; T, n, u ν ).
To obtain a physically realizable result, in addition we have to check the condition for entropy increase: where, for both equilibrium and nonequilibrium FO distributions 18 This condition is not necessary to obtain a solution of the freeze out problem, but it should always be checked to exclude non-physical solutions. We have to note at this stage, that the post FO distribution must not be an equilibrium (or stationary) solution of the Boltzmann Transport Equation, and consequently on the post FO side the energy flow-, baryon flow-and entropy flow-velocities may all be different.
We can now remind briefly what the problem of negative contributions to the Cooper-Frye formula is, and a possible way out. For a FO surface with time-like normal, both p µ and dσ µ are time-like vectors, thus, p µ dσ µ > 0, and the integrand in the integral (4) is always positive. For a FO surface with space-like normal, p µ is time-like and dσ µ is space-like, thus, p µ dσ µ can be both positive and negative. (Note that p µ may point now both in the postand pre-FO directions.) Thus, the integrand in the integral (4) 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-FO side p µ is unrestricted and p µ dσ µ may have both signs, because we are supposing that pre-FO phase is in the thermal equilibrium. However, in the zero width limit of the FO front, it is difficult to understand such a situation. What happens actually is that internal rescatterings occur inside the finite FO domain and feed particles back to the pre-FO side to maintain the thermal equilibrium there. On the post-FO side, however, we do not allow rescattering and back scattering any more. If a particle has passed the freeze out domain it cannot scatter back. In other words, the post-FO distribution should have the form 25,26 , where Θ(x) is the step function. Consequently, this distribution cannot be an ideal gas distribution. (On the pre-FO side, the distribution may or not be ideal). The conservation laws across a small element of the freeze out front with space-like normal take the form: 3. The allowed momentum region for space-like FO Let us assume that the FO process happens in the positive x direction, in other words we go in the positive x-direction from the pre FO domain to the post FO domain. The FO hyper-surface has a space-like normal (dσ µ dσ µ = −(dσ) 2 ), so that dσ µ is orthogonal to the hyper-surface (i.e. to the time-like tangent vector of the surface t µ , t µ dσ µ = 0) and points into the post FO (positive x-) direction (while dσ µ points in the pre FO direction).
Depending on the reference frame, the space-like FO front can propagate both in the positive or negative xdirection, or it can be Lorentz transformed into its own rest frame, the Rest Frame of the Front (RFF), where dσ µ = (0, 1, 0, 0)dσ. In other reference frames t 0 is frame dependent and may be both positive and negative.
In order to use the parameter v in the following discussion, we have to select a given reference frame and fix the value in that frame. Let us choose the frame comoving with the peak of the post FO invariant momentum distribution. Following ref. 26 let us denote this frame as the Rest Frame of the Gas (RFG). Note, however, that, contrary to what its name seems to suggest, this frame is not the Local Rest frame of the post FO matter, since the post FO distribution is not spherically symmetric in momentum space! Thus, for the following discussion we define v in the RFG frame as , and this means that in the RFF the peak four-velocity has the form u µ RF G = γ(1, v, 0, 0)| RF F . This velocity v can also be both positive and negative (!) in RFF, i.e. the peak velocity may point to the post FO direction (as we would expect), and also in the pre FO direction in RFF which seems to be confusing. This is, however, not a problem in itself because both flow velocities (Eckart and Landau) are always positive on the post FO side. As we will see later, in special cases it is possible that we obtain negative post FO peak velocity for positive pre FO flow velocities. This indicates we have to discuss the importance of the cut by Θ(p µ dσ µ ), otherwise one might be tempted to believe that this cut is not affecting the post FO distribution too much, and the correct treatment causes only a few percent cut which is negligible. We show in the following that this is not the case.
The p µ dσ µ > 0 requirement in the RFG frame means that only momenta with component contribute to the post FO momentum distribution. This means that the boundaries of the allowed domain in the [p x , p ⊥ ]-plane are hyperbolas in the post FO RFG and the domains of the positive x-side of the corresponding hyperboloids in the 3-dimensional momentum space may contribute to the FO distribution (Fig. 1). In the massless limit the hyperboloids become cones around the x-axis and centered at the origin.  a. Dominant case. For v = 0 the hyperboloid becomes a plane boundary at p x ≥ 0 in RFG, i.e. the boundary cuts the Jüttner distribution in the middle (because in the RFG the peak of the distribution is centered at p = 0). If v > 0, the peak velocity points to the post FO direction and the cut-hyperboloid is fully in the negative p x region, thus for large v values only a smaller fraction is excluded from the full Jüttner distribution. Most of the particles described by a full Jüttner distribution freeze out in this case. We will consider this as the dominant FO case.
On the other hand, as we can see in Fig. 1, for the v < 0 case the peak flow points back into the pre FO direction, and only a small fraction of a full Jüttner distribution will freeze out. The cut eliminates the major part of the particles. This situation leads to less frozen out particles, but it yields an unusual post FO distribution.

Conserved currents for cut Jüttner distribution
We now study the particular case where f F O is a Jüttner (or relativistic Boltzmann 32,18 ) distribution and so f * F O is a cut Jüttner distribution: where µ is the chemical potential related to the invariant scalar densityn, of the non-cut Jüttner distribution as µ = T ln n(2πh) 3 / 4πm 2 T K 2 m T , and u µ RF G ,n and T are parameters of the distribution f * F O originating from the full Jüttner distribution. These are not the flow velocity, proper density and temperature of the cut Jüttner distribution! The cut Jüttner distribution is not a thermal equilibrium distribution, e.g. it does not have a temperature at all.
This distribution for massless particles was considered in part 2 of ref. 26 (also in 25 ). The cut selects particles with momenta p µ dσ µ > 0.
We can evaluate the baryon four-current, N µ , by inserting the cut Jüttner distribution into the definition, eq. (3), and we get a time directed-, N 0 , as well as a spatial component, N x , (where x is the direction of the spatial component of the FO-normal, dσ µ , in RFG). In ref. 26 the spatial component, N x , of the four current was not evaluated, so seeing only the zeroth component, N 0 , the unsuspecting reader might have believed falsly, that RFG (Rest Frame of the Gas) is the Local Rest frame of the gas. Performing the calculation, in the post FO RFG frame the baryon current reads as where j = sign(v),ñ = 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 Jüttner gas i.e. K n (z, z) = K n (z). Just as in case of the non-cut distributions the cut Jüttner distribution yields few modified Bessel functions in the expression of the four currents, while the relativistic Fermi and relativistic Bose distributions lead to a series of these functions. When evaluating the limits we used the relation K n (a, b) a=b −→ K n (a) Note that the proper density of the cut Jüttner distribution, n, is reduced compared to the proper density of the complete spherical Jüttner distribution,n. The energy momentum tensor in the post FO RFG is where B = (1 + b + b 2 /2 + b 3 /6)e −b and T zz = T yy . This energy-momentum tensor may then be Lorentz transformed into the Landau Local Rest (LLR) frame of the post FO matter, which moves with u µ L in the RFG, or into the Rest Frame of the Freeze out Front (RFF) where dσ µ = (0, 1, 0, 0)dσ. Alternatively both can be transformed to the frame where we want to evaluate the conservation laws, eq. (5), and the parameters of the post FO, cut Jüttner distribution can be determined so, that it satisfies the conservation laws. In the massless limit the energy momentum tensor in the RFG is: and T zz = T yy = (T 00 − T xx )/2.
In addition we have to check the entropy condition. In the RFG frame the entropy current reads Note that in the m = 0 limit the vectors S µ and N µ are parallel to each other. This is explained by Figure 1, which shows that in the RFG the cut in the m = 0 limit becomes a central cone, and since the distribution is centrally symmetric in this frame, the integrals will be proportional to each other.

A. Solubility of the freeze out problem
The situation is non-trivial and we have to take into account the possible directions of the flow and of dσ µ . Note: we must not assume that the flow is parallel to the freeze out direction. Let us start on the pre FO side labeled by '0'. Here in the LR frame u µ 0 = (1, 0, 0, 0)| pre−LR and we can choose the x-direction in this frame to point into the FO direction, so that dσ µ = γ 0 (v 0 , 1, 0, 0)dσ| pre−LR . We assume that we know the FO hypersurface, i.e. we know v 0 . Then, in this frame the conservation laws have three nonvanishing components yielding three known parameters N µ 0 dσ µ , T 0µ 0 dσ µ , and T xµ 0 dσ µ . To find the solution we need these values in the RFG frame. However, the 3-dimensional direction of the xaxis will not change because the front is assumed to be isotropic in its own [y, z]-plane. Thus, in the RFG the peak flow parameter is u µ RF G = (1, 0, 0, 0)| RF G , and the normal of the FO front is n µ = γ(v, 1, 0, 0). Note that v = v 0 ! Furthermore, let us recall that the parameter v determines the post FO peak flow parameter in RFF, Consequently the conservation laws (9,10) yield three nonvanishing equations in the RFG frame, While the first equation is an invariant scalar, the remaining two are components of a 4-vector, so they should be transformed into the same reference frame, i.e., to RFG. Since we evaluated the quantities based on the cut-Jüttner distribution in the RFG, we also need the pre FO quantities in the RFG. These can be determined by using the standard fluid dynamical form of T µν as seen from the RFG. From this frame the pre FO flow velocity is given by the difference of the pre-and post-FO flow velocities: In the general case the solution can be obtained numerically. In the m = 0 limit the solution is simpler and gives an interesting insight into the problem. The continuity equation leads to the equation where Q −1 = Q −1 1 (µ, T ) =ñ(µ, T )/(4n 0 γ 0 v 0 ) , which leads to a 3rd order equation and can be solved for v alnalytically. The energy equation, [T 0µ dσ µ ] = 0, leads to the same equation but with another coefficient Then, dividing the equation [T 0µ dσ µ ] = 0, by the equation [T xµ dσ µ ] = 0, yields another 3rd order equation for v: where R 0 = e 0 v 0 /p 0 . This equation can be solved analytically and yields one physical root, v = 2 − 3/R 0 (and two unphysical ones v = −1). Inserting v then into eq. (18), we can obtain the resulting chemical potential, µ, also. The possibility of this simple analytic solution is a consequence of the fact that in the m = 0 limit the cut of the Jüttner distribution is made along central cones in the RFG, which then divide the energy and the baryon charge exactly in the same proportions.
As an illustration we studied the freeze out of quark-gluon plasma (QGP) to cut Jüttner gas, in the massless limit. The pre FO side QGP is described by the most simple bag-model EoS (eqs. (5.28-33) in ref. 18 ), thus local equilibrium is assumed and all pre FO parameters are assumed to be known including the baryon, energy-momentum and entropy currents.
On the post FO side these currents were evaluated earlier in this section, and the equations arising from the conservation laws, (5), were solved as presented above. Figure 2 indicates the change of flow velocity during freeze out. Physical solution exists only for positive initial velocities, v 0 ≥ 0. The velocity parameter (!) of the post FO cut Jüttner distribution varies from −1 to +1, but the post FO Eckart flow velocity is of course always positive in RFF. Thus, the post FO baryon current is also positive in RFF (this is obvious since we do not allow any particle to cross the front backwards), and consequently, the pre FO current and v 0 should also be positive because of the continuity equation. For small initial velocities, v 0 −→ 0, the post FO velocities approach zero also, but for moderate velocities, deduced recently from experiments, v = 0.3 − 0.7, the difference between the post and pre FO flow velocities may be essential.
In order to show the effect of these modifications compared to the original Cooper-Frye treatment (where the increase of the flow velocity is ignored) we can consider case (a) in Figure 2. The cut Jüttner distribution always leads to an exponential p t spectrum, but according the new modified treatment starting from v 0 = 0.2c the post FO flow velocity increases to v f low = 0.4c, while the post FO parameter velocity (which determines the p t spectrum) increases to v = 0.6c. This corresponds to an increase of the slope parameter, T slope , by 60%! This is due to the large latent heat arising from the large value of the bag constant taken in case (a). In case (b) the the same effect is present but it is weaker. This change of the flow velocity is a basic feature of the correct freeze-out treatment, and it is a consequence of the conservation laws and not of the positivity requirement of p µ dσ µ in spacelike FO. Thus, the flow velocity change occurs both in spacelike and in timelike freeze-out. This effect can cause for example the conversion of latent heat to collective kinetic energy and not to heat if the freeze-out coincides with an exotherm phase transition 29 . Figure 3 shows that the baryon density, (15), decreases in the freeze out process. This is connected to the fact that the post FO flow velocities are above the pre FO ones, as shown in Figure 2.
We should mention that the post FO temperature parameter of the cut Jüttner distribution becomes rather high, about an order of magnitude higher than the pre FO temperature. However, we have to recall that the term temperature is not applicable for a non-equilibrium distribution, therefore this result has no physical significance, it just illustrates the parametrization of the distribution of the assumed cut Jüttner shape.
Finally we have to check the entropy condition for these solutions. As we know 29,20,33 QGP can freeze out to hadronic matter with entropy production only if the QGP is supercooled or considerably supercooled. This remains valid for the cut Jüttner assumption as post FO distribution also. With most parametrizations only low temperature QGP is able to freeze out. For the cut Jüttner gas we cannot speak of a critical temperature, because this gas is not in equilibrium and consequently cannot be in phase equilibrium either. Still this distribution can be attributed an entropy current by its kinetic definition, and the entropy condition can be checked (Fig. 4).  In reality the entropy condition is not so stringent as Figure 4 indicates. In this illustrative study the post FO EoS had relatively few degrees of freedom to accomodate the high entropy content of QGP. By including many post FO mesons and other hadronic degrees of freedom in our post FO EoS, the entropy condition can be satisfied in a much wider range of parameters.

Freeze out distribution from kinetic theory
We have seen that taking the cut Jüttner distribution as an ansatz for the post FO distribution, we can solve the freeze out problem formally. Although we can satisfy all requirements, the obtained parameter values make it questionable whether the cut Jüttner ansatz is an adequate assumption. The shape of the distribution with the sharp cuts is also a rather unphysical feature of the distribution.
To obtain a more realistic, and physically better applicable FO distributions, we should evaluate the distribution in more physical, microscopic non-equilibrium models. Kinetic theory is a straightforward candidate for this task.
A first very simplified attempt to solve the freeze out problem dynamically in one dimensional kinetic model 27 returned the cut Jüttner distribution also, but only in highly unrealistic situations: only when the model yielded incomplete freeze out. Thus, further work is needed to find physically realistic post freeze out distributions in kinetic models or in other dynamical microscopic models.

Conclusions
The importance of taking into account conservation laws in the description of the freeze out process is pointed out. For freeze out across hypersurfaces with space-like normals the approach suggested by Bugaev, 26 assuming cut Jüttner distribution as post freeze out distribution is worked out, and the freeze out problem was solved as an example for QGP freezing out into a cut Jüttner gas. This calculation indicates that results including the Cooper-Frye freeze out procedure should be reconsidered and new emphasis should be given to the precise evaluation of the post freeze out particle distributions.
The deviation from the earlier Cooper-Fry approach (where changes of flow velocity, density and temperature were ignored) is apparent if the pre FO matter has large energy content in form of compressional energy, latent heat, or in any other way, which is not present in the post FO, noninteracting matter. As this post FO matter is not necessarily in thermal equilibrium, we cannot consider it as a thermal phase with equilibrium thermodynamical parameters. Thus, this idealized approach assuming a FO surface is always assuming a discontinuity irrespective of what was the phase of the pre FO matter. Nevertheless, this treatment leads to the strongest modifications in cases when a first order phase transition with large latent heat is coupled to the freeze out process.
Here we have considered an idealized transition as a discontinuity across a hypersurface. In as much as the flow across the surface is stationary our results are valid irrespective of the surface thickness, because we used only conservation laws. On the other hand in heavy ion reactions the flow across the surface can be considered stationary only if it is 1-2 fm wide. With purely kinetic freeze out this is not a very realistic assumption. 1 On the other hand rapid hadronization from supercooled QGP may satisfy the required conditions and the sharp surface approximation is then realistic. 34