Mapping out the QCD phase transition in multiparticle production

We analyze multiparticle production in a thermal framework for 7 central nucleus nucleus collisions, $e^+$+ $e^-$ annihilation into hadrons on the Z resonance and 4 hadronic reactions (p+p and p+$\bar{p}$ with partial centrality selec tion), with center of mass energies ranging from $\sqrt{s}$= 2.6 GeV (per nucleon pair) to 1.8 TeV. Thermodynamic parameters at chemical freeze-out (temperature and baryon and strangeness fugacities) are obtained from appropriate fits, generally improving in quality for reactions subjected to centrality cuts. All systems with nonvanishing fugacities are extrapolated along trajectories of equal energy density, density and entropy density to zero fugacities. The so obtained temperatures extrapolated to zero fugacities as a function of initial energy density $\epsilon_{in}$ universally show a strong rise followed by a saturating limit of $T_{lim}$ = 155 $\pm$ 6 $\pm$ 20 MeV. We interpret this behaviour as mapping out the boundary between quark gluon plasma and hadronic phases. The ratio of strange antiquarks to light ones as a function of the initial energy density $\epsilon_{in}$ shows the same behaviour as the temperature, saturating at a value of 0.365 $\pm$ 0.033 $\pm$ 0.07. No distinctive feature of 'strangeness enhancement' is seen for heavy ion collisions relative to hadronic and leptonic reactions, when compared at the same initial energy density.


Introduction
Hadronic reactions involving copious production of secondary particles have been associated with an underlying thermodynamic behaviour since the earliest observations in cosmic rays [1]. The observable energy regime ranging up to 140 TeV per nucleon pair corresponding to incident cosmic ray nucleons of E ≤ 10 10 GeV demands elaborate simulations of the observed extended air shower development , in order to extract definite multiplicity distributions of the elementary hadronic and nuclear reactions [2]. Thermodynamic models are widely and successfully used to describe identified particle ratios in hadronic and especially heavy ion collisions [3] , [4], [5], [6], [7], but were also extended to e + e − annihilation into hadronic final states [8] , [9]. In heavy ion collisions two main parameters -energy per nucleon pair and centrality -are varied and their influence on thermodynamic variablestemperature and chemical potentials -are studied. We derive and discuss thermal properties in search of the phase boundary between quark gluon plasma and condensed hadrons. We consider as one characteristic signature of this boundary the critical dependence of kaon number densities on the initial energy density. The kaon multiplicities observed in Pb+Pb collisions at √ s=17 GeV [10], and in other nucleus+nucleus and p+p collisions at √ s ∼ 5-19 GeV [11,12], serve as motivation for our thermal description extending the previous analysis of one of us (SK) [11,12] (see figure  8).
In section 2 we discuss the dependence of the critical temperature on the vaccum pressure in QCD, limited to the case of vanishing chemical potentials, extending the work of reference [13]. In section 3 we seek to assign each thermodynamic state established at chemical freeze-out for any given reaction for which chemical potentials for baryon number and strangeness do not vanish, an equivalent state at zero fugacities. We do this extrapolating along curves of equal entropy density, energy density and number density. The extrapolation of the temperature of systems with finite fugacities to zero fugacities, and the investigation of their dependence on the initial energy density has been proposed in [12]. We find that a more universal parameter as a measure of positive and negative strangeness production is λ s = 2 s u + d . This parameter is also extrapolated to the zero fugacity systems. The details of this procedure are worked out in the first four subsections of section 3. The results at zero fugacities are contained in subsection 3.5 and represented in figures 2 -8 in relation with the energy density initially achieved in each reaction. Despite large errors, the phase boundary between quark gluon plasma and hadronic phases is clearly mapped out.
We use in the hadronic phase the approximation of noninteracting hadron resonances to describe in this sense global ratios of hadrons produced in the following reactions : 1) central Au+Au collisions at RHIC √ s = 130 A GeV [14] , [15] .

Outline of basic points
We consider first the thermodynamic potential Φ of a grand canonical ensemble of hadron resonances without further interactions among them. The thermodynamic variables are i.e. volume, temperature, baryon and strangeness fugacity. The fugacity of the third component of isospin is set to zero, neglecting isospin asymmetries. Φ = g V ; g = g ( β , χ b , χ s ) = α g α + g 0 β = 1 / T ; g = β p ; p : pressure (2) In eq. (2) the sum extends over hadron resonances denoted by α, where we only include the pseudoscalar and vector u,d,s meson nonets as well as the spin 1/2 baryon octet and spin 3/2 decuplet and their antiparticles as well as the f 0 (400 − 1200) or σ, interpreted as scalar glueball [27] for simplicity. g 0 = β p 0 takes into account the nonzero vacuum pressure characterizing the hadronic phase in QCD, which we restrict to the three light flavors.
The temperature variation of the vacuum pressure is an approximation with T cr = T cr ( χ b , χ s ). The quantities g α in eq. (2) for noninteracting resonances are then given by In eq. (5) ( I , Sp , B , S ) α denote isospin, spin, baryon number and strangeness of the hadron resonance α. The -,+ signs apply to bosons and fermions respectively. If we take into account the variation of the vacuum pressure with temperature then the masses m α → m α ( T ) of quasiexcitations become temperature dependent quantities. In all subsequent calculations we neglect these temperature dependent effects, which however do set in dramatically, when the temperature deviates from the crirtical one by less than 10 %.
We state here, that a thermal description cannot account for any azimuthal ( ϕ -) dependence of inclusive cross sections. Thus azimuthal anisotropies provide a measure for non-equilibration. The potential Φ in eq. (2) and the entropy denoted by S give rise to the differentials In eq. (7) ̺ s denotes the entropy density, ε the energy density and ̺ b,s the net baryon and strangeness (number) density respectively.
As is apparent from eq. (7) the entropy density obtains no contribution from vacuum energy density and pressure, as long as the vacuum pressure is taken as independent of temparature as indicated in eq. (6). Energy density and total number density then take the form The approximation of hadronic interactions by free hadron resonances is of course at best a consistent approximate thermodynamic description. We will use it here to describe in this sense global ratios of hadrons produced in the reactions listed in the introduction. We will present a detailed analysis of these collisions, as seen at the time of chemical freeze-out in the next section. The key parameter characterizing not this freeze-out instance but the time of first onset of thermal conditions is the initial energy density ε in . It determines the thermal parameters of the hadronic system as seen at chemical freeze-out and distinguishes those systems which have hadronized after transcurring the quark gluon plasma phase from those which remained throughout in the hadronic phase. The phase boundary for vanishing fugacities In the remainder of this section we analyze the equilibrium condition for both phases to coexist in the case of vanishing fugacities or equivalently chemical potentials. This latter choice is just for clarity of argument and also because for limiting center of mass energies and finite baryon number, strangeness and charge of the initial system this is the limiting case. This condition corresponds to equating the Gibbs densities of the two phases In the same approximative vein as applied to the hadronic phase (H) we approximate the various quasiexcitations in the quark gluon phase by free gluons, quarks and antiquarks analogous to the form of g α in eq. (5) In eq. (11) the color/spin weights, masses and boson/fermion signs are specified. The quark masses represent a renormalization group invariant 'best' choice. It turns out that near the critical temperature T cr = O ( 200 MeV ) the (anti)quark contributions to the Gibbs density do not deviate substantially from the limit of vanishing quark masses including m s . The crossings of g H and g QGP as displayed in figure 1 seem to suggest a first order phase transition with respect to energy density. This is however due to the approximation of fixed masses for hadron excitations and of free quark and gluon flavors in the region of T cr = 194 ± 18 MeV. We note that the above estimate of the critical tempeature at zero fugacities confirms within the theoretical error the previous estimate of one of us [13].
Within the approximations and errors we obtain the critical temperature and critical hadronic energy density on the side of the hadronic phase for vanishing fugacities Obviously the critical boundary curve could be extended to arbitrary chemical potentials, but again the 'melting' of masses of hadronic quasiexcitations modifies the Gibbs density on the side of the hadronic phase, whereas the approximation of free quark and gluon modes does not warrant the extrapolation to arbitrary chemical potentials, especially for small temperature.
3 Reduction of chemical freeze-out parameters to zero fugacities in multiparticle production In this section we apply the thermal model introduced in the last section (equation 5) to extract the intensive thermodynamic parameters defined in eq. 1. We then extrapolate these parameters to zero fugacities along states with equal entropy density, energy density or density. The goal of this analysis is to compare the temperature at zero fugacities with the initial energy density achieved in the collision in order to reveal a boundary reflecting the QCD phase transition.
One serious well known problem when comparing models to experimental data is due to decays of resonances. This effect is called 'feeding'. We compare calculated particle ratios to experimental data, taking all (strong, electromagnetic and weak) decays with a branching ratio not below 1% into account. We try to account for experimental acceptance for K + and K − , assuming a 50% feeding to pions, due to their long decay length. The other weak decays (e.g. Λ, K 0 s ) have a much shorter decay length and are assumed to fully feed into secondary particles. The K 0 l decay is not considered. Contrary to references [4,5] we do not take all hadronic resonances below 2 GeV mass into account. The theoretical error inherent to the free resonance approximation allows in our opinion our reduced set. This is supported by the finding that we derive similar thermodyamic parameters with similar accuracy as e.g. the ones found in reference [16]. Though the high mass resonances are definitely produced at sufficiently high energy, it is not clear that the thermal description becomes more accurate by including them because the quasiexcitations may not correspond to them, especially if these are quark and gluon modes. All similar thermal descriptions suffer from replacing interactions by noninteracting resonances. Furthermore, the isotropic angular distributions characterizing the models are far from reality [30,31] (e.g. flow phenomena, diffractive phenomena).
We enforce strangeness conservation, to obtain the best fit, exept in the case of p+p collisions at √ s=17 GeV. We find the error on the fitted parameters, at ± 1 unit of the ratio of χ 2 /DOF away from the best fitted value. Our quoted error reflects the inherent theoretical error, despite the difficulty to quantify it. Note that other analyses [4,5] finding a smaller error seem to use as error the fitted parameters at ± 1 unit of the χ 2 itself.
In the following sections we note if the data we use are taken with a minimum bias trigger (that means no trigger bias is imposed) or a central respectively a peripheral trigger. Central is here understood as the lowest impact parameter region. It is usually selected in the experiment by considering only collisions with the largest particle multiplicity (e.g. CERN-WA97, Fermilab-E735) or the largest transverse energy (e.g. CERN-NA52), or the smallest forward going energy (which reflects the spectator nucleon number) (e.g. CERN-NA49).
Respectively, peripheral collisions are a selection of collisions with the largest impact parameters using the same means as described above (charged multiplicity, energy).
In the following analysis we also give the quantity: λ s = (2s) (u+d) which we take as a measure of the strangeness suppression factor defined and used in the literature as λ s = 2(s+s) (u+u+d+d) [5]. We use antiquarks to consider in a simplified way only the newly produced valence antiquarks. We give also the equivalent λ s at zero fugacities. From quark counting rules one expects λ s to be in the range 0.3-0.5 [32]. In principle one should include in any definition of λ s the newly produced sea quarks however the relevant proportion of the latter inside hadrons, is still subject of experimental and theoretical investigations.
Generally we find that thermal fits to minimum bias or peripheral nuclear and hadronic reactions are not satisfactory, in contrast to central collisions. We conjecture that this is due to the presence of at least two thermal sources which we attribute to diffractive versus pomeron induced subprocesses. Diffractive processes at finite √ s feed back into the midrapidity region, whereas pomeron induced subprocesses result in particle production at midrapidity dominantly. Following this conjecture, we conclude that the relative importance of diffractive to pomeron induced subprocesses decreases with increasing centrality of the collision. In this respect it would be helpful firstly, to increase experimental coverage in the target and/or the projectile fragmentation regions. This study has been proposed by [33,34]. Secondly, it is important that experiments apply centrality cuts, in nucleus+nucleus as well as in p+nucleus and in elementary particle collisions.

Nucleus nucleus reactions
In the following 4 subsections we analyse data from nucleus nucleus collisions at 1. RHIC, 2. SPS, 3. AGS and 4. SIS energies.

Central Au+Au collisions at √ s=130 GeV
We used 3 measured ratios at midrapidity from [14,15] : Λ/Λ, p/p and h − / Ncharged and imposed strangeness conservation. We exclude weak decay products as is done in the experiment [14]. The data from reference [15] are preliminary. The predicted and the experimental ratios are shown in table 1. The resulting χ 2 /DOF is (1.41/1) (CL 23%) for a temperature of 168 ± 40 MeV. When adding a systematic error of 15% linearly to account for the experimental feeding uncertainties and the fact that the ratios are measured at midrapidity only, the resulting χ 2 /DOF is 0.27/1 (CL 60%). The errors on the temperature and λ s are not well estimated in this analysis due to the insensitivity of the ratios used, to the temperature at fixed fugacities. We estimate the error on the temperature searching for a variation of the χ 2 /DOF by 1, while changing the fugacities and fixing each time the T by imposing strangeness conservation. In order to improve the fit quality the addition of other ratios more sensitive to temperature variations is needed, like K/π etc., which will soon be measured by the experiments at RHIC. After defining the (T, µ b , µ s ) values describing the particle ratios produced in central Au+Au collisions at √ s=130 GeV at the chemical freeze-out we extrapolate to the T at zero fugacities (table 2) Table 1: RHIC Au+Au at √ s= 130 GeV.
Predicted versus experimental particle ratios for the best fit of our model.  Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities.

Central nucleus nucleus collisions at ∼ 200 A GeV
Here we show results from central Pb+Pb collisions at 158 A GeV, and from S+S and S+A collisions at 200 A GeV.

Central Pb+Pb collisions at 158 A GeV
For central Pb+Pb collisions at 158 A GeV, we use the particle ratios from table I in [16] to compare with our model predictions. We don't consider a chemical potential for isospin [16], therefore we don't consider e.g. the π − /π + ratio. We did not use the ratios in [16] for which feeding correction is included. We did not use the d/d ratio, since there is strong experimental evidence that d and d are formed in the thermal hadronic freeze-out through coalescence [10,35]. Furthermore, the p may not freeze out in the chemical hadronic freezeout due to the high cross section for p + p annihilation to e.g. π as noted in [36]. We introduce a systematic error of 14% (quadratically added) in the ratios including WA97 data, because they are measured at midrapidity only. We add also quadratically) a 10% systematic error to ratios including h − due to uncertainty in the feeding correction. The predicted ratios are shown in table 3 together with the experimental data.
The χ 2 at which strangeness is conserved is 16.1 over 12 degrees of freedom (CL=20 % χ 2 /DOF =1.34) and corresponds to T=162 +10 -28 MeV. The asymmetry of the error (defined as the T deviation 1 unit of χ 2 /DOF away) is due to the fact that the 16.1 is not at the minimum. The minimum of the χ 2 is 11.9 over 11 DOF (CL=50%, χ 2 /DOF =1.08) at T=154 +14 -16 MeV. There is a probable physics process which may induce an imbalance of strangeness in experiment, if the measurements do not include the target and beam fragmentation regions: in particular particles with negative S especially hyperons may be produced near beam rapidity. In ref. [37] it is conjectured that strangeness is not balanced (with missing particles with S negative). However we will demand in the following exact strangeness conservation. After defining the (T, µ b , µ s ) values describing the particle ratios produced in central Pb+Pb collisions at 158 A GeV at the chemical freeze-out we extrapolate to T at zero fugacities. The resulting equivalent temperatures for equal entropy density (ρ s ), energy density (ρ e ), and number density (ρ n ) are shown in table 4. ratio model data Table 3: SPS Pb+Pb at √ s= 17 GeV, most central events.
Predicted versus experimental particle ratios for the best fit of our model.  Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities.

Central S+A collisions at 200 A GeV
In table 5 we give the same quantities for central S+A collisions at 200 A GeV with A=Au,W,Pb. The data used are taken from reference [17]. Here we mainly used the p/p and K + /K − ratios to find (T, µ b , µ s ) while we also suppose exact strangeness conservation. The resulted predicted ratios are p/p = 0.131 and K + /K − = 1.59 to be compared with the experimental data p/p = 0.12 ± 0.02 and K + /K − = 1.59 ± 0.15. The resulting temperature is T=166 +10 -28 MeV. The errors are taken to be percentually the same as in the Pb+Pb system. We do not perform a full fit, as the resulting temperature is within the errors in agreement with reference [17] where a full fit to many ratios is performed. The calculated values in [17]   Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities.

Central S+S collisions at 200 A GeV
In table 6 the results for the central S+S collisions at 200 A GeV are shown, data are taken from references [18]. We use K 0  Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities.

Central nucleus nucleus collisions at ∼ 10 A GeV
Here we show results from central Si+Au collisions at 14.6 A GeV, and from Au+Au collisions at 11.6 A GeV.

Central Si+Au collisions at 14.6 A GeV
In table 7 the results for central Si+Au collisions at 14.6 A GeV are shown. We use the ratios K/π, K + /K − , Λ/Λ and φ/π from [5] as well as strangeness conservation. We did not use p in view of the large annihilation expected at low energy. The resulting χ 2 /DOF =6.8/2=3.4 (CL ∼ 3.3%) imposing strangeness conservation is at T=120 +8 -35 MeV. We note however that the best χ 2 ignoring strangeness conservation (χ 2 /DOF =0.092/2 with CL 63%) is at µ b /T =4.05, µ s /T =0.9 and T=102 MeV.
Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities.
Central Au+Au collisions at 11.6 A GeV The results for central Au+Au collisions at 11.6 A GeV are shown in tables 8, 9. We used K + /K − , K/π, p/π + and K/Λ ratios from [5]. The resulting  Table 8: AGS Au+Au at √ s= 4.9 GeV, most central events.
Predicted versus experimental particle ratios for the best fit of our model. ± 0.08 ± 0.037 Table 9: AGS Au+Au at √ s= 4.9 GeV, most central events.
Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities.
Predicted versus experimental particle ratios for the best fit of our model.  Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities.

Electron positron reactions
The primary qq production in e + e − collisions is not thermal due to the hard nature of the primary γ, Z couplings. However as long as the final particle multiplicity is much higher than 2, it is conceivable that this fact does not matter and the subsequent fragmentation of quark-antiquark jets into hadrons is thermal. We used the K/π, ρ/π, π/p and ∆/π ratios.
Our results for e + e − collisions at √ s=91 GeV are shown in table 12. We used the data from reference [9] in particular the initially produced particles from the Pei model (table 1 in [9]). We used as systematic error of the feeding correction performed in this model, the difference between the Pei and the Jetset model (also in table 1 in [9] [38] and [9]. Our best value of the temperature is however nearer to the results of reference [9].  Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities.

Hadronic reactions 3.3.1 Proton proton reactions
Here we show results from p+p collisions at √ s=17 and 27 GeV.
The bad quality of the fits of p+p collisions at √ s=27 GeV renders the resulting temperatures questionable. Nevertheless we point out that the two analyses yield incompatible temperature values within the errors. The analysis of [23] does not introduce systematic errors as we do. However, it allows for an arbitrary weighting factor (called γ s ) acting only on strange particles. This factor is similar to an increase of the systematic experimental errors of strange particle multiplicities only.  Table 13: p+p collisions at √ s= 27 GeV.
Predicted versus experimental particle ratios for the best fit of our model. The χ 2 /DOF of this fit is not acceptable.  Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities. The χ 2 /DOF of this fit is not acceptable.
Proton proton reactions at √ s=17 GeV We used 6 measured ratios from references [24,25] namely K/π, p/p, π/p, K + /K − , φ/π, K/p and imposed strangeness conservation. We add a systematic error of 15% linearly to account for the experimental feeding uncertainties.
The predicted ratios in table 15 are overestimating K + , K − and φ yields, therefore suggesting that if thermal conditions prevail, then at least two thermal reservoires are present. This may be due to the influence of diffractive processes in minimum bias triggers.
We therefore repeat the fit with only ratios which are not much influenced by diffractive processes namely K/π, K + /K − , and φ/π. Note that the temperature at the minimum of χ 2 does not satisfy strangeness conservation. This is combatible with the hypothesis of the importance of diffractive processes showing up in leading baryons (e.g. p, Λ). Imposing strangeness conservation the fit results in a temperature of 144 MeV with a bad χ 2 /DOF =169/1 (CL ∼ 10 −39 ). At the minimum of the χ 2 the temperatures is 96 +8 -9 MeV with a better χ 2 /DOF =5. 13   Predicted versus experimental particle ratios for the best fit of our model. Only three ratios are fitted (φ/π, K + /K − , K/π). Thermodynamic parameters for the best fit and temperatures and λ s extrapolated to zero fugacities. Only three ratios are fitted (φ/π, K + /K − , K/π).

Proton antiproton reactions
Here we show results from p + p collisions at √ s=900 GeV and central and peripheral p + p collisions at √ s=1.8 TeV .
Central proton antiproton reactions at √ s=900 GeV For p + p collisions at √ s=900 GeV we used 5 measured ratios out of the multiplicities from reference [23]. in particular K 0 s /charged and N/charged,  Table 18: p+p collisions at √ s= 900 GeV.
Predicted versus experimental particle ratios for the best fit of our model.

Peripheral proton antiproton reactions at √ s=1.8 GeV
For the most peripheral p + p collisions at √ s=1.8 TeV we used 2 measured ratios from reference [26] namely K/π and p/π (table 19). We use the only 2 available ratios from experiment, which with zero fugacities and temperature as the only parameter leaves one degree of freedom. The most peripheral collisions are defined as the ones with the lowest measured charge multiplicity. The temperature from the fit is with a χ 2 /DOF of 18.04/1 having a CL of ∼ 2.0E-4, therefore the fit is very bad. We add a systematic error of 0.004 linearly to account for the experimental feeding uncertainties. This error has been estimated from deviations between the ratios found in different experimental runs with p + p collisions at √ s=1.8 TeV [39]. The predicted and the experimental ratios are shown in p/π 6.87E-2 7.40E-2 ± 9.00E-3 We did not add a systematic error.

The initial energy density estimation
To estimate the initial energy density achieved in each collision we use several methods. The first method (α) is based on Bjorken's formula [40]: where (dE T /dy) ycm is the transverse energy (E T = p 2 T + m 2 ) per unit rapidity at midrapidity, R T is the transverse radius of the particle source after a formation time τ taken to be 1 fm/c. We estimate the (dE T /dy) ycm when not available from experiment, as follows. We first estimate the total E T as (Method (α)) where the sum runs over particles of type α with multiplicity per rapidity unit at midrapidity (ycm) (dN/dy) α , mean transverse momentum p T α and mass m α . The mean transverse momentum is also taken at midrapidity, when available. We take three types of particles: π, K and nucleons, as well as antinucleons when available. To calculate the differencial (dN/dy) ycm , when not available, we first find the mean (dN/dy), dividing the total particle multiplicity by the total rapidity interval and then we use an extrapolation factor A to extract the value at midrapidity. This factor is found from measured rapidity distributions of particles at or as near as possible to the √ s considered.
This may vary from nuclear to hadronic reaction and with √ s. In the cases where the dN/dy of the particles are measured at midrapidty we use these values.
The formation time is taken in the literature usually in all reactions 1 fm/c, and this is what we use. In the following we give some examples of our ǫ in estimation with method (α).
1. e + + e − For the e + e − collisions at √ s=91 GeV we used the multiplicities from [9]. For two jet events in e + e − collisions the jet axis is the longitudinal axis defining rapidity and p T refers to it. We used the dN/dy at midrapidity of π, K, p in quark jets produced in e + e − collisions at √ s=91 GeV by the DELPHI collaboration [42].
We used the mean transverse momenta from hadron hadron collisions at high √ s. In particular we use the mean transverse momenta from p + p collisions at 1.8 TeV from [26] at charged multiplicity N∼ 45 (to represent minimum bias values) namely p T π = 0.34 GeV, p T K = 0.5 GeV and p T p = 0.59 GeV.
In hadron collisions the mean transverse momenta do not change much with √ s above 10 GeV (see e.g. figure 2.7 in [22]). The transverse size of the initial hadronic system of quarks and antiquarks produced, is given by the uncertainty principle, as 1 over the average p T . For the transverse radius of e + e − collisions we take therefore the inverse mean p T of 0.34 GeV, giving R T =0.6 fm. The resulting initial energy density is: The systematic error (which in this case is dominated by the uncertainty in the R T determination) is about 50%.
2. p + p at √ s=1.8 TeV For p + p we use the same method, while we use as transverse radius the radius of the nucleon of 0.8 fm. We use again pions, kaons, nucleons and antinucleons, and the mean transverse momenta measured in [26]. We derive the dE T /dy(ycm from the measured charged total multiplicity per unit rapidity at midrapidity, using the ratios and mean transverse momenta measured as a function of charged multiplicity in [26]. No extrapolation factor to midrapidity is needed here. We find 772 GeV/f m 3 . In reference [26] the ǫ in for dN c /dη = 15 of 3 GeV/fm 3 is calculated using pions only, and therefore is lower than our estimate 2 .

p + p
√ s=900 GeV For p + p at √ s=900 GeV we use the multiplicities from [38] and the mean transverse momenta for p+p at √ s=1.8 TeV from [26].
The initial energy density has been estimated by experimenters for some reactions shown here, using an other method (β), namely the Bjorken formula (equation 13) and the measured transverse energy : where E is the total energy measured with e.g. calorimeters and θ is the angle to the incident beam direction. The resulting values, (used in figures 1, 2 and 5 in the following section), are: 6 GeV/f m 3 [43,44], [43,44], ǫ in (Au+Au √ s=4.9 GeV) = 1.
We estimate the maximal initial energy density with a third method (γ), taking the nuclear energy density of two overlapping nuclei 2 ǫ A times the γ factor of the colliding particles in the center of mass minus one: (method (γ)) with γ = ( √ s/2)/m nucleon , and ǫ A,small =0.179 GeV/fm 3 for small nuclei and ǫ A,big =0.138 GeV/fm 3 for large nuclei is the normal nuclear matter density. The value in equation 16 multiplied by the stopping power gives an estimate of the initial energy density.
We apply this method to the Si+Au and Au+Au collisions at √ s ∼ 5 GeV, Ni+Ni collisions at √ s=2.6 GeV, and Pb+Pb collisions at √ s=17 GeV, to calculate the maximal achieved initial energy density, yielding: Few comments: 1. A problem with method (α) is, that in general it underestimates slightly the E t , since not all particles are taken into account.

A general problem is with equation 15 of method (β):
For particles produced near midrapidity E × sinθ is approximately p T not E T , for γ cm >> 1, whereas E × sinθ is approximately m T for particles moving nonrelativistically in the lab frame.
3. A problem with method (γ) is, that it assumes geometrical and not dynamical compression, possibly slightly underestimating the ǫ in .

5.
We find however different ǫ in values for the low √ s reactions Si+Au and Au+Au at 5.4 and 4.9 GeV, for each one of the above methods with a scatter of a factor of two or more. Also the maximal estimated ǫ in with method (γ) is lower than the results from methods (α) and (β). We therefore conclude that the Bjorken estimate may not be adequate for low energy reactions. This estimate was not meant to be used in the nonrelativistic regime. The method (γ) seems more adequate for these systems. Nevertheless we show two distinct cases in the following figures.
Using method (γ) the ǫ in is slightly higher for Si+Au at 5.4 GeV than for Au+Au at 4.9 GeV √ s, unlike the results of methods (α) and (β). The stopping power which is missing in the equation 16, can however hardly be very different for these two reactions, but may have an influence on the ordering of the ǫ in in the above two reactions.
Systematic error estimation summary Comparing results from methods (α) and (γ) we estimate the systematic error in Si+Au and Au+Au collisions at √ s=5. 4 and 4.9 GeV to be about 50%, mainly resulting from the difficulty in applying the Bjorken formula. For the higher energy nucleus+nucleus collisions the systematic error is smaller ∼ 30% [11] resulting from the estimation of E T . The systematic error for e + + e − is about 50%, resulting from uncertainties in the transverse area definition. The systematic error for p + p colllisions at √ s=900 GeV results from the E T definition and the extrapolation to midrapidity and is estimated to be about 50% . For p + p collisions at the Tevatron, the uncertainty comes mainly from not taking all particles into account with method (α) and is about 30%. The systematic error on the ǫ in estimate for Ni+Ni at √ s=2.6 GeV is 50%, resulting from comparison of our ǫ in value to model calculations [46].

3.5
The combined results: T, λ s at zero fugacities Figures 2 and 3 show the resulting temperature and the λ s factor extrapolated to zero fugacities along isentropic paths, as a function of the initial energy density. The latter is taken from the experimentally measured transverse energy, when available, that is, for Si+Au √ s=5.4 GeV, Au+Au √ s=4.9 GeV, S+S √ s=19 GeV, S+A √ s=19 GeV and Pb+Pb √ s=17 GeV. For the remaining colliding systems, where no initial energy density estimation is available, we use our method (α) based on equation 14, exept for the Ni+Ni system where we use the method (γ). The fact that the initial energy density estimation in figures 2 and 3, was not done with the same method, introduces an additional systematic error on the ǫ in scale. To reduce this uncertainty and show the systematic error of the initial energy density estimation we use in the following the ǫ in estimated for 1.) high √ s (>10 GeV) by equations 13 and 14, and for 2.) low √ s (<10 GeV) by equation 16. Figures 4 and 5 show therefore the temperature and the λ s factor again, extrapolated to zero fugacities along isentropic paths, as a function of the initial energy density, with a different estimation of the latter. In particular, the initial energy density is not taken from the experimentally measured transverse energy, but is estimated using our method (α) based on equation 14, with the exeption of the data at low √ s (< 10 GeV), that is, for Ni+Ni at √ s=2.6 GeV, Si+Au at √ s=5.4 GeV and Au+Au at √ s=4.9 GeV, where we used equation 16, which as discussed in the previous section, seems more adequate for the low energies as the Bjorken estimate. Figure 8 (a) shows the number density of kaons as a function of the initial energy density from reference [12], while in figure 8 (b) the temperature extrapolated to zero fugacities along isentropic paths, is shown as a function of the initial energy density. In figure 8 (a), the collision systems p+p at √ s=17 GeV, S+S at √ s=19 GeV, Pb+Pb at √ s=17 GeV, Au+Au at √ s=4. In particular, the initial energy densities for all heavy ion systems is taken from the experimental calorimetric measurements (method β above) when these are available. For Au+Au at √ s=130 GeV we use our estimate with method (α). For p+p at √ s=17 GeV we used the estimate from reference [12] of ǫ in (p+p at √ s=17 GeV) from [12]=0.85 GeV/fm 3 This is found in [12] by a) estimating the dependence of the ǫ in for Pb+Pb collisions at √ s=17 GeV on the number of participant nucleons N and b) extrapolating this function to N=2.
The rise and subsequent saturation seen in the kaon number density below and above ǫ=1.3 GeV/fm 3 ( figure 8 (a)), shows a clear relation to the same behaviour seen in the temperature ( figure 8 (b)) as a function of the initial energy density.    Figure 4: Temperature at chemical freeze-out and for zero fugacities as a function of the initial energy density for several nucleus+nucleus, hadron+hadron and lepton+lepton collisions. The initial energy density has been estimated in a different way as the one in figures 2 and 3 (see text). We demand for the fits confidence level > 1%.  Figure 6: Temperature at chemical freeze-out and for zero fugacities as a function of the initial energy density for several nucleus+nucleus, hadron+hadron and lepton+lepton collisions. The initial energy density has been estimated in a different way as the one in figures 2 and 3 (see text). We demand for the fits confidence level > 10%.  : (a) Kaon yield per interaction over the source volume at the thermal freeze-out as a function of the initial energy density (ǫ in ) from reference [11,12].
(b) Temperature at chemical freeze-out and for zero fugacities as a function of the initial energy density for several nucleus+nucleus, hadron+hadron and lepton+lepton collisions. All collision systems shown in figure (a), are displayed in figure (b) using the same initial energy density as defined in [11,12]. We demand for the fits confidence level > 1%.

Conclusions
We have discussed in section 2 the dependence of the critical temperature in the QCD phase transition on the vacuum pressure including in addition modes of noninteracting hadron resonances for u,d,s flavors of light quarks as defining the hadronic phase. Equating the pressure of this hadronic phase to the pressure of the plasma phase, represented by noninteracting u,d,s quarks and antiquarks and gluons as shown in figure 1, we obtain for zero fugacities T cr = 194 ± 18 MeV . We have performed a thermal analysis of yields in multiparticle production for 13 reactions summarized in table 21 and discussed in detail in section 3. The intensive thermal parameters ( temperature and baryon as well as strangeness fugacity ) together with an error estimate are used to extrapolate the state associated with the chemical freeze-out of each reaction, studied along curves of equal entropy-, energy-and number-density, to zero fugacities. We represent the so obtained states by two intensive parameters T and λ s = 2 s u + d , i.e. temperature and the ratio of antistrange to nonstrange (valence) antiquark abundances. These two quantities as functions of the initial hadronic energy density achieved in each reactions are displayed in figures 2 -5. The initial energy densities are estimated in subsection 3.4. The resulting error (horizontal error in figures 2 -8) is of the order of 50% around ǫ crit , for the e + +e − collisions at √ s=91 GeV and for p + p collisions at √ s=900 GeV, and approximately 30% at higher ǫ in .
The reactions with an estimated confidence level above 10% were retained, which fall within the errors into two groups :

central Au+Au collisions at
√ s = 4.9 GeV 9. Ni+Ni collisions at GSI √ s = 2.8 A GeV Both quantities T and λ s saturate to constant values T lim = 155 ± 6 ± 20 MeV (syst) ; λ lim = 0.365 ± 0.033 ± 0.07 (syst) (17) above the dividing energy density ε in ∼ 1 GeV /f m 3 . Within the errors this is compatible with the critical energy density ε crit ∼ 1GeV /f m 3 obtained in lattice QCD [47] as well as with the results on the critical parameters obtained in section 2 ( eq. 12 ). This latter agreement and the extension of the work in [13] to include excited hadronic degrees of freedom is a new result.
This limiting temperature is expected to be somewhat below the critical one, because hadronization is not an instantaneous process. The mean values show a difference in temperature of 40 MeV but the errors mainly due to the approximations of free quasiexcitations and experimental errors are large adding linearly also amount to ∼ 40 MeV.
We estimate the systematic error on the limiting T , λ s values, fitting the high ǫ in region of the figures 4, 5 using once a line fit and once a horizontal line fit, varying in both cases the number of bins taken, and estimating the deviations of the mean values resulting from each of the two fits. The error is ∼ 20 MeV for the temperature and 0.07 for the λ s .
We interpret this saturation as characteristic for all reactions which before chemical freeze-out have passed through the quark gluon plasma phase, as opposed to those which remained throughout in the hadronic phase.
We note that due to the large errors the dividing line as drawn at ε in ∼ 1GeV /f m 3 is subject to a corresponding uncertainty. The saturation phenomenon observed here is in our interpretation confirmed by the same phenomenon derived for the kaon number density as a function of the initial energy density by one of us [11,12] (figure 8).
The agreement of the theoretical thermodynamic model with experimental ratios of identified particle multiplicities is quantified in the errors as discussed in section 3. We emphasize that agreement cannot be perfect due to several obvious and less obvious facts : a ) The rapidity and transverse momentum distributions do not agree with the spherically symmetric distributions in the thermodynamic description.
b ) There may well exist kinematically distinct thermal systems, as is indicated by the phenomenon of longitudinal and transverse flow [31].
c) The target and projectile diffractive regions may form distinct further thermal subsystems characterized by a different temperature than the particles produced near midrapidity [48].
The results derived here do indicate in our interpretation a high degree of thermalization of the systems studied. This implies the clear indication, that the quark gluon plasma phase is part of the hadronization process characteristic for hadronic, leptonic and nuclear reactions above the critical energy density. The universal excitation of the quark gluon plasma phase in nuclear and hadronic as well as leptonic multiparticle production is a new result.
The behaviour of strangeness production as reflected by the parameter λ s also reveals a new aspect: the saturation phenomenon corresponds to a strangeness enhancement relative to states with low initial energy density. No 'strangeness enhancement' is seen for heavy ion collisions relative to hadronic and leptonic reactions, when compared at the same initial energy density and zero fugacities. This enhancement arises only if systems with very different thermal properties in particular different fugacities are compared to one another at the same √ s, therefore at different energy density when nucleus nucleus collisions are compared to elementary particle collisions.
LHC results for central p+p and Au+Au collisions at √ s= 14, 5.5 TeV will be important for the confirmation of the above picture. We conclude by proposing further experimentation in the following areas : i ) study of p+ p collisions at the Tevatron, with centrality selection on transverse energy, similar to nuclear collisions, and possibly with an improved particle identification option including the fragmentation regions.
ii ) in conjunction with i) an extension to compare open and closed c+ c production.
iii) within the program of heavy ion collisions at RHIC and at the SPS, allready under way, a dedicated study in the neighborhood of ǫ in,crit = 1 GeV/fm 3 , e.g. for ǫ in =(0.6-1.5) GeV/fm 3 , in nucleus+nucleus, p+A and p+p collisions, using centrality selection.