Covariant description of kinetic freeze out through a finite space-like layer

The problem of Freeze Out (FO) in relativistic heavy ion reactions is addressed. We develop and analyze an idealized one-dimensional model of FO in a finite layer, based on the covariant FO probability. The resulting post FO phase-space distributions are discussed for different FO probabilities and layer thicknesses.


I. INTRODUCTION
The hydrodynamical description of relativistic particle collisions was first discussed more than 50 years ago by Landau [1] and nowadays it is frequently used in different versions for simulations of heavy ion collisions. Such a simulation basically includes three main stages. The initial stage, the fluid-dynamical stage and the so-called Freeze Out (FO) stage when the hydrodynamical description breaks down. During this latter stage, the matter becomes dilute, cold, and non-interacting, the particles stream toward the detectors freely, their momentum distribution freezes out. Thus, the freeze out stage is essentially the last part of a collision process and the source of the observables.
The usual recipe is to assume the validity of hydrodynamical treatment up to a sharp FO hypersurface, e.g. when the temperature reaches a certain value, T F O . When we reach this hypersurface, all interactions cease and the distribution of particles can be calculated.
In such a treatment FO is a discontinuity where the properties of the matter change suddenly across some hypersurface in space-time. The general theory of discontinuities in relativistic flow was first discussed by Taub [2]. That description can only be applied to discontinuities across propagating hypersurfaces, which have spacelike, (dσ µ dσ µ = −1), normal vector. The discontinuities across hypersurfaces with time-like, (dσ µ dσ µ = 1), normal vector were considered unphysical. The remedy for this came only 40 years later in [3], generalizing Taub's approach for both time-like and space-like hypersurfaces. Consequently, it is possible to take into account conservation laws exactly across any surface of discontinuity in relativistic flow.
As it was shown recently in [4,5], the frequently used Cooper-Frye prescription [6] to calculate post FO particle spectra gives correct results only for discontinuities across time-like normal vectors. The problem of negative contributions in the Cooper-Frye formula was healed by a simple cut-off, Θ(p µ dσ µ ), proposed by Bugaev [4]. However, this formulation is still based on the existence of a sharp FO hypersurface, which is a strong idealization of a FO layer of finite thickness [7]. Thus, by assuming an immediate sharp FO process, the questions of final state interactions and the departure from local equilibrium are left unjustified.
The recent paper [8,9] formulates the freeze out problem in the framework of kinetic transport theory. The dynamical FO description has to be based on the Modified Boltzmann Transport Equation (MBTE), rather than on the commonly used Boltzmann Transport Equation (BTE). The MBTE abandons the local molecular chaos assumption and the requirement of smooth variation of the phase-space distribution, f (x, p), in space-time. This modification of BTE, makes it even more difficult to solve the FO problem from first principles. Therefore, it is very important to build phenomenological models, which can explain the basic features of the FO process.
The present paper aims to build such a simple phenomenological model. The kinetic approach presented, is applicable for FO in a layer of finite thickness with a space-like normal vector. It can be viewed as a continuation and generalization of [10,11,12]. The kinetic model for FO in time-like direction was discussed in a recent paper [13], however, the fully covariant model analysis and the treatment are presented in [14].
In present work we use stationary, one-dimensional FO models for the transparent presentation. Such models can be solved semianalytically, what allows us to trace the effects of different model components, assumptions and restrictions applied on the FO description. We do not aim to apply directly the results presented here to experimental heavy ion collision data, instead our purpose is to study qualitatively the basic features of the FO process. We want to demonstrate the applicability of the proposed covariant FO escape rate, and most im-portantly, to see the consequences of finishing FO in a finite layer. Up to now, two extreme ways of describing FO were used: (i) FO on an infinitely narrow hypersurface and (ii) infinitely long FO in volume emission type of model. To our knowledge this is the first attempt to, at least qualitatively, understand how FO in finite space-time domain can be simulated and what will be its outcome. In such stationary, one-dimensional models the expansion cannot be realistically included, therefore it is ignored.
In realistic simulations of high energy heavy-ion reactions the full 3D description of expanding and freezing out system should be included. This work is under initial development.

II. FREEZE OUT FROM KINETIC THEORY
A kinetic theory describes the time evolution of a single particle distribution function, f (x, p) = f (t, x, p 0 , p), in the 6D phase-space. To describe freeze out in a kinetic model, we split the distribution function into two parts [10,15,16]: The free component, f f , is the distribution of the frozen out particles, while f i , is the distribution of the interacting particles. Initially, we have only the interacting part, then as a consequence of FO dynamics, f i gradually disappears, while f f gradually builds up. In this paper we convert the description of the FO process from a sudden FO, i.e. on a sharp hypersurface, into a gradual FO, i.e. in some finite space-time domain. Freeze out is known to be a strongly directed process [17], where the particles are allowed to cross the FO layer only outwards, in the direction of the normal vector, dσ µ , of the FO hypersurface. Many dynamical processes happen in a way, where the phenomenon propagates into some direction, such as detonations, deflagrations, shocks, condensation waves, etc. Basically, this means that the gradients of the described quantity (the distribution function in our case) in all perpendicular directions can be neglected compared to the gradient in the given direction dσ µ , i.e. ▽f ≈ dσ µ ∂ µ f . In such a situation these can be effectively described as one-dimensional processes, and the space-time domain, where such a process takes place, can be viewed as a layer.
Therefore, we develop a one-dimensional model for the FO process in a layer of finite thickness, L. We assume that the boundaries of this layer are approximately parallel, and thus, the thickness of the layer does not vary much. This can be justified, for example, in the case when the system size is much larger than L. At the inside boundary of this layer there are only interacting particles, whereas at the outside boundary all particles are frozen out and no interacting particles remain. Note that the normal to the FO layer, dσ µ , can be both spacelike or time-like. The gradual FO model for the infinitely long onedimensional FO process was presented in recent works [10,11,12]. We are going to build a similar model, but now we make sure that FO is completely finished within a finite layer.

A. Freeze out in a finite layer
In kinetic theory the interaction between particles is due to collisions. A quantitative characterization of collisions is given by the mean free path (m.f.p.), giving the average distance between collisions. The m.f.p., λ mf p , is inversely proportional to the density, λ mf p ∼ 1/n(x). If we have a finite FO layer, the interacting particles inside this domain must have a finite m.f.p. During the FO process, as the density of the interacting particles decreases, they are entering into a collisionless regime, where their final m.f.p., tends to infinity, or at least, gets much larger than the system size L. The realistic FO process for nucleons in a heavy ion collisions happens within a finite space-time FO domain, which has a thickness of a few initial mean free paths [18]. Hence, one must realize that the FO process cannot be fully exploited by the means of the m.f.p. concept, since we have to describe a process where we have on average a few collisions per particle before freeze out. Therefore, this type of processes should be analyzed by having also another characteristic length scale different from the m.f.p. In our case it should be related to the thickness, L, of the FO layer.
Recent conjectures based on strong flow and relatively small dissipation find that the state where collective flow starts is strongly interacting and strongly correlated while the viscosity is not large [19]. This indicates a small m.f.p., in the interacting matter, while at the surface λ mf p → ∞. Several indications point out that in high energy heavy ion reactions freeze out and hadronization happens simultaneously from a supercooled plasma [20,21,22]. This could be modeled in a way that prehadron formation and clusterization starts gradually in the plasma, and this process is coupled to FO in a finite layer. The FO is finished when the temperature of the interacting phase drops under a critical value and all quarks cluster into hadrons, which no longer collide. This is the possible qualitative scenario with well defined finite thickness L of the FO layer. Now, let us recall the equations describing the evolution in the simple kinetic FO model [10,11,12]. Starting from a fully equilibrated Jüttner distribution, f J (p), i.e. f i (s = 0, p) = f J (p) and f f (s = 0, p) = 0, the two components of the momentum distribution develop in the direction of the freeze out, i.e. along dσ µ , according to the following differential equations: where W esc (s, p) is the escape rate governing the FO development and s = x µ dσ µ . Here x µ is a 4-vector having its origin 1 at the inner surface, S 1 , of the FO layer, see Fig. 1. In order to obtain the probability to escape, for a particle passing from 0 till s, P esc (s), we have to integrate the escape rate along a trajectory crossing the FO layer: The definition for the escape probability was previously given in [16], in terms of collision or scattering rates, where the FO process was lasting infinitely long. In our finite layer FO description the quantity that defines the escape probability is the escape rate. To have a complete physical FO finished at a finite distance/time, we require: P esc → 1, when s → L. In usual cascade models the probability of collision never becomes exactly zero, and correspondingly P esc never becomes exactly one, and the FO process lasts ad infinitum. This is due to the fact that the probability of collision is calculated based on the thermally averaged cross section, which does not vanish for thermal, e.g. Gaussian, momentum distributions. In reality the free or frozen where the cut-off factor, Θ(p s ) ≡ Θ(p µ dσ µ ), forbids the FO of particles with momenta not pointing outward [4]. This FO parameter, λ ′ (s, p), is not necessarily an average distance in space or duration in time between two subsequent collisions, like the m.f.p. The m.f.p., tends to infinity as the density decreases, while the FO just becomes faster in this limit. Actually, the FO scale behaves in the opposite way to the m.f.p. This can be seen, for example, in a simple, purely geometrical freeze out model, which takes into account the divergence of the flow in a 3D expansion [23]. Both this and the phase transition or clusterization effect described at the beginning of this section, lead to a finite FO layer L, even if the m.f.p., λ mf p ∼ = 1/nσ is still finite at the outer edge of this layer. We consider the thickness of the layer L to be the "proper" thickness of the FO layer, because it depends only on invariant scalar matter properties like cross section, proper density, velocity divergence, phase transition or clusterization rates. These should be evaluated in the Local Rest frame (LR) of the matter, and since the layer is finite, around the middle of this layer. The proper thickness is analogous to the proper time, i.e. time measured in the rest frame of the particle, hence the proper thickness is the thickness of the FO layer measured in the rest frame attached to the freeze out front, that is, the Rest Frame of the Front (RFF). Some of the parameters like the velocity divergence and the phase transition rate describe the dynamical changes in the layer, so these can determine the properties, e.g. the thickness, of the finite layer. However, calculating L from the above mentioned properties is beyond the scope of this paper, and L is treated as a parameter in the following.
Let us consider the Rest Frame of the Front, where the normal vector of the front points either in time, t, or in space, x, direction, introducing the following notations 2 . Indeed, if dσ µ is space-like the resulting equations can be transformed into a frame where the process is stationary (here dσ µ = (0, 1, 0, 0) and correspondingly s ≡ x), while in the case of a time-like normal vector the equations can be transformed into a frame where the process is uniform and time-dependent (dσ µ = (1, 0, 0, 0), s ≡ t). For sake of transparency and simplicity we will perform calculations only for FO in a finite layer with space-like normal vector in this paper, but many intermediate results can and will be obtained in Lorentz invariant way. Inside the FO layer particles are separated into still colliding or interacting and not-colliding or free particles. The probability not to collide with anything on the way out, should depend on the number of particles, which are in the way of a particle moving outward in the direction p/|p| across the FO layer of thickness L, see Fig. 1. If we follow a particle moving outward form the beginning, (x µ = 0), i.e. the inner surface of the FO layer, S 1 , to a position x µ , there is still a distance L − s cos θ p ahead of us, where θ p is the angle between the normal vector and p/|p|. As this remaining distance becomes smaller the probability to freeze out becomes larger, thus, we may assume that the escape rate is inversely proportional to some power, a, of this quantity [9,24].
Based on the above assumptions we write the escape rate as: where this newly introduced parameter, λ, is the initial, i.e. at S 1 , characteristic FO length of the interacting matter, λ = λ ′ (s = 0, cos θ p = 1). The power a is influencing the FO profile across the front. Indeed, calculating the escape probability, P esc , eq. (3), with the escape rate, given by eq. (5), we find , for a = 1, and for a = 1. Thus, we see that for different a-values we have different FO profiles: a = 1: power like FO, a > 1: fast, exponential like FO, a < 1: no complete FO within the finite layer, since P esc does not tend to 1 as s approaches L.
In papers [10,11,12] the authors were using a = 1, and were modeling FO in an infinite layer. In order to study the effects of FO within a finite space-time domain, we would like to compare the results of our calculations with those of earlier works, therefore we shall also take a = 1 in further calculations. It is easy to check that our escape rate, eq. (5), equals the earlier expression in L → ∞ limit. Thus, the model discussed in this paper is a generalization of the models for infinitely long FO, described in [10,11,12], and allows us to study FO in a layer of finite thickness. The angular factor, cos θ p , maximized the FO probability for those particles, which propagate in the direction closest to the normal of the FO layer. For the FO in timelike directions, studied in [13], the angular factor was 1. This factor, and correspondingly the escape rate, eq. (5), are not covariant. Furthermore, this earlier formulation does not take into account either that the escape rate of particles should be proportional to the particle velocity (the conventional non-relativistic limit of the collision rate contains the thermal average < σv >). Let us consider the simplest situation, when the Rest Frame of the Front is the same as Rest Frame of the Gas (RFG), where the flow velocity is u µ = (1, 0, 0, 0). If freeze out propagates in space-like direction, i.e. dσ µ = (0, 1, 0, 0), as shown in Fig. 1, then cos θ p = p x /|p|. Therefore, a straightforward generalization of the escape rate, based on the above arguments, is where the r.h.s. of this equation is an invariant scalar in covariant form. Now, we assume that this simple generalization is valid for any space-like or time-like FO direction, even when RFG and RFF are different [9,17,25]. Based on the above arguments, we can write the total escape rate from eq. (5) in a Lorentz invariant 3 form: which now opens room for general study of FO in relativistic flow in layers of any thickness. Former FO calculations in [10,11,12] were always performed in RFF. Aiming for semianalytical results and transparent presentation, as well as in order to compare our results with former calculations, we will also study the system evolution in RFF, but now this is only our preference. In principle calculations can be performed in any reference frame. In more realistic many dimensional models, which will take into account the system expansion simultaneous with the gradual FO, it will be probably more adequate to work in RFG or in Lab frame, and our invariant escape rate, eq. (7), can be directly used as a basic FO ingredient of such models.

B. The Lorentz invariant escape rate
In this section let us study this new angular factor, in more detail. We will take the p-dependent part of the escape rate, eq. (7), and denote it as: In RFG, where the flow velocity of the matter is u µ = (1, 0, 0, 0) RF G by definition, W (p), is given as and it is smoothly changing as the direction of the normal vector changes in RFG. This will be discussed in more detail in the rest of this section.
In the following, we will take different typical points of the FO hypersurface, A, B, C, D, E, F, see Fig. 2. At these points, the normal vectors of the hypersurface, dσ µ = (h, i, j, k) RF G , are given below. To calculate the normal vector for different cases shown in Fig. 2, we simply make use of the Lorentz transformation. The normal vector of the time-like part of the FO hypersurface may be defined as the local t ′ -axis, while the normal vector of the space-like part may be defined as the local x ′ -axis. As dσ µ is normalized to unity, its components may be interpreted in terms of γ σ and v σ , . So, we have: The resulting phase-space escape rates are shown in Fig.  3 for the six cases described above. Similar calculations can be done in RFF, where dσ µ = (1, 0, 0, 0) for A, B, C and dσ µ = (0, 1, 0, 0) for D, E, F, leading the following values for W (p): For these cases, A, B, C, D, E, F, in RFF the resulting phase-space escape rates are shown in Fig. 4.
Figs. 3 and 4 show that the momentum dependence of the escape rate, uniform in point A, becomes different at different points of the FO hypersurface, but this change is continuous, when we are crossing the light cone, from point C to point D. Although in RFF, Fig. 4, it seems that there is a principal difference between space-like and time-like FO directions, due to the cut-off Θ(p µ dσ µ ) function, but this is only the consequence of the chosen reference frame, i.e. RFF is defined in a way to stress the difference between these two cases, since going from C to D, the normal vector has a jump, i.e. dσ µ = (1, 0, 0, 0) goes over to dσ µ = (0, 1, 0, 0). Nevertheless, W (p) is a continuous function as we change dσ µ , and in other frames, for example in RFG, Fig. 3, we can see this clearly.  C. The updated simple kinetic model Now, using the new invariant escape rate, eq. (7), we can generalize the simple model presented in [10,11,12], i.e. eqs. (2), for a finite space-time FO layer: Solving the first equation we find for the interacting component: Now, inserting this result into the second differential equation, from eqs. (10), we obtain the FO solution, which describes the momentum distribution of the frozen out particles: As s tends to L, i.e. to the outer boundary of the FO layer, this distribution, depending on the direction of the normal vector (space-like) or time-like will tend to the (cut) Jüttner distribution, f J (p). This means that (part of) the original Jüttner distribution survives even when we reach the outer boundary of the FO surface. To remedy this highly unrealistic result, in [10,11,12,13], rethermalization in the interacting component was taken into account via the relaxation time approximation, i.e. we insert into the equation for the interacting component a new term, which describes that the interacting component approaches some equilibrated (Jüttner) distribution, f eq (s), with a relaxation length, λ 0 : Let us concentrate on the equation for the interacting component. Here the first term from eq. (13), related to FO, moves the distribution out of the equilibrium, and decreases the energy-momentum density and baryon density of the interacting particles. The second term from eq. (13), changes the distribution in the direction of the thermalization, while it does not effect the conserved quantities. The relative strength of the FO and rethermalization processes is determined by the two characteristic lengths, λ and λ 0 .
In general the evolution of the interacting component can be solved numerically or semianalytically, at every step of the integration. Then, the change of conserved quantities due to FO should be evaluated using the actual distribution, f i (s, p) at the corresponding point s. For the purpose of this work, namely for the qualitative study of the FO features, it is enough to use an approximate solution, similarly as it was done in [11,12,13]. This would also allow us to make a direct comparison with results of these older calculations. Thus, the evaluation of the change of the conserved quantities is done analytically, i.e. f i (s, p) is approximated with an equilibrium distribution function f eq (s) with parameters, T (s), n(s), u µ (s).
This approximation is based on the fact that in most physical situations the overall number of particle collisions vastly exceeds the number of those collisions, after which a particle leaves the system or freezes out. This allows us to take that rethermalization 4 happens faster than the freeze out, i.e. that λ 0 < λ or λ 0 ≪ λ. Of course, this argument is true only at the beginning of the FO process, when the density of the interacting particles is still large. When s is close to L, i.e. near the outer hypersurface, the first term in eq. (13) becomes more important than the rethermalization term because of its denominator, but as we shall see in the results section, particles freeze out exponentially fast and for large s, when say 99% of the matter is frozen out, the error we introduce with our approximate solution can not really affect the physical situation. For illustration, let us take a test function, f eq (s) = e −s 2 , (we ignore p dependence for the moment), which is a smoothly and fastly decreasing function 5 . In Fig. 5 4 The words "immediate rethermalization" used in a few earlier publications, were badly chosen, misleading and inappropriate. 5 In real calculations the s-dependence of feq(s, p) is calculated from the energy, momentum and baryon charge loss of the in-we show the numerical solutions for the interacting component for λ 0 = λ and λ 0 = 0.1λ. The results show that for the latter case we can safely take for the f i the approximate solution [11,12,13]:

D. Conservation Laws
The goal of the freeze out calculations is to find the final post FO momentum distribution, and then the corresponding quantities defined through it, starting from the initial pre FO distribution. On the pre FO side we can have equilibrated matter or gas. Its local rest frame defines the RFG, see Fig. 6. We can also define the reference frame, which is attached to the freeze out front, namely the RFF, see Fig. 7. These choices are usually advantageous, but other choices are also possible. Furthermore, the conservation laws and the nondecreasing entropy condition must be satisfied [10]: where [A] = A−A 0 . The pre FO side baryon and entropy currents and the energy-momentum tensor are denoted by N µ 0 , S µ 0 , T µν 0 , while the post FO quantities are denoted by N µ , S µ , T µν .
The change of conserved quantities caused by the particle transfer from the interacting matter to the free matter can be obtained in the following way. For the conserved particle 4-current we have teracting component, where these losses are determined by the momentum dependence of the escape rate and the actual shape of f i (s, p), as discussed here.
Then, using the kinetic definition of the particle current, together with eqs. (14,15), we obtain Similarly the change in the energy-momentum is The parameters of the equilibrium (Jüttner) distribution, f eq (s), have to be recalculated after each step, ds, from the conservation laws as in [11,12,13]. The change of flow velocity, du µ i (s), can be calculated using Eckart's or Landau's definition of the flow, i.e. from dN µ i (s) or dT µν i (s) correspondingly. Then the change of conserved particle density is given by and for the change of energy density we have The change of the temperature of interacting component can be found from this last equation, eq. (21), and from the Equation of State (EoS). This closes our system of equations.
If we fix the FO direction to the x-direction, then eqs. (18,19) can be rewritten as: where, the four momentum of particles is p µ = (p 0 , p), p = |p|, p x = p cos θ p , the flow velocity of the interacting matter is u µ

E. Changes of the conserved current and energy-momentum tensor
In this section we show new analytical results for the changes of the conserved particle current and energymomentum tensor. The formulae are analogous to those from ref. [10,11], but now they are calculated with the Lorentz invariant angular factor from eq. (7). We show results for both massive and massless particles.

III. RESULTS AND DISCUSSION
In this section we calculate the post FO distributions and compare the results to former calculations presented in [10,11,12]. The effect of two main differences due to the new Lorentz invariant escape rate, eq. (7), is to be checked: -the infinite, (∞), FO layer or finite, (L), FO layer, -the simple angular escape rate, P , or the covariant escape rate, W .
We performed calculations for a baryonfree massless gas, where we have used a simple EoS, e = σ SB T 4 , where σ SB = π 2 10 . The change of temperature is calculated based on this EoS and eq. (21). There are no conserved charges in our system, consequently we use Landau's definition of the flow velocity [11]: where ∆ µν i (x) = g µν − u µ i (x) u ν i (x) is a projector to the plane orthogonal to u µ i (x), while e i (x) and P i (x) are the local energy density and pressure of the interacting component, i.e. T µν i (x) = (e i (x) + P i (x))u µ i (x)u ν i (x) − P i (x)g µν . A detailed treatment of Eckart's flow velocity can be found in [10,11]. For such a system we finally obtain the following set of differential equations: We will present the results for four different cases: P ∞ : We use the simple, but relativistically not invariant angular factor, cos θ p , in the escape rate, The system is characterized by an infinite FO length (up to x max = 300λ in calculations). The results are shown in Figs. (8,10,12). This is the same model as in [10,11,12].
P L : Next, we are using the simple angular factor, but in this case inside a finite FO layer, L = 10λ, The results are shown in Figs. (9,11,14).
W ∞ : Then, we are dealing with the new Lorentz invariant angular factor in the escape rate and with an infinite FO length, x max = 300λ, The results are shown in Figs. (8, 10, 12, 13).
W L : Finally we present the primary results of this paper, using both our new improvements, i.e. the covariant escape rate of eq. (7), The results are shown in Figs. (9,11,14,15).
In the case of FO in the infinite layer the factor, L/(L−x), was replaced by 1. We presented the situation at a distance of x max = 300λ, where the amount of still interacting particles is negligible. For the particular cases when we are dealing with an infinite FO, i.e. P ∞ and W ∞ , or with finite layer FO, i.e. P L and W L , the results are plotted together. Thus, in one figure the focus is on the consequences caused by the different angular factors.
Thin lines always denote the cases with simple relativistically not invariant angular factor, these correspond to P ∞ and P L . Thick lines always correspond to cases with covariant angular factor, W ∞ and W L .
All the figures are presented in the RFF.
A. The evolution of temperature of the interacting component The first set of figures, Figs. 8, 9, shows the evolution of temperature of the interacting component, in fact the gradual cooling of the interacting matter, for the different cases, P ∞ , P L , W ∞ and W L .
First, on all figures matter with larger (positive) flow velocity, v 0 , cools faster. This is caused by the momentum dependence of the escape rate, which basically tells that faster particle in the FO direction, will freeze out faster. Thus, the remaining interacting component cools down, since the most energetic particles freeze out more often than the slow ones. Of course, for larger initial flow velocity, v 0 , in the FO direction, there are more particles moving in the FO direction with higher momenta in average, than for a smaller flow velocity. Now, comparing Fig. 8 with Fig. 9, we can see the difference between finite and infinite FO dynamics. In a finite layer the cooling of interacting matter goes increasingly faster as FO proceeds, while for FO in infinite layer the cooling gradually slows down as x increases. The reason is the factor, L/(L − x), which speeds up FO as L − x decreases, and forces it to be completed within L.  The difference between P and W escape rates comes from the denominator, which is p 0 in case of P and p µ u µ in case of W . This difference leads to a stronger cooling for the escape rate P which is bigger than W , if v 0 = 0. This can be seen well at later stages of infinitely long FO, Fig. 8, particularly for the positive initial flow velocity. In all other cases the difference between old and new an-gular factors is insignificant, what supports our "naive" generalization of the angular factor.  Comparing Fig. 10 with Fig. 11, we can see again that in a finite layer the flow velocity decreases faster and faster as FO proceeds, while for FO in infinite layer the velocity change gradually slows down as x increases. The reason is the L/(L − x) factor, as discussed above.
The difference between the evolution of the flow velocity, due to the different angular factors, is again not significant, supporting its generalization.
C. The evolution of the transverse momentum and contour plots of the post FO distribution The next set of figures, Figs. 12, 13, shows the evolution of the transverse momentum distribution, while Figs. 14, 15, present the contour plots of the post FO momentum distribution, for W ∞ and W L . We have presented a one-dimensional model here, but we assume that it is applicable for the direction transverse to the beam in heavy ion experiments. The presented plots should be qualitatively compared to the transverse momentum distributions of measured pions. The local transverse momentum (here px) distribution for a baryonfree massless gas at (py = 0), calculated with the two escape rates: P∞ (thin lines) and W∞ (thick lines), for an infinitely long FO, (xmax = 300λ). The initial parameters are v0 = 0 and T0 = 170 MeV. The transverse momentum spectrum is obviously curved due to the freeze out process. The slope of the transverse momentum distribution increases as we are approaching infinity.
What we see is that all the final post FO momentum distributions are essentially the same. This is very important outcome from our analysis, which we will discuss below. Also, one can see that resulting post FO distributions are non-thermal distributions, as it has been shown already in [10,11,12], they strongly deviate from exponential form in the low momentum region. The increase in the final FO spectra over the thermal distribution for low momenta is connected to the fact that at late stages of the FO process, the interacting component is cold and its flow velocity is negative. So, it contributes only to the low momentum region of the post FO spectra. These results were obtained in a stationary onedimensional model with a single flow velocity. In reality different space-time sections of the overall FO layer are moving with respect to each other with considerable velocities, i.e. v ≈ 0.2 − 0.7. Therefore, the superposition of these parts of the FO layer wash out the very sharp peaks at small momenta, while the curvature at higher momenta, although it is smaller, may persist even after superposition. There are several effects mentioned in the literature, which can cause such a curvature. The effects discussed in this section, arising from kinetic description, may contribute to the curvature of the spectra, but we need a more realistic full scale, nonstationary 3dimensional model to estimate the expected shape of the p t spectra in measurements. Consequently, both the contributions of space-like and time-like sections of the FO layer have to contribute.

D. Freeze out in layers of different thickness
In this section we show the results of calculations performed with W L escape rate, for different finite FO layer thicknesses. Some results of such an analysis have also been presented in [26].
In Figs. 16, 17, we present the evolution of the temperature and flow velocity of the interacting component for L = 2λ, 5λ, 10λ, 15λ. We plot the resulting curves as function of x/L, what allows us to present them all in one figure. We clearly see, and this agrees also with our previous comparison to infinitely long FO, that by introducing and varying the thickness of the FO layer, we are strongly affecting the evolution of the interacting component.
We can also study how fast the energy density of the interacting component is decreasing, see Fig. 18. Since there is no expansion in our simple model, the evolution of the energy density is equivalent to the evolution of the total energy of the remaining interacting matter. We can see that the decrease of the energy density of the interacting component is exponentially fast, what justifies our way of getting approximate an solution for the interacting component, see section II C.  seen already in the previous section, the final post FO distribution, which is in fact the measured quantity!

IV. CONCLUSIONS AND OUTLOOK
In this paper we presented a simplified, but still nontrivial, Lorentz invariant freeze out model, which allows us to obtain analytical results in the case of a massless baryonfree gas. In addition the model realizes freeze out within a finite freeze out layer.
We do not aim to apply directly the results presented here to the experimental heavy ion collision data, instead our purpose was to study qualitatively the basic features of the freeze out effect, and to demonstrate the applicability of this covariant formulation for FO in finite length.
As it has been indicated in the previous publications [10,11,12], the final post FO distributions are nonequilibrated distributions, which deviate from thermal ones particularly in the low momentum region. The final spectra have a complicated form and were calculated here numerically. In large scale (e.g. 3-dim CFD) simulations for space-like FO the Cancelling Jüttner distribution [27], may be a satisfactory analytical approximation.
Our analysis shows that by introducing and varying the thickness of the FO layer, we are strongly affecting the evolution of the interacting component, but the final post FO distributions, even for small thicknesses, e.g. L = 2λ, look very close to our results for an infinitely long FO, first obtained in [10,11,12].
The results suggest that if the measured post FO spectrum is curved, as shown in Fig. 19, then it doesn't matter how thick FO layer was, and we do not need to model the details of FO dynamics in simulations of collisions! Once we have a good parametrization of the post FO spectrum (asymmetric, non-thermal), it is enough to write down the conservation laws and non-decreasing entropy condition with this distribution function [7], (and probably with some volume scaling factor to effectively account for the expansion during FO). This Cooper-Frye type of description can be viewed from two sides. From experimental side, when we know the post FO spectra, we can extract information about the conditions in the interacting matter before FO. In theoretical, e.g. fluid dynamical, simulations such a procedure would allow us to calculate parameters of the final post FO distributions to be compared with data. In this way our results may justify the use of FO hypersurface in hydrodynamical models for heavy ion collisions, but with proper nonthermal post FO distributions.
At the same time, while the final distribution, f (p), is not sensitive to the kinetic evolution, other measurables, especially the two particle correlation function may be more sensitive to the details and extent of the FO process.
The model can also be applied to FO across a layer with time-like normal. While several of the conclusions can be extended to the time-like case, it also requires still additional studies [14].
For realistic simulations of high energy heavy-ion reactions the full 3D description of expansion and FO of the system should be modeled simultaneously. We believe that our invariant escape rate, can be a basic ingredient of such models. with some difference in the low momenta region, which is shown in more detail in the "zoom in" subplot. The two thick lines correspond to some effective thermal distributions, with the corresponding parameters displayed in the plot legend. These are shown to illustrate the difference between obtained post FO distributions and thermal distributions.