Abstract

Topics on hydrodynamic model of nucleus-nucleus collisions

Y. Hama^I; T. Kodama^II; O. Socolowski Jr.^I

^IInstituto de Física, Universidade de São Paulo, C.P. 66318, 05315-970 São Paulo-SP, Brazil

^IIInstituto de Física, Universidade Federal do Rio de Janeiro, C.P. 68528, 21945-970 Rio de Janeiro-RJ , Brazil

ABSTRACT

A survey is given on the applications of hydrodynamic model of nucleus-nucleus collisons, focusing especially on i) the resolution of hydrodynamic equations for arbitrary configurations, by using the smoothed-particle hydrodynamic approach; ii) effects of the event-by-event fluctuation of the initial conditions on the observables; iii) decoupling criteria; iv) analytical solutions; and others.

1 Introduction

Hydrodynamic Model has been proposed by Landau [1] in 1953 as an improvement over the Fermi statistical model [2] for the multiple particle production phenomena in high-energy nuclear collisions. At that time, these phenomena were observed in cosmic rays. Although the Fermi model offered an ingenious insight into the mechanism of the high-energy nuclear collision processes and gave a prediction for the energy dependence of the multiplicity, which was verified by the data, it was known that it had troubles in reproducing particle spectra and relative abundance of K over p. This was because, in this model, particles were assumed to be emitted directly from the hot and dense, thermally equilibrated matter formed in high-energy nuclear collisions, which was supposed to be at rest, so that the model predicted isotropic momentum distribution which did not agree with the observed spectra. Furthermore, because of high temperature (T >> m_K) reached in the process, multiplicity ratio depended only on the isotopic-spin statistical weights, namely K/p = 4/3. This conclusion of the model was also not in agreement with data.

These problems were solved naturally by letting the hot and dense matter expand before particle emission takes place, reducing thus heavy-particle multiplicities, because of the Boltzmann factor, and giving at the same time alongated momentum spectra, due to a violent longitudinal expansion caused by a large pressure gradient in the beam direction. A nice feature of this model is that, since the entropy is conserved in the ideal case Landau studied, the energy dependence of the total particle multiplicity predicted by the Fermi model, and verified experimentally, is preserved.

When accelerator data on multiparticle production began to appear, first in pp collisions at CERN ISR, and later in

According to Hydrodynamic Model, the description of high-energy nuclear collisions goes as follows. At the beginning, two Lorentz contracted (in the c.m. frame) nuclei collide and it is assumed that, after a complex process involving microscopic collisions of nuclear constituents, a hot and dense matter is formed, which would be in local thermal equilibrium. The description of this initial thermalization process is out of the scope of hydrodynamics. In hydrodynamics, we simply assume that the local thermal equilibrium is attained and these states of matter are specified by some appropriate initial conditions (IC) in terms of distributions of fluid velocity and thermodynamical quantities for a given time-like parameter. Then, it follows a hydrodynamical expansion, described by the conservation equations of energy-momentum, baryon number and other conserved numbers, such as strangeness, isotopic spin, etc.

where

is the energy-momentum tensor, n_B, n_S, e, p are, respectively, the baryon number density, the strangeness density, the energy density and the pressure, all of them given in the proper frame of reference of the fluid element, and u^m is the four-velocity of the fluid. Moreover, we have to specify some equations of state (EoS), which depend on the nature of the hot matter produced.

As the expansion proceedes, the fluid becomes cooler and cooler and more rarefied, occurring finally the decoupling of the constituent particles, that is, they don't interact any more until their detection. However, long-lived resonances and other unstable particles may decay after this instant of time. The observable quantities such as dN/dy, ds/dm_T, < v₂ > , ... are then computed by using these decoupled or free particles.

The main object of studies by using Hydrodynamic Model is to investigate, through comparison of its predictions with data, properties of the matter formed during high-energy nuclear collisions, specified by the initial conditions, equations of state and freeze-out or decoupling conditions. We emphasize that these properties are not known a priori. It should also be stressed that even the basic assumption of ''local equilibrium'' is not granted for a priori. We expect that experimental and theoretical studies of some appropriate observables may respond these questions. Therefore, it is fundamental to find what are these ''most appropriate observables''.

In this survey, we shall discuss some aspects of this model, by focusing mostly on those ones, like development of hydrodynamic code capable to treat problems with highly asymmetrical configurations, effects of the initial-condition fluctuaions and improvement of the description of decoupling process. These are features which have been investigated and developped within the São Paulo - Rio de Janeiro Collaboration in the last ~ 15 years. For a review of other aspects of recent developments, see for instance Ref. [4,5,6].

In the following, in the next Section, we discuss the initial conditions, by emphasizing the importance of the event-by-event fluctuations as shown by some event generators. In Section, we describe several equations of state, usually employed in these studies. Section is devoted to the resolution of hydrodynamic equations. There, we begin describing some analytic solutions, turn to the variational formulation and, finally, an application of this approach to develop a numerical code, using algorithm of smoothed-particle hydrodynamics. Then, we consider the decoupling mechanisms in Section, by stressing that, although the commonly used Cooper-Frye presciption [7] is convenient and can give many good results, more realistic treatment of decoupling is needed in order to correctly extract information on the hot matter formed in the collision process. In Section, we give some results obtained with the methodology described here. Finally, conclusions are drawn in Section 7.

2 Initial conditions

In usual hydrodynamic approach of high-energy nuclear collisions, one customarily assumes some highly symmetric and smooth IC, parametrized in a convenient way, which would correspond to the mean distributions of hydrodynamic variables averaged over several events [5,6,8]. However, our systems are not large, so large fluctuations varying from event to event are expected, even under the same initial conditions of colliding objects, such as the incident energy and the impact parameter of the nuclei. What are the effects of the event-by-event fluctuation of IC? Are they sizable? Do they depend on EoS? These are some questions which arise regarding this subject.

As mentioned in the Introduction, IC are determined by a complex process involving microscopic collisions of nuclear constituents not accounted for by hydrodynamic model, so when we want to introduce fluctuations in the IC of a hydrodynamic system, we must go beyond the hydrodynamic degrees of freedom. Just to see whether such event-by-event fluctuations of IC give sizeable effects, so merit a more detailed study, in [9], Paiva et al. used the Interacting Gluon Model [10] (IGM) to generating fluctuating IC and, using Khalatnikov 1-dimensional solution [11], showed that the rapidity distribution obtained by averaging over results starting from fluctuating IC is quite different from that obtained starting from the averaged IC.

There are some other simulations, which try to incorporate, in hydrodynamic computations, fluctuating IC given by more elaborate microscopic models: with a use of some event generator, e.g. HIJING [12], VNI [13], URASiMA [14], NeXuS [15], or some effective theory such as string model [16], perturbative QCD + saturation of produced partons [17] or color glass condensate [18]. In principle, one could test each of these different microscopic models, by connecting them to some hydrodynamic code and computing several observables to see which are the differences among them and which are more suitable for describing experimental data, provided the other ingredients of the hydrodynamic model are well known, that is not the case. Here, we shall instead discuss not the details of such models, but more or less model-independent consequences of such fluctuations. Anyhow, we have to adopt some microscopic model. In the following, we shall mainly discuss the recent works of São Paulo-Rio de Janeiro Collaboration, using NeXuS event generator, coupled to hydrodynamic code SPheRIO¹ 1 Smoothed Particle hydrodynamic evolution of Relativistic heavy IOn collisions. .

NeXuS is a microscopic model based on the Regge-Gribov theory and can give, in the event-by-event basis, detailed space distributions of energy-momentum tensor, baryon-number, strangeness and charge densities, at a given initial time t = ~ 1fm, for any given pair of incident nuclei or hadrons. One important point when we use a microscopic model to create a set of IC for hydrodynamics is that the energy-momentum tensor produced by the microscopic model does not necessarily correspond to that of local equilibrium. For example, NeXuS generates, as its output, the energy-momentum tensor T^mn(x) and the current densities of conserved quantum numbers, (x), (x) and (x), where B, S and Q refer to baryon number, strangeness and electric charge, respectively. However, the four-velocities corresponding to these currents usually do not coincide, and more importantly, do not coincide with that of the frame where T^mn becomes diagonal. Furthermore, the space components of the diagonalized T^mn are not necessarily identical (anisotropic stress). These facts mean that the matter is not in local equilibrium. In order to transform the energy-momentum tensor to that of the equilibrated one, we adopt the following procedure. First, following Landau [19], we identify the normalized time-like eigenvector of T^mn as the four-velocity of the fluid and the eigenvalue as the energy density,

Using this four-velocity, we calculate the proper baryon number density as

and analogously for the other densities. Once e and n's are obtained, all the other thermodynamical quantities are calculated using the equations of state. By this procedure we force the system into a local thermal equilibrium, conserving the proper energy density of the system.

The IC thus generated are used as inputs for SPheRIO code. We show in Figs.1 and ² an example of such a fluctuating event, produced by NeXuS event generator, for central Au+Au collision at 130A GeV, compared with an average over 30 events. As can be seen, the energy-density distribution for a single event (left), at the mid-rapidity plane, presents several blobs of high-density matter, whereas in the averaged IC (right) the distribution is smoothed out, even though the number of events is only 30. Similar bumpy event structure was also shown in calculations with HIJING [12]. So, the main question here is whether the averages over the observables computed starting from such fluctuating IC, like the one at the left panel of Figs.1 and ² are similar or sizeably different from the correspondent ones computed from averaged smooth IC like that at the right panel there, or symbolically,

where < f > is the value of some relevant quantity f, obtained by averaging over N total number of events and (IC® f)_j is the value of the same quantity in the j-th event, with some event-dependent IC, whereas ( < IC > ® f) represents the same quantity f given with the average IC. This is a crucial point in data analyses, because the left-hand side is closer to the data point experimentalists obtain, whereas the right-hand side is the quantity usually computed by theorists for that data point, in order to extract properties of the matter formed in nuclear collisions. If the bumpy structure shown by event generators effectively exists in experimental situations, then how do these hot spots manifest themselves in the observables?

A general conclusion one can draw about this question is that the total entropy of the system becomes always smaller when one takes such fluctuations into account, in comparison to the case without fluctuations, which means with average over the event-by-event fluctuating IC taken before the expansion. This can be seen by observing that, in ideal hydrodynamics, both energy and entropy are conserved. Then, considering for simplicity an ideal gas so that S_i = a(E_i)^3/4, with a = const. > 0, for each random event,

Here, < E > and < S > in the left-hand sides mean the averaged energy and entropy over the fluctuating events, whereas the right-hand side of < S > is the entropy corresponding to the averaged initial conditions, with the averaged energy < E > . The linear terms of the expansion in DE_i/ < E > are cancelled out when the summation is performed. If one recalls that particle multiplicity is proportional to the entropy for each particle species, one would expect that also the multiplicity becomes, in general, smaller when one takes such fluctuations into account.

Other possible manifestations of this inhomogeneity of IC that we can expect intuitively are: enhancement of high-p_T components due to more violent expansion in the surface region, smaller HBT radii, due to concentrations of matter in small spots, azimuthal asymmetry even in central collisions. However, computations are needed to obtain quantitative conclusions, whether such discrepancies are meaningful or not. We shall discuss this question in Sec.VI.

3 Equations of state

As mentioned already, the basic assumption in hydrodynamical models is the local thermal equilibrium. Once this condition is satisfied, all the thermodynamical relations should be valid in each space-time point² 2 Recently it has been suggested that the thermodynamical relations can be satisfied without having the thermal equilibrium in the sense of Boltzmann distribution[20,21]. . Thus, the energy, pressure and temperature are given as functions of baryon number and entropy densities, specifying the properties of the matter. In this Section, we discuss how to obtain simple phenomelogical equations of state (EoS) for the hydrodynamical description of relativistic nuclear collisions [22].

3.1 Hadronic gas

The strong interactions among hadrons are very complicated and difficult to be incorporated into the EoS for practical use. However, for very high energy, we may consider that the hadronic gas may be approximated as an ideal gas, although the degree of approximation can not be evaluated theoretially. The recent thermal model for the description of chemical abundances [23] show that such an approach can reproduce quite well the observed multiplicity ratios of produced hadrons. Here we assume that all the particles can be treated as constituents of a quantum ideal gas, except for a correction due to the excluded volume. We also include a main part of observed resonances in Particle Data Tables. The inclusion of resonances can be considered as an effective way to consider the interactions among hadrons as explained later.

First, we recall that, in a grand canonical ensemble for an ideal gas of quantum particles, the thermodynamical potential per volume (the pressure) is given by

were q = ±1 (+ for fermions, - for bosons), b = 1/T is the inverse of the temperature T, m the chemical potential, g the degeneracy factor and e(k) = with m the mass of the particle. The number density n and the energy density e can be obtained by the usual thermodynamical relations, n = (¶p/¶m)_V,T, e = (¶p/¶b)_l, where l = e^bm is the fugacity. The entropy density of the gas can be calculated as s = b(p+e-mn).

For example, in Landau's model [1], the equation of state was taken as that of the massless pion gas. For bosons with m = m = 0, Eq.(10) can be integrated analytically to give

and accordingly

For the pion gas, due to the isospin factor, we can take g = 3.

For more realistic equations of state, we should include all the resonance particles in the gas. Furthermore, we should also take into account more than one type of conserved quantum numbers, such as electric charge (equivalently the 3^rd component of isospin), baryon number and strangeness. In this case, the chemical potential must be written as m = Bm_B+Sm_S+T⁽³⁾m₃ where B,S,T⁽³⁾ are baryon, strangeness and the thrid component of isospin quantum numbers, respectively, and m_B, m_S and m₃ are the corresponding chemical potentials. Thus, for a mixture of particles with these conserved quantum numbers, Eq.(10) should be generalized to

where the sum refers to the particle species (including resonances) and m_i = B_im_B+S_im_S+T_i⁽³⁾m₃, with B_i,S_i,T_i⁽³⁾ the quantum numbers of the i-th particle species. We verify that the baryon number density of the mixture is

where n⁽ⁱ⁾ = (¶p_i(T,m_i)/¶m_i) is the number density of the i-th particle species.

Except for pions, most of hadrons and resonances can be well approximated by the Boltzmann limit. In this case, we have

and

where K₂ is a modified Bessel function. From these relations, we can see immediately the usual ideal gas relation,

When the widths of the resonances are taken into account, Eq.(13) must be modified. For an interacting gas, the power series expansion of the pressure in terms of fugacity,

is known as the cluster expansion (which is intimately related with the virial expansion) and b_n are called ''virial'' coefficients. Here, p^id is the pressure of the corresponding ideal gas and, roughly speaking, the index n in the sum represents the order of multiple particle interactions. For n = 2, the contribution to the pressure comes from the 2-body interactions. Beth and Uhlenbeck [24] showed that the second virial coefficient can be expressed in terms of the scattering phase-shift of constituting particles. This approach was generalized to the relativistic Boltzmann gas by Dashen, Ma and Bernstein [25] and the result for b₂ is

where d_l(W) is the phase shift for the l-th partial wave.

Consider the case of a gas of particles with mass M and suppose there exists a resonance in the two particle collision,

When the resonance has a width G and spin J, only l = J dominates the sum and

so that we have the Breit-Wigner formula,

Therefore, the pressure of the system can be written as

p = p^id+p_R,

where

with

g_R = 2S+1,

m_R = 2m.

For extremely narrow resonances, G®0,

which is exactly the pressure of the ideal relativistic Boltzmann gas made of resonances with mass M_R. Eq.(23) suggests that, for more general case, the effect of resonance width can be obtained by a convolution of the normalized mass spectrum f(M) of the resonance, with the pressure p^id(M) of the ideal gas of mass M, as

p_R = dM_Rf(M_R)p^id(M).

3.2 Effect of resonance width as a function of temperature

Equation (23) shows that the effect of the resonance width on the pressure of the gas is temperature dependent. To see this effect, let us introduce the quantity,

so that

Another way to see the effect of the width, we may introduce the effective mass of the resonance M_eff defined

by

Using this effective mass, we can write the resonance pressure as

that is, as if it were the pressure of ideal gas of particles with mass M_eff. In Figs.3, ⁴ , ⁵ and ⁶ , we show the temperature dependence of F and M_eff for two tipical cases, one for the resonance r (light and narrow width) and the resonance f (heavy and large width) resonances. As we see, the ideal gas approximation deviates substancially for low temperatures, especially for large width resonances. The ideal gas approximation is only valid for T » G.

3.3 Excluded-volume correction

From the analysis of thermal models [23], it became clear that the ideal gas description requires a modification to adjust the size of the system. The volume to fit the particle abundances is found to be too small. To avoid this problem, the correction due to the excluded volume effect, like a Van der Waals hard core correction is introduced [26]. According to this prescription, Eq.(13) is modified to the following coupled equations

where as before m_i = B_im_B+S_im_S+m₃ is the chemical potential and v_i is the excluded volume of the i-th hadron species. The superscript id refers to the ideal gas case. The above equations constitute an implicit equation for p_HG so that these two equations are solved iteratively to obtain p_HG for a given set of parameters, T,m_B,m_S and m₃. The number density of the i-th hadron is given by

where is the ideal gas expression of the particle density for the i-th particle species.

3.4 Gas of Quarks and Gluons

The simplest way to introduce the phase of quarks and gluons in the equations of state is the use of the MIT Bag models. The effect of gluon and quark condensate in the physical vacuum is expressed as the energy density of the vacuum (or vacuum pressure). Thus the energy density and pressure of an ideal quark-gluon gas calculated in the QCD vacuum should be modified according to the rule,

e ® e+B,

p ® p-B,

where B is the vacuum pressure. Note that this vacuum pressure has a property analogous to the cosmological constant L of Einstein. Now, when we consider just the u and d quarks and neglect their masses, we have

with

g_q= 2×2×3,

g_G = 2×8,

the statisfical factors of quarks and gluons. For quarks, we have m_q = m_B/3. For m_B = 0, we have

or effectively g_qgp = 37. To include the strangeness and also charge conservation, we proceed in the same way as for the hadronic gas and we have

where g_l = 2×3 and

3.5 Construction of equations of state for the practical use

The expressions Eqs.(27,28) or Eq.(3.4) are, however, not convenient for the use in hydrodynamical calculations. This is because the variables in such calculations are the conserved quantum numbers and the entropy density and not the chemical potentials and the temperature. So, we need to invert the relations,

to get

However, this is a formidable task even numerically. We are thus forced to reduce the degrees of freedom for the practical application to hydrodynamics. For this purpose, we set the isospin and strangeness densities to null everywhere. That is, we impose the conditions,

These conditions together with Eqs.(33) and (34) determine m_S and m₃ as functions of T and m_B. Therefore, n_B and s in Eqs.(32) and (35) become now functions of two variables T and m_B,

which can be inverted numerically and give

The above inversion process allows us to write any thermodynamic quantity as function of n_B and s, both in hadronic gas and quark-gluon plasma. Then, for a given pair of density parameters n_B and s, we should determine which phase the physical system assumes. When two phases are in equlibrium, we must have [27]

so that it determines the phase boundary line in the (m_B,T) plane and separates this plane into two domains. The domain where p_HG > p_QGP is the hadron gas phase and the other p_HG < p_QGP is the QGP phase. These two domains are in contact on the phase boundary line Eq.(48).

However, the above two domains in the (m_B,T) plane are mapped into two separate domains in the (n_B,s) plane and there appeas a new third domain between them. That is, the phase boundary line in the (m_B,T) plane spreads into a domain in the (n_B,s) plane. This domain is the mixed phase. In oder to determine thermodynamical quantities in this mixed phase as functions of (n_B,s), we should introduce another criterion in addition to the phase boundary condition of Gibbs. As mentioned above, any point (m_B,T) on the phase boundary line corresponds to two points in the (n_B,s) plane: one for the hadron phase, (, s^H), and the other, (, s^Q), for the QGP phase. In the mixed phase, any density of extensive quantity should be a linear function of the qgp/hadron concentration ratio. Thus for a given value of the baryon number density n_B in the mixed phase, (, s^H) and (, s^Q) should satisfy

From this equation, we determine the two points on the phase boundaries, (, s^H) and (, s^Q). All the other extensive quantities, say e, can then be obtained as

Finally we can construct the equations of state in the whole region of (n_B,s). The parameters of the final equations of state are then,

• Number of resonances included in the hadronic gas: Here, we take all the mesons with mass smaller than 1.5 GeV, and baryons smaller than 2 GeV. Resonance widths are not included.

• Quark masses: We may safely take m_u = m_d = 0, but for the strange quark, we take m_s = 120 MeV.

• Size of excluded volume: In the example shown in the figures below, n₀ = (4pr/3) with r₀ = 0.7 fm for baryons and r₀ = 0 for mesons.

• Bag constant: We take B = 380 MeV/fm³.

Figure 7 shows the line of constant temperature in (n_B, s) plane. Fig.8 shows the lines of constant entropy per baryon in (m_B,T) plane.

Note that the present equations of state have, by construction, the first order phase transition between hadronic gas and quark gluon plasma. Recent lattice calculations [28] indicate that there exists a critical point so that the phase transition in the region starting from that point to zero m_B axis is not the first order but may be either a higher order transition, or a smooth cross over [29]. Effects of such possibilities should be investigated within the hydrodynamical models.

4 Resolution of hydrodynamic equations

In general, exact analytic resolution of relativistic hydrodynamics is a difficult task due to the highly non-linear nature of these equations. So, usually one resorts to numerical computations. However, since the analytical studies are more transparent, it would be useful to find analytical soluions, even though they correspond to highly ideal cases. Khalatnikov's one-dimensional analytical solution [11] to Landau's initial conditions [1], namely, ideal fluid at rest in a Lorentz contracted thin spatial region, gave rise to a new approach in high-energy physics. The boost-invariant solution [30] found 20 years later, is frequently utilized as the basis for estimations of initial energy densities in ultra-relativistic nucleus-nucleus collisions [31]. Below, we shall first describe families of analytical solutions we obtained in collaboration with T. Csörgo [32,33].

Considering that it is not trivial to analytically solve the equations of hydrodynamics, it would be nice if we could develop a method to obtain an approximate but analytical solution of hydrodynamics. This has been done with variational formulation [34], although we had not develop it further by applying it to practical problems of high-energy nuclear collisions. However, we did apply the variational method to develop a numerical code SPheRIO, based on the so-called smoothed-particle hydrodynamics (SPH) [35,36] tecnique, to high-energy nuclear collisions [37], which is flexible enough, capable to treat systems with configurations without any symmetry and also exploding in time.

In the following, after presenting some analytical solutions, we shall describe in the subsection the variational formulation of hydrodynamics, showing how it could used to get approximate solutions. Then, in the subsection , we shall apply this method to adapt the SPH hydrodynamics for relativistic heavy-ion collisions.

4.1 Analytic solutions

After Landau's initial proposal [1], the first analytic solution obtained is due to Khalatnikov [11], considering 1-dimensional expansion of ideal gas. A more simpler boost-invariant solution has been obtained later [30] and applied to estimate the initial energy densities in ultra-relativistic nucleus-nucleus collisions [31]. Both of these have been frequently used in the study of nuclear collisions, showing the usefulness of such simple analytical solutions.

However, the boost-invariant solution has some shortcomings: i) it is scale invariant, having a flat rapidity distribution, corresponding to the extreme relativistic collisions, which has never been seen; ii) it contains no transverse flow. In [32], we tried to overcome these shortcomings. We started by assuming a simple equations of state, corresponding to a gas containing massive conserved quanta, namely,

having two free parameters, m and k. Non-relativistic hydrodynamics of ideal gases corresponds to the limiting case of m >> T, v² << 1 and k = 3/2.

Then, we looked for similarity flows, i.e.,

where X(t), Y(t), Z(t) are the scales of length, in three orthogonal directions, and f(x^m) is any of the thermodynamical quantities, such as n(x^m), T(x^m), p(x^m), ..., and

is a scaling variable.

Using this parametrization for the one-dimensional expansion, it can be easily verified, by direct substitution into eqs.(1,2) with eqs.(4,51,52,53 and 55), that a family of solutions can be written as

with p₀ = n₀T₀ and where V(s) is an arbitrary non-negative function of s = z²/(

Next, considering cylindrically symmetric flows, with boost invariance along z direction (collision axis), we could find the following family of solutions, for transverse flows:

Here, p₀ = n₀T₀ and V(s) is an arbitrary non-negative function of s ( = r/(t) in this case), with r_t = and ₀ = and t_z = . So, this is a generalization of the one-dimensional scale-invariant solution, including a class of transverse flows.

More recently [33], we extended these solutions still further, considering less symmetrical flows, but still keeping the same EoS (1,2) and similarity flow, with constant velocity, as it appears in eqs.(56,57).

4.2 Variational formulation

As shown above, even for very simple equations of state, analytic solutions are limited to special cases. For realistic situations, even the equations of state themselves are available only in the form of numerical tables. Therefore, the numerical resources are essential for realistic studies of hydrodynamical behavior of ultra-relativistic collisional processes. However, it is well-known that any numerical method for partial differencial equations requires highly sophisticated techniques to avoid numerical instabilities, and usually it needs a very large scale computation, especially when we want to describe correctly explosive processes such as relativistic heavy-ion collisions. However, as emphasized in the Introduction, in the hydrodynamic approach of high-energy nuclear collisions, its main ingredients, i.e., the equations of state, the initial conditions and the freezeout conditions are not quite well known. In such a situation, we actually don't need a very precise solution of the hydrodynamic equations, but a general flow pattern which characterizes the final configuration of the system as a response to given equations of state, initial conditions and the decoupling procedure. So, we prefer a rather simple scheme of solving the hydrodynamic equations, not unnecessarily too precise but robust enough to deal with any kind of geometry. From this point of view, we stressed in Ref. [34] the advantage of a variational approach to relativistic hydrodynamics.

Although not commonly found in general textbooks, the variational formulation of hydrodynamics has been studied by several authors [38,39]. In Ref. [34], starting from the action

where e is the proper energy density, the relativistic hydrodynamics was derived from the variational principle.

Here, we show its derivation, generalizing to include the rotational flow. To do this, we introduce the two variables, the proper baryon density, n, and entropy density, s, which satisfy the conservation laws,

where u^m is the four velocity of the fluid, satisfying the normalization,

The functional form of the energy density,

e = e( n, s ),

specifies the thermodynamical properties of the fluid. The pressure, temperature and chemical potential are obtained by the usual thermodynamic relations,

The hydrodynamical equations of motion for the fluid is given by the variational principle with respect to n, s and u^m under the constraints, Eq.(59), and the normalization of the four-velocity, Eq.(60). These constraints can be incorporated in the variational principle in terms of Lagrangian multipliers to write

where l,z and w are Lagrangian multipliers and arbitrary functions of x. Equivalently, the fluid dynamics is given by the effective Lagrangian,

where now all of n,s,u^m,l,z,w are independent variational variables.

The variations with respect to n,s and u^m lead immediately to

Variations with respect to l,z and w give simply the constraints, Eqs. (59) and (60). Multiplying the both sides of Eq.(65) by u^m, and using Eqs. (60,63,64), we get

where p is the pressure. Eq.(66) shows that w is the enthalpy density. Substituting back this w into Eq.(65) and multiplying by u_n, we have

wu_mu_n = -(nu_n)(¶_ml) -(su_n)(¶_mz).

Taking the divergence and using the continuity equations ¶ⁿ( nu_n) = 0 and ¶ⁿ(su_n) = 0, we get

But

and analogously

so that Eq.(67) becomes

Here, we have used Eq.(65) and the Gibbs-Duhem relation,

and the property,

uⁿ(¶_mu_n) = 0.

Finally we arrive at the standard form of the relativistic hydrodynamic equation (1),

where

is the usual energy-momentum tensor of the fluid.

It is important to observe that, the effective Lagrangian Eq.(62) evaluated in the proper comoving frame of the fluid motion is

which is nothing but the negative of the thermodynamic potential of the fluid element at x^m.

Now, interesting application of this approach appears, if we can parametrize possible solutions of continuity equations (59) in terms of a certain number of time-dependent parameters (t) = {a_i(t),i = 1,...,N}, such that

together with the velocity field,

then the action, Eq. (58), may be written as a time integral of an effective Lagrangian

The constraint terms vanish for these ansatz. In this case, the equations of motion for the variables a_i(t) are obtained as the Euler-Lagrange equations. This method could be applied to relativistic heavy-ion collisions, trying to describe them in a simple analytic and effective way in terms of few parameters. However, the general parametrization which solve exactly the continuity equations is not easy. An approximate way to solve the continuity equation is proposed in a numerical method called Smoothed Particle Hydrodynamics.

4.3 Smoothed particle hydrodynamics

The SPH algorithm was first introduced for astrophysical applications [35,36]. In [37], we extended this numerical method to heavy-ion collisions by the use of the variational approach discussed in the preceding subsection.

4.3.1 SPH representation of densities

The basic idea of the SPH method is to parametrize the continous density distribution of any extensive physical quantity in terms of sum of base functions with finite support. This procedure introduces two types of approximations of different nature. To see this, let us suppose that A is the physical extensive quantity and a(r,t) the corresponding density distribution. We start with the identity,

a(r,t) = a(r',t)d(r-r') d³r'.

Now, let us introduce the first approximation. Substitute the Dirac d-function by a smooth, normalized function W with finite support, say, h, and transform the density a to ã as

where as mentioned, W is normalized,

and having the property of finite support,

W(r-r';h)® 0,for|r-r'| > h.

At this stage, the new density ã(r,t) describes the smoothed part of the original density a(r,t). From the Fourier transform, we can see that for this smoothed density, the Fourier components with large wave numbers, corresponding to

vanish rapidly. In other words, the kernel function W serves as the short wavelength cut-off filter. Physically, it is useful to introduce such a filter, since we very often want to eliminate very short scale part in order to extract the global feature of the dynamics of the system. In 60's, similar idea has been used to smooth out the spectrum density of the nuclear shell model to extract the collective liquid drop potential by Strutinski [40].

Now, we introduce the second approximation. This is rather to do with the practical aspect, that is, to reduce the degrees of freedom involved in the calculation. We replace the integral Eq.(74) by the sum over a finite and discrete set of points,{r_i,i = 1,..,N}.

where the weight A_i should be chosen appropriately to minimize the difference between ã(r,t) and a_SPH(r,t) everywhere. The above expression means that we are representing the continuous density as due to a sum of a finite number of dynamic units carrying the quantity A_i. These units are centered at the position r_i.

Finally, the correspondence,

can be considered as an approximation ansatz for the density a(r,t) with a finite number of parameters, {A_i,r_i ,i = 1..,N}. Due to the normalization of the kernel W, we have

so that we should choose

as the total value of the quantity A of the system.

4.3.2 Solution of continuity equations

For the application in Hydrodynamics, we can use these parameters as the variational variables so that they are time dependent. When we deal with one or more extensive quantities, we usually choose one conserved quantity as the reference density, say r, and represent it by the SPH form, choosing appropriately the weights {n_i} to get

and take n_i constant in time. Other extensive quantity, say A, is calculated as in Eq.(76) with weights,

The quantity (a/r)_i then represents the quantity A for the unit reference quantity r at the position r = r_i(t). Note that the time dependence of the density in Eq.(77) comes from those of {r_i(t)}. In this sense, Eq.(77) can be understood as if {r_i(t)} were Lagrangian coordinates attached to small volumes of the order of h³, with some conserved quantity, such as baryon number or entropy in the adiabatic expansion. From now on, we refer to these dynamic units characterized by their coordinates {r_i} and carrying the quantity A_i as ''SPH particles''. From Eq.(78), the quantity A_i can be interpreted as the part of A carried by the i-th SPH particle.

The most powerful point of the above scheme of SPH representation is that we can solve the continuity equation in a very simple manner. Suppose M is a conserved quantity. Then the corresponding density r should satisfy the continuity equation,

where v is the velocity field. The SPH expression for the current j = rv is

so that

On the other hand, from Eq.(77),

By inspection, if we identify

then we can see that Eq.(79) is automatically satisfied.

For the application to the relativistic heavy ion collisions, we can take the entropy and baryon number as the basic conserved quantities. Then, their densities (in the space-fixed frame) are parametrized as

where n_i and b_i are the entropy and baryon number attached to the i-th ''particle''. The total entropy and baryon number are then given by

The proper densities of entropy and baryon number are related with these space-fixed frame quantities as

s = g-1s*,

n= g-1n*,

where g = u⁰ is the Lorentz factor asscociated with the fluid velocity. For the hydrodynamical description of nuclear and hadronic collisions at ultra-relativistic energies, we prefer to use the entropy than the baryon number as the reference conserved number to write the SPH representation of any other extensive quantities. This is because, the baryon density may become zero but the entropy density never vanishes in the physically interesting region within a hydrodynamical description.

4.3.3 SPH action and SPH equations

In the variational derivation of this method, the set of time-dependent variables {r_i,i = 1,...,n} are taken as the variational degrees of freedom and their equations of motion are determined by minimizing the action for the hydrodynamic system. Thus, SPH may be considered as an effective description, in which the coordinates {r_i(t)} associated with particlesare the optimal dynamical parameters which minimize the model action. We observe that {n_i} or {b_i} are not dynamical variables and are determined by the initital conditions together with the constraints for the variational procedure.

The effective Lagrangian, Eq.(71), is rewritten in SPH representation as

where E_i is the ''rest energy'' of the i-th ''particle''. Then, the equations of motion are obtained from the usual variational procedure. This leads to the following coupled equations

4.3.4 General coordinate system

The variational procedure can readily be extended to coordinate system, with a non-Cartesian metric. The use of generalized coordinate systems is particularly important when we consider realistic initial conditions for simulations of RHIC processes. As we know, in a relativistic heavy-ion collisional process, the initial state is a cold, quantum nuclear matter. Just after the collision, the hadronic matter stays at a highly off-shell state and the materialization occurs only after ~ 1fm/c in the proper time. Therefore, the local thermodynamical state would emerge for some local proper time and not for the global space-fixed time t. Thus, it is important to choose a convenient coordinate system for the description of the relativistic heavy-ion collisions. For example, one often uses the hyperbolic time and longitudinal coordinates to be described later.

Let us consider a general coordinate system, with

However, in order to unambiguously define the conserved quantity, we consider only the case where the time-like coordinate is orthogonal to the space-like coordinates,

The action principle for the relativistic fluid motion can be written as [34]

together with the constraint for the conserved entropy current,

or

where

and we use the notation

t = x⁰,g = u⁰.

The generalized gamma factor g is related to the velocity

where -g is the 3×3 space part of the metric tensor. That is

Let us now introduce the SPH representation. We may, for example, express the entropy density by the ansatz

or by

as well. These two possibilities, besides others, are simply different ways to parametrize the variational ansatz in terms of a linear combination of given functions W(-_a(t)). The most important property of an ansatz should be that W satisfies the normalization condition imposed by the basic conserved quantity. Since the total entropy is expressed as

the normalization of W should be taken to be

for the parametrization Eq.(94) and

for the parametrization Eq.(95). In the usual SPH calculations, it is not desirable to introduce in W the space-time dependence through its normalization condition. In this respect, the most natural way to introduce the SPH representation is Eq.(94). With this choice, the SPH action is given by

The variational principle leads to the following equation of motion,

where

and the operator Ñ is just the simple derivative operator with respect to the coordinate variables in use.

For ultrarelativistic heavy-ion collisions, a useful set of variables is

As mentioned above, the initial conditions for RHIC processes are specified in terms of the proper time rather than of the fixed-frame time t. The variable t is not exactly the physical proper time of the matter, but in general it is a good approximation in ultra relativistic collisions.

The metric tensor for this coordinate system is given by

Since the metric is space independent, we can use the parametrization

where

and W is normalized as

The SPH equation becomes

where the h component of the momentum is related to the velocity dh/dt as

whereas in the transverse direction, we have

The Lorentz factor is given by

4.3.5 Landau Model

In order to show the efficiency of the method and also to show the correct choice of the coordinate system, let us investigate the Landau model in the SPH scheme, using the ordinary Cartesian coodinates and h-t coodinates. Since the analytical solution is known, we can compare the numerical solutions to it.

We thus solve the hydrodynamical evolution of a system of one-dimensional relativistic massless baryon-free gas initially at rest. The equation of state of a relativistic massless boson gas is

where

To apply the SPH method, we introduce the discrete one-dimensional space variable x_i(t),i = 1,..,n (and similarly for h-t coodinates). The relation between the momentum and velocity is then

where, in this case, v can be solved analytically with respect to p.

In Figs.9a and ^b, we show the results of our SPH calculation together with the exact solution [11,1]. In these examples, we took only 100 particles with equally spaced x_i (or h_i). As we see from this example, in spite of rather small number of particles, the SPH solution is quite satisfactory. In particular, when we use the h-t coordinates with an appropriate distribution of s (^Fig.9b), an excellent agreement with the analytical solution can be obtained. The computation time needed to get these solutions is even less than that needed to numerically evaluate the analytical solution.

4.3.6 Transverse expansion on longitudinal scaling expansion

As a further test, closer to a realistic situation than that of Figs.9, we calculated the transverse expansion of a cylindrically symmetric homogeneous massless pion gas, undergoing a longitudinal scaling expansion, and initially at rest in transverse directions. Such a problem has been discussed by several authors as a useful base to understand the transverse expansion. In Fig.10, we compare our results (a full 3D calculation without assuming cylindrical symmetry) with (2+1) numerical results, obtained by the use of the method of characteristics [41]. In this example, we used also 50×50×50 particles. The result is quite satisfactory. If we decrease the accuracy by 10%, we can reduce the particle number almost by one order of magnitude.

4.3.7 Shock formation and Neumann-Richtmyer pseudo viscosity

As seen in the previous examples, our entropy-based relativistic SPH method works quite well for the adiabatic dynamics of the massless pion gas. However, for the application to realistic problems, it is fundamental to see how this scheme works for non-adiabatic cases, too. This is because, whenever a piece of fluid matter flows into another region of the fluid with a speed exceeding the sound velocity of the fluid, there appears a shock wave, and this is essencially a nonadiabatic process. Thus, except for a really quasi-static dynamics, there should be an entropy production mechanism. This becomes especially important in a domain close to the phase transition region, because there the velocity of sound tends to zero. In the following, we study some examples of one dimensional shock problems in the scheme of the SPH methods.

The shock front manifests as a discontinuity in thermodynamical quantities in a hydrodynamic solution. Mathematically speaking, the shock front should be treated as a boundary connecting two distinct hydrodynamic solutions. The smoothed particle ansatz excludes such a possibility from the beginning. Since short-wavelength excitation modes do not exist in the SPH ansatz, the energy and momentum conservation required by the hydrodynamics results in very rapidly oscillating motion of each SPH particle. Such a situation occurs, for example, when a very high-energy density gas is released into a low density region. This kind of shock, for the case of a baryon gas, is discussed in [42] and also, in the SPH context, in [43]. Here, we appliy our entropy-based SPH approach to the massless pion gas.

Figure 11 gives the typical behavior of SPH solution for such a situation, if entropy production is not taken into account. As discussed above, there appear in fact rapid oscillations in thermodynamical quantities just behind the shock front. Actually, such oscillations always appear in any numerical approach if entropy production is not included.

In order to avoid these unphysical oscillations, von Neuman and Richtmyer [44] introduced the concept of pseudoviscosity. The idea is to set the dissipative pressure where the shock wave discontinuity is present. To do this, Neuman and Richtmyer proposed to replace the pressure by

p® p+Q,

where Q is the pseudoviscosity and they took the following ansatz,

The above formula is for nonrelativistic one-dimensional hydrodynamics. Here, r is the mass density, Dx is the space grid size and a is a constant of the order of unity. In order to generalize the above pseudoviscosity for relativistic SPH case, we replace the quantity /r by -q = -¶_mu^m and Dx by h, where h is as before the width of the smoothing kernel W. More precisely, we take the following form which is a slightly modified expression suggested by Ref. [43],

where

Actually, Q is equivalent to the bulk viscosity and therefore there is no heat flow associated with it. What this artificial viscosity does is to convert the collective flow energy into the microscopic thermal energy. As a consequence, the total energy, that is, the sum of the collective flow energy and the internal thermal energy is still conserved. In order to incorporate the internal energy conservation in the SPH scheme, we substitute all the pressure p_i by p_i+Q_i, and we add the following equation for the entropy production,

Figure 12 is the solution of the same problem as in Fig.11, but with the entropy production taken into account. In this calculation, the parameters have been chosen as

a = 2,b = 4

and h = 0.5fm for 1000 SPH particles. As we see, the rapid oscillations have been smoothed out (and in turn, the numerical calculation became much more efficient).

It is known that the energy- and momentum-flux conservations through a shock front relate the ratio s₂/s₁ of entropy densities after and before the shock to the velocity v_s of the shock front as (Hugoniot-Rankine relation)

In Fig.13, we show the velocity of the shock front obtained in our SPH calculations as function of the entropy ratio(dots). Each point corresponds to the different initial condition. They are compared with the Hugoniot-Rankine relation Eq.(109) (curve). The accordance shows that our SPH calculation reproduces faithfully the conservation of kinetic energy and momentum of the flow through the shock front.

In the usual hydrodynamic computations using space grids, the symmetry of the problem is often a crucial factor to perform a calculation of reasonable size. The SPH method cures this aspect and furnishes a robust algorithm particularly appropriate to the description of processes where rapid expansions of the fluid should be treated. The equations of motion are derived by a variational procedure from the SPH model action with respect to the Lagrangian comoving coordinates. This guarantees that the method furnishes the maximal efficiency for a given number of degrees of freedom, keeping strictly the energy and momentum conservation. For this reason, solutions can be obtained with a very reasonable precision, with a relatively small number of SPH particles. This is the basic advantage of the present method, when we analyze the event-by-event dynamics of the relativistic heavy-ion collisions.

On the other hand, the precision of this method increases rather slowly with the number of SPH particles. Therefore, a relatively large number of particles is required if one wants a very precise numerical solution. However, for the application to the RHIC physics, we may need rather crude precision especially if we consider the dubious validity of the rigorous hydrodynamics. For a calculation with typically 10% errors, the SPH algorithm presented here furnishes a very efficient tool to study the flow phenomena in the RHIC physics.

A fundamental difficulty of the relativistic hydrodynamics for viscous fluid [45,46] is that the dissipation term causes an intrinsic instability to the system. This instability basically comes from the fact that the dissipation term contains q = ¶^mu_m (see Eqs.(106,108)), so that it necessarily introduces the third time-derivative into the equation. This means that we have to specify, at least, a part of the acceleration as the initial condition. Even we specify the initial acceleration, the requirement of the internal self-consistency among the equations above leads to intrinsically unstable solutions. Israel proposed [45,46] to cure these difficulties by introducing higher-order thermodynamics with respect to deviations from the equilibrium. Recently this ''second order thermodynamics'' formalism was discussed in the context of Bjorken type solution [47]. In the examples presented in the present paper, we did not address this question and simply estimated the quantity q from the quantities one time step before. In practice, this will cause no numerical instability and the behavior of the solution is quite satisfactory.

In spite of the above conceptual difficulties when nonadiabatic process is involved, the SPH approach has a nice feature as its flexibility, allowing the treatment of problems with initial conditions without any symmetry, as happens in small systems as those resulting in relativistic nuclear collisions. It should be stressed that, due to the use of Lagrangian coordinates, the method is most suitable for explosive processes like the relativistic heavy ion collisions. (Lagrangian coordinates have been used for treatment of relativistic nuclear collisions also by Nonaka, Honda and Muroya [14]). Furthermore, the variatioal approach guarantees that the SPH equations (85) give the optimal description of motions for a given total number of ''particles'' {r_i(t)}, which are our parameters. In this approach, no numerical instabilities will occur, since the whole system is a Lagrangian system. A numerical code, called SPheRIO has been developped by us on the basis of this algorithm. We shall discuss, in Sec., some results obtained using this code.

5 Decoupling criteria

5.1 Cooper-Frye prescription

As mentioned in the Introduction, the decoupling process is customarily described using the Cooper-Frye prescription [7], which gives the invariant momentum distribution as

This description of decoupling introduces a sharp freezeout hypersurface s, usually characterized by a constant temperature T_f.o.. Before crossing it, particles have a hydrodynamical behavior and, when they cross it suddenly decouple, free-streaming toward the detectors, keeping memory of the conditions (flow, temperature) of where and when they crossed the three dimensional surface.

In SPH representation, we write

where the summation is over all the SPH particles, which should be taken where they cross the hyper-surface T = T_f.o. and n_jm is the normal to this hyper-surface.

Another often used procedure is to take such a freezeout temperature not only constant for a given energy but also energy-independent.

Though operationally simple, and actually useful for obtaining a nice comprehension of several aspects of the phenomena, such a concept of sharp freezeout hypersurface and also of a constant freezeout temperature are clearly highly idealized when applied to finite-volume and finite-lifetime systems as those formed in high-energy heavy-ion collisions.

5.2 Finite-size effect

Before going further, let us for a moment assume that such a freezeout temperature is meaningful. At least, as an average temperature, it should exist. Then, how can we estimate it based on the properties of the system? A simple and natural criterion has already been given by Landau [1], by which a particle decouples when its mean free-path l in the medium becomes larger than the system size L,

This means that T_f.o. is not an intrinsic thermodynamic property of the fluid, but it depends also on the size of the system. In [48], we applied this idea to estimate T_f.o. as function of the incident energy both for pp (p) and nucleus-nucleus collisions, obtaining approximately

Here, the energy dependence appears as a consequence of the increase in the initial energy density, which implies longer expansion time, so larger size L of the system (both longitudinally and tranversally) at the moment of decoupling, requiring lower density, so lower decoupling temperature, too. Fig.14 shows the comparison made in [48] of an estimate of the incident-energy dependence of T_f.o., using Landau's criterion mentioned above, with T_f.o. obtained in a data analysis of p and K transverse-momentum spectra in pp (and p) collisions, in terms of a hydrodynamic parametrization of transverse velocity distribution and temperature. The data were taken from [49,50,51].

An independent estimate of T_f.o.(s) for pp was made by Navarra et al. [52], based on somewhat different but related freeze-out criterion

t_hydro ~ t_col,

where t_hydro and t_col are, respectively, characteristic time for hydrodynamic expansion and particle collision, obtaining a similar result.

When Au + Au collisions data began to appear, Xu and Kaneta [53] reported that T_f.o. seems to decrease with also in nuclear collisions at high energies, and we could verify that the reported results are consistent with T_f.o. ~ ()^-1/12. It is interesting that, more recently [54], in trying to understand the experimentally observed anomalous behavior of the inverse slope parameter T^* of kaon transverse-momentum spectra in central Pb + Pb (Au + Au) collisions [55], we could succeed to obtain the increase in T^* when going from SPS to RHIC domain only with decreasing freezeout temperature T_f.o.(s) with increasing as shown in Fig.15 and Table1. Here, SPheRIO code has been used with averaged IC. This is because, as will be shown in Sec., Figs.19, ²⁰ , ²¹ , ²² , the slope parameter is not sensitive to IC fluctuations.

5.3 Continuous emission

The decoupling temperature discussed in the preceding paragraphs should actually be interpreted as an average value of such a temperature. When applied to finite-volume and finite-lifetime systems as those formed in high-energy heavy-ion collisions, a sharply defined freeze-out hypersurface is clearly too idealized and for a more precise analysis of data, one needs a more elaborate description of the decoupling process. In [57], we introduced a method that we call Continuous Emission (CE) and, as compared to the usual Cooper-Frye one, we believe closer to what happens in the actual collisions. The essential point is the introduction of momentum-dependent escape probability

of a particle from a space-time point x without collision in the medium, so that the emission may occur from any point of the fluid and at any time, according to this probability. This means that we are interpreting probabilistically the Landau condition, (112), as should be and also giving the system size L a more precise meaning, namely, the quantity of matter the escaping particle encounters in his trajectory. The integral above is evaluated in the proper frame of the particle. Then, the distribution function f(x,k) of the expanding system has two components, one representing the portion of the fluid already free and another corresponding to the part still interacting, i.e.,

We may write the free portion as

The inclusive one-particle distribution, for instance, is then written as

where the surface term corresponds to particles already free at the initial time.

As particles can be emitted in different stages of the fluid expansion, it is natural that the appearance of the observables in CE becomes different from that of the usual Cooper-Frye prescription at a constant temperature [7]. In general, in CE, the large-p_T particles are mainly emitted at early times when the fluid is hot and mostly from its surface, whereas the small-p_T components are emitted later when the fluid is cooler and from a larger spatial domain. We shall discuss, in Sec.6, the manifestation of this effect in two-pion interferometry.

More detailed description of this method and several of its predictions are accounted for by Grassi in a separate paper of this issue [58].

6 Applications

In the following, we shall present some applications of NeXuS+SPheRIO code [59,60,61], which has been constructed by coupling the NeXuS event generator to SPheRIO code, described in Subsection4.3.

As mentioned in the Introduction, the main ingredients of hydrodynamic calculations are IC, EoS and the decoupling procedure. None of these are well known at the present moment.

Since we are mainly concerned with the effects of IC fluctuations and consequences of different decoupling descriptions, here we just take the IC created by NeXuS event generator and use the equations of state described in Section 3.5. The strangeness conservation has not been taken into account. As for the decoupling procedure, we adopted both the conventional sharp freeze-out prescription and the continuous emission description, because one of our purpose is to see the differences resulting from these two descriptions.

We do not expect that these options will reproduce all the experimental data. In fact, we found that the use of the same version of NeXuS code for high- and low-energy nuclear collisions caused some discrepancies in reproducing the rapidity distribution of charged particle. Therefore, we have introduced additional parameters to ajust at least the overall rapidity distributions.

6.1 Effects of fluctuating initial conditions

In Sec.2, we stressed that, due to the finite size of the systems, large fluctuations are expected in the initial stage of actual nuclear collisions, even for a fixed impact parameter, and that IC generated by realistic event generators do show such effects [12,15]. Let us see in this Section what are the effects of fluctuating IC on some of the observables.

In [12], by solving the hydrodynamic equations with longitudinal boost-invariance, it has been shown that the bumpy IC i) develop azimuthally asymmetric flows, even for central collisions, because there is no symmetry in each fluctuating event; and also ii) enhance high-p_T direct photon yields due to the high temperature in the blobs.

Since our interest here is just to show the effects of fluctuating IC, we chose the simplest decoupling criterion, namely the usual Cooper-Frye sudden decoupling with freezeout temperature T_fo in this Subsection.

6.1.1 Rapidity and m_T distributions

Let us first consider Pb+Pb collisions at SPS. In Fig.16, we show the rapidity distributions for negative particles and p-, respectively, for the most central Pb+Pb collisions at = 17.3 A GeV. Each event, computed from randomly generated IC like the one shown in Figs.1 and ² (left), is represented by a thin curve (for each type of particle, either negatives or p-). The thick solid lines represent the averages over 50 events, with the corresponding dispersions. We compare the average distributions with the distributions computed starting from the averaged IC like the one shown in Figs.1 and ² (right), represented by dashed lines here. Here, we took T_fo = 140MeV. The data points are shown for comparison [62].

As is seen in this Figure, the rapidity distributions show large fluctuations from event to event and, although similar, there is a non-negligible difference between the average distributions and the ones obtained from the average IC, especially in the case of negative particles, constituted mostly of pions. For the same average initial energy, the multiplicity decreases about 6 ~ 7% if fluctuations exist, confirming what we showed in Sec. 2. Finally, the data points in Fig.16 are closer to the averaged distribution over fluctuating IC.

What happens with other centralities is similar as seen in Figs.17 and 18. Some differences are: i) as the collisions become more peripheral, naturally the dispersions increase; ii) the differences between the average distributions and the ones produced by average IC seems to decrease.

In Fig.19, we show the corresponding m_T distributions, for the most central Pb+Pb collisions at = 17.3 A GeV. One can see that apparently the fluctuation effects on the transverse-momentum spectra are very small and the averaged spectra are in good agreement with those computed with averaged IC, and also with data. This smallness of the fluctuation effects on the transverse-momentum spectra is notin contradiction with the large fluctuations we saw above of rapidity distributions. In part this is due to the logarithmic scale used here. One can conclude, however, that the IC fluctuations affect very little the slope of the m_T spectra. This conclusion is valid for other centralities, as seen in ^Figs.20 and ²¹ , except for the most peripheral case, where the dispersions are very large. If we look more carefully, we can perceive some difference between the averaged spectra and those computed with the averaged IC, which appears in all the centralities, that is the tail of the averaged spectra is more concave and the distributions become higher, probably because the expansion is more violent in this case, due to the high density spots in the IC.

At RHIC energies the results are similar. In ^Fig.22 , we show the p_T distributions of charged particles in the most central Au+Au collisions at 130A GeV, calculated both with NeXuS fluctuating IC and with averaged IC. Like in the case of Pb+Pb collisions at SPS, we used sudden freezeout, but with lower freezeout temperature, T_fo = 128MeV, the same value used in [54] to fit the inverse slope parameter of kaons, in the same collisions, as shown in Fig.15 and Table1 of Section 5. As is seen, the fluctuation effects are small in p_T distribution, both curves agree each other and with data, and the average over fluctuating events gives a slightly more concave shape, just like in Pb+Pb collisions at SPS.

We show in Fig.23, the pseudorapidity distribution of charged particles in the most central Au+Au collisions at 130A GeV, calculated with NeXuS fluctuating IC and the corresponding averaged IC, and compared with data [63]. Qualitatively, the behavior is similar to the case of Pb+Pb collisions at SPS, namely the average h distribution is close to the distribution with the averaged IC, being the latter slightly higher than the former. The difference here is smaller than at the lower energy case. Probably this is due to the decrease of DE_i/ < E > in Eq.(9) as the incident energy increases.

6.1.2 Elliptic-flow parameter v₂

Let us turn to the elliptic-flow parameter v₂, defined as the second Fourier coefficient of the azimuthal distribution [65]

Thus,

where the bracket denotes the average value and y gives the event-dependent collision plane.

Usually, v₂ parameter is interpreted as indicating a flow asymmetry caused by the initial-condition asymmetry associated with the non-zero impact parameter. In this case, as the produced matter in the collision is likely to be flattened in the impact-parameter direction, the pressure gradient would be larger in this direction and so would be the flow. However, as mentioned at the top of this section, it was shown in Ref. [12] that, even for central collisions, fluctuating IC develop azimuthally asymmetric flows, because in this case there is no symmetry in each event and, experimentally, the impact parameter cannot be determined unambiguously. So, one expects larger v₂ for fluctuating IC as compared to averaged IC case. Remark that the fluctuation we are talking about is not the often discussed [65,66] finite-multiplicity fluctuation at the end of the process. In Fig.24, results of NeXuS+SPheRIO for v₂ of pions produced in Pb+Pb collisions at 17.3A GeV are shown, both with fluctuating IC and averaged IC and compared with data [62]. Large discrepancy is seen between the two ways of calculating this parameter, and the data are closer to < v₂ > , calculated with fluctuating IC.

6.1.3 Two-pion interferometry

As is well known, the identical-particle correlation, also known as Hanbury-Brown-Twiss effect (HBT effect) [67] is a powerful tool for probing geometrical sizes of the space-time region from which they were emitted. If the source is static like a star, it is directly related to the spatial dimensions of the particle emission source. When applied to a dynamical source, however, several non-trivial effects appear [68,69], reflecting its time evolution as happens in high-energy heavy-ion collisions. Being so, the inclusion of IC fluctuations may affect considerably the so-called HBT radii, because, as discussed in Sec. 2, the IC in the event-by-event base often show small high-density spots in the energy distribution, and our expectation is that such spots manifest themselves at the end when particles are emitted, giving smaller HBT radii.

We shall discuss here only the recent application of NeXuS+SPheRIO code for HBT effect with IC fluctuations, for Au+Au collisions at 130AGeV [61]. A more detailed account of two-particle correlations is given by Padula [70] in a separate paper in this issue.

For studying the fluctuation effects of IC on HBT correlation, first we assume Cooper-Frye sudden freezeout (FO) at T_f.o. = 128MeV. This freezeout temperature is the same one previously found by studying the energy dependence of kaon slope parameter T* [54], discussed in Subsection 5.2 and appears in Table1. We showed in Subsection 6.1.1 that indeed this choice of T_f.o. reproduces the charged p_T-distribution data (including pions) at 130AGeV, both with averaged and fluctuating IC. We also neglect the resonance decays. It is argued [71] that, since resonance decays contribute to the correlations with very small q values (q q_min, where q_min is the minimum measureable q), the experimentally determined HBT radii are essentially due to the direct pions. Then the two-particle correlation function is expressed in terms of the distribution function f(x,p) as

where P = (p₁+p₂)/2 and q = (p₁-p₂) and p_i is the momentum of the ith pion. Usually

In SPH representation, we write I(q,P) as

where the sumation is over all the SPH particles. In the Cooper-Frye freezeout, these particles should be taken where they cross the hyper-surface T = T_f.o. and n_jm is the normal to this hyper-surface. Notice that, if we put p₁ = p₂ = P, so that q = 0, Eqs. 117 and 118 are reduced, respectively, to Eqs. 106 and 107, that is, to the inclusive one-particle distribution.

In Fig.25, we compare the correlation function C₂ averaged over 15 fluctuating events with those computed starting from the averaged IC (so, without fluctuations). One can see that the IC fluctuations are reflected in large fluctuations also in the HBT correlations. When averaged, the resulting correlation < C₂ > are broader than those computed with averaged IC, so giving smaller radii as expected. Also the shape of the correlation functions changes.

We plot the m_T dependence of HBT radii, with Gaussian fit of C₂, in Fig.26, together with RHIC data [72,73] and results with CE, which will be discussed in Sec. 6.2.2. It is seen that the smooth IC with sudden FO makes the m_T dependence of R_o flat or even increasing, which is in agreement with other hydro calculations [8] but in conflict with the data. The fluctuating IC make the radii smaller, especially in the case of R_o, however without changing the m_T-dependence.

6.2 Effects of continuous emission

As discussed in Sec. 5, it is more likely that the decoupling occurs not suddenly, but continuously from every point of the fluid and at every instant of time, according to some escaping probability (x,p) given by Eq.(114). There are several nice predictions of the Continuous Emission Model (CEM) as discussed in [58]. However, although more realistic, this description is not handy because, depends on the momentum of the escaping particle and, moreover, on the future of the fluid as seen in Eq.(114). In order to make the computation practicable, in [61] we first took on the average, i.e.,

(x,p)Þ < (x,p) > º (x).

Then, approximated linearly the density

r(x') = as(x')

(where, for example in the ideal massless pion gas case, a = (45z(3))/2p⁴ = .278) in the integral of Eq.(114). Thus,

where < sv > has been included in k = 0.5a < sv > . Although approximately, now we can compute in each space-time point.

In Fig.27, we show the time evolution of the escaping probability (x) given by Eq.(123) in the mid-rapidity plane for the most central Au+Au collisions at 130A GeV. The IC have been computed at t = 1fm and averaged over 30 NeXuS events. Here and in the computations of observables below, the parameter k has been estimated to be .3, corresponding to < sv > = 2fm² in the zero temperature limit and some 20% larger at T = m_p. As seen, the probability remains 0.1 < (x) < 0.8 in a quite large domain, especially for t > 10fm, indicating that both the emission zone and duration are expected to be large, in opposition to the standard sudden freezeout case. For comparison, we show also the temperature distribution in Fig.27, with some isotherms.

To calculating the spectra, now Eq.(117) is translated into SPH language, which is given by Eq.(111), but where the sum is computed in the present case not over T = T_f.o. hypersurface but picking out SPH particles according to the probability , given by Eq.(123) with the normal n_jm pointing to the 4-gradient of . Since our procedure favors emission from fast outgoing SPH ''particles'', because r decreases faster there and so does s in this case making larger, we believe that the main feature of CEM is preserved in our approximation.

6.2.1 m_T distributions

In Fig.28, we show the charged m_T distribution at mid-rapidity in the most central Au+Au collisions at 130A GeV, computed by using CEM, with and without fluctuations. Compare this figure with ^Fig.22 , where the same data are compared with calculations using sudden freezeout. One sees that, although the emission mechanisms are quite different, the results are similar, for the choice of the parameters, T_fo in one case and k in the other. However, the origins of the spectrum shape are different. Whereas, in the freezeout case, all the particles are emitted at the same temperature T_fo and the concave shape of the spectra is due to different transverse velocities of the fluid at different instants of time, in the continuous emission case, the large-p_T particles are in general emitted earlier and at higher temperature, the small-p_T particles are emitted later, when the fluid is cooler.

6.2.2 Two-pion interferometry

Let us now consider the effect of continuous emission on HBT effect. Since HBT is sensitive to the space-time geometry of the fluid, and CEM produces important modification of the emission zone, we expect considerable changes in the so-called HBT radii. As mentioned above, according to this picture, the large-p_T particles are mainly emitted at early times when the fluid is hot and mostly from its surface, whereas the small-p_T components are emitted later when the fluid is cooler and from larger spatial domain.

In [74], we considered this effect, using the boost-invariant solution [30], and showed that whereas the so-called side radius is independent of the average p_T, the out radius decreases with < p_T > , because of the reason mentioned above. This behavior is expected to essentially remain in the general 3-dimensional expansion, described by SPheRIO code.

To compute the correlation function C₂(q,P) in CEM, we first rewrite the integral(121) as

which is similar to Eq.(117). Then, translate it into SPH language by using Eq.(122), but with the sum evaluated by picking out SPH particles according to the probability , given by Eq.(119) with the normal n_jm pointing to the 4-gradient of , exactly in the same way as done in calculating the inclusive spectrum.

We plot in Fig.26, some results for the m_T dependence of HBT radii computed in this way, both with fluctuating IC and averaged IC, for Au+Au collisions at 130AGeV. For comparison, we show also the results with sudden freezeout and data [72,73]. Comparing the averaged IC case (with CEM (CE1)), with the corresponding freezeout (FO1), one sees that, while R_L remains essentially the same, R_s decreases faster and as for R_o, it decreases now inverting its m_T behavior.

In Fig.26, one can also see the combined effect of fluctuating initial conditions together with continuous emission (curve CE2). All the radii are smaller than CE1, with averaged IC (but using CEM), as happend in the sudden freezeout case. The agreement with data is excelent both for R_L and R_s, and improved considerably the results for R_o with respect to the usual hydrodynamic description.

7 Conclusions and Outlook

In this survey, we discussed on several aspects of applications of hydrodynamic model to nucleus-nucleus collisons, giving especial emphases on i) method of solving the hydrodynamic equations for arbitrary configurations; ii) accounting for the probable event-by-event fluctuations of the initial conditions; and iii) the decoupling criteria for obtaining the observables.

The Smoothed-Particle Hydrodynamic approach, using a special hyperbolic coordinate system, is particularly interesting tool for describing rapidly expandig matter as that in high-energy nucleus-nucleus collisons, because of its efficiency and flexibility, as demonstrated through examples considered in Sec. 4.3.5, 4.3.6 and through applications, shown in Sec. 6, where the systems, in general, do not present any symmetry.

The initial conditions for hydrodynamic expansion of matter formed in relativistic heavy-ion collisions are likely neither smooth nor symmetrical, because of small size of these systems. This property has also been shown by some event generators, which take the microscopic dynamics into account. The fluctuations in the initial conditions may produce large discrepancies in the computed results for some observables, in comparison with those obtained with the usual averaged, smooth and symmetrical initial conditions. So, it is our opinion that, to extracting correct information from certain experimental data, the inclusion of the fluctuations mentioned above is mandatory.

We find that, although Cooper-Frye formula can give good results in many cases, a more realistic decoupling description is needed for obtaining many other observables. A description we proposed is the Continuous Emission Model. One typical example we showed here, in which this model improve considerably the description of data is the m_T dependence of HBT radii, observed at RHIC.

In addition to the points summarized above, there are several basic questions to be addressed in the hydrodynamical scenary. First of all, although for the bulk properties of observables so far studied the ideal hydrodynamical description works fairly well, we still do not know how the non-equilibrium processes affect the results of these studies. This question refers both to the initial conditions and to the final particle decoupling process. For example, the present study of continuous emission showed that the real particle emission mechanism seems to be very far from the conventional Cooper-Frye scenary of sharp freeze-out surface. In Fig.27, we showed that the escape probability varies rather slowly both in space and time, and there is no indication of a sharp freezeout hypersurface. Another way of showing the same is given in Fig.29, where we plotted the entropy vs. temperature of the SPH particles at the emission point. The temperature values corresponding to the particle emission spread widely. This shows that, even for a large system (Au+Au) there exists a substantial part of the system which remains out of equilibrium for a long time. Up to now, no dynamical reaction of these non-equilibrium components on the hydrodynamical evolution is studied. In this sense, it is very important to develop transport theoretical investigations [75,76] related to hydro description.

Another interesting type of phenomena related to the hydrodynamical model is the ocurrence of instabilities associated to the phase transition and surface tension (finite-size effects), such as spinodal instabilities [77]. Some studies in this direction is in progress.

Acknowledgments

The authors are grateful to all the SPheRIO team, namely, F. Grassi. C.E. Aguiar, T. Osada, B.M. Tavares, L.L.S. Portugal, G. Grunfeld, R. Andrade, for their collaboration and stimulating discussions. Especially, the most of the hard task in the development of the code is due to T. Osada. The main part of the consideration on hadronic resonance is due to C.E. Aguiar. B.M. Tavares, L.L.S. Portugal, A.G. Grunfeld helped the construction of equations of state. We are indebted to K. Werner for valuable discussions on NeXuS event generator and for providing us with the correspondent code. Analytical solutions have been obtained in a fruitful collaboration with T. Csörgo. We are thankful to Y.M. Sinyukov and L.P. Csernai for profitable discussions on transport theoretical description of continuous emission. We also appreciated the discussions and the encouragements received from F. Pottag, S.S. Padula, F.S. Navarra, S.B. Paiva, M. Makler, M. Gorenstein, M. Ga\'zdzicki, J. Rafelski, H.T. Elze and E. Fraga.

We acknowledge financial support by FAPESP (1990/4074-5, 1993/2463-2, 1995/4635-0, 1998/2249-4, 2000/04422-7 and 2001/09861-1), CAPES/PROBRAL, CNPq, FAPERJ and PRONEX.

References

[1] L. D. Landau, Izv. Akad. Nauk SSSR 17, 51 (1953); S. Z. Belenkij and L. D. Landau, Usp. Fiz. Nauk 56, 309 (1955); in ''Collected papers of L. D. Landau" (ed. D. Ter-Haar, Pergamon, Oxford, 1965) 569 and 665.

[2] E. Fermi, Prog. Theor. Phys. 5, 570 (1050).

[3] P. Carruthers, Annals N.Y. Acad. Sci. 229, 91 (1974).

[4] T. Hirano, J. Phys. G 30, S845 (2004).

[5] P. Huovinen, nucl-th/0305064.

[6] P. F. Kolb and U. Heinz, nucl-th/0305084.

[7] F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).

[8] K. Morita, S. Muroya, C. Nonaka, and T. Hirano, Phys. Rev. C66, 054904 (2002).

[9] S. Paiva, Y. Hama, and T. Kodama, Phys. Rev. C 55, 1455 (1997).

[10] G.N. Fowler, F. Navarra, M. Plümer, A. Vourdas, R.M. Weiner, and G. Wilk Phys. Rev. C 40, 1219 (1989); See also F.S. Navarra, Braz. J. Phys. 35, 3 (2005).

[11] I.M. Khalatnikov, Zh. Eksp. Teor. Fiz. 27, 529 (1954); see also [1].

[12] M. Gyulassy, D.H. Rischke, and B. Zhang, Nucl. Phys. A613, 397 (1997).

[13] B.R. Schlei and D. Strotman, Phys. Rev. C59, 9 (1999).

[14] C. Nonaka, E. Honda, and S. Muroya, Eur. J. Phys. C17, 663 (2000).

[15] H.J. Drescher, F.M. Liu, S. Ostapchenko, T. Pierog, and K. Werner, Phys. Rev. C65, 054902 (2002).

[16] V.K. Magas, L.P. Csernai, and D. Strottman, Phys. Rev. C64, 014901 (2001).

[17] K.J. Eskola, P.V. Ruuskanen, S.S. Räsänen, and K. Tuominen, Nucl. Phys. A 696, 715 (2001).

[18] T. Hirano and Y. Nara, nucl-th/0403029.

[19] L.D. Landau and E.M. Lifshitz, Fluid Mechanics, second ed. Butterworth-Heinemann Ltd. Oxford (U.K.).

[20] J. Berges, S. Borsanyi and C. Wetterich, hep-ph/0403234.

[21] T. Kodama, H.T. Elze, C.E. Aguiar, and T. Koide, cond-mat/0406732.

[22] T. Kodama, Introduction to Relativistic Gases, in AIP Conf. Proc. 631, 3 (2003).

[23] See for example, P. Braun-Munzinger, K. Redlich and J. Stachel, nucl-th/0304013.

[24] E. Beth and G.E. Uhlenbeck, Physica 3, 729 (1936).

[25] R. Dashen, S. Ma, and H. Bernstein, Phys. Rev. 187, 345 (1969).

[26] D.H. Rischke, M.I. Gorenstein, H. Stöcker, and W. Greiner, Zeit. Phys. 51, 485 (1991); G.D. Yen, M.I. Gorenstein, W. Greiner, and S.N. Yang, Phys. Rev. C56, 2210 (1997).

[27] J. Sollfrank and U. Heinz, arXiv:nucl-th/9505004.

[28] Z. Fodor and S.D. Katz, arXiv:hep-lat/0106002.

[29] K. Rajagopal, arXiv:hep-ph/0009058.

[30] R.C. Hwa, Phys. Rev. D10, 2260 (1974); C.B. Chiu and K.-H. Wang, Phys. Rev. D12, 272 (1975); C.B. Chiu, E.C.G. Sudarshan, and K.-H. Wang, Phys. Rev. D12, 902 (1975).

[31] J.D. Bjorken, Phys. Rev. D27, 140 (1983).

[32] T. Csörgo, F. Grassi, Y. Hama, and T. Kodama, Phys. Lett. B 565, 107 (2003).

[33] T. Csörgo, L.P. Csernai, Y. Hama, and T. Kodama, Heavy Ion Phys. 21, 73 (2004); nucl-th/0306004.

[34] H.-T. Elze, Y. Hama, T. Kodama, M. Makler, and J. Rafelski, J. Phys. G25, 1935 (1999); Heavy Ion Phys. 10, 275 (1999); Proceedings of the IVth. Workshop on Quantum Chromodynamics, Paris 1998, eds. H.M. Fried and B. Müller (World Scientific, Singapura, 1999) 227, hep-ph/9809570.

[35] L.B. Lucy Astrophys. J. 82, 1013 (1977).

[36] R.A. Gingold and J.J. Monaghan, Mon. Not. R. Astro. Soc. 181, 375 (1977).

[37] C.E. Aguiar, T. Kodama, T. Osada, and Y. Hama, J. Phys. G27, 75 (2001); T. Kodama, C.E. Aguiar, T. Osada and Y. Hama, J. Phys. G27, 557 (2001).

[38] A.H. Taub, Phys. Rev. 94, 1468 (1954); J.W. Herivel, Proc. Camb. Phil. Soc. 51, 344 (1955); B.F. Schutz Jr., Phys. Rev. D2, 2762 (1970); J.R. Ray, J. Math. Phys. 13, 1451 (1972); G.Yu. Bogoslovskii, Sov. Phys. Dokl. 31, 879 (1986); T. Chiueh, Phys. Rev. E49, 1269 (1994).

[39] L. Mittag, M.J. Sephen and W. Yourgrau, in Variational Principles in Dynamics and Quantum Theory, ed. W.W. Yourgrau and S. Mandelstam (Dover, New York, 1968).

[40] V.M. Strutinski, Nucl. Phys. A95, 420 (1967).

[41] Y. Hama and F. W. Pottag, Resolution of Hydrodynamical Equations for Transverse Expansions - II, University of São Paulo print IFUSP/P-481, Aug./84.

[42] V. Schneider, U. Katscher, D. H. Rischke, B. Waldhauser, J. A. Maruhn, and C. D. Munz, J. Comput. Phys. 105, 92 (1993).

[43] S. Siegler and H. Riffert, astro-ph/9904070.

[44] J. Von Neumann and R. D. Richtmyer, J. Appl. Phys. 21, 232 (1950).

[45] W. Israel, Ann. Phys. 100, 310 (1976); W. Israel and J.M. Stewart, Ann. Phys. 118, 341 (1979).

[46] W.A. Hiscock and L. Lindbolm, Phys. Rev. D 31, 725 (1985); Phys. Rev. D35, 3723 (1987).

[47] A. Muronga, Phys. Rev. Lett. 88, 062302 (2002); A. Muronga, Phys. Rev. C69, 034903 (2004).

[48] Y. Hama and F.S. Navarra, Z. Phys. C53, 501 (1992).

[49] British-Scandinavian Collab., B. Alper et al., Nucl. Phys. B100, 237 (1975); F.W. Büsser et al., Phys. Lett. B46, 471 (1973) and Nucl. Phys. B 106, 1 (1976); K. Eggert et al., Nucl. Phys. B98, 49 (1975); British-Scandinavian-MIT Collab., K. Guettler et al., Phys. Lett. B64, 111 (1976); Nucl. Phys. B116, 77 (1976).

[50] UA1 Collab., C. Albajar et al., Nucl. Phys. B 335, 261 (1990); UA2 Collab., M. Banner et al. Phys. Lett. B 115, 59 (1982); Phys. Lett. B122 322 (1983); Z. Phys. C27, 329 (1985); UA5 Collab., G.J. Alner et al., Nucl. Phys. B258, 505 (1985); UA5 Collab., R.E. Ansorge et al., Phys. Lett. B199, 311 (1987).

[51] T. Alexopoulos et al., Phys. Rev. Lett. 64, 991 (1990).

[52] F.S. Navarra, M.C. Nemes, U. Ornik, and S.M. Paiva, Phys. Rev. C45, R2552 (1992).

[53] Nu Xu and M. Kaneta, Nucl. Phys. A698, 306 (2002).

[54] M. Ga\'zdzicki, M.I. Gorenstein, F. Grassi, Y. Hama, T. Kodama, and O. Socolowski Jr., Braz. J. Phys. 34, 322 (2004); Y. Hama, F. Grassi, O. Socolowski Jr., T. Kodama, M. Ga\'zdzicki, and M.I. Gorenstein, Acta Phys. Polon. B35, 179 (2004).

[55] M.I. Gorenstein, M. Ga\'zdzicki, and K.A. Bugaev, Phys. Lett. B567, 175 (2003).

[56] E866 Collab., L. Ahle et al., Phys. Lett. B490, 53 (2000); NA49 Collab., S.V. Afanasiev et al., Phys. Rev. C66, 054902 (2002); BRAHMS Collab., D. Ouerdane et al., Nucl. Phys. A715, 478 (2003); STAR Collab., C. Adler et al., nucl-ex/0206008.

[57] F. Grassi, Y. Hama, and T. Kodama, Phys. Lett. B355, 9 (1996); F. Grassi, Y. Hama, and T. Kodama, Z. Phys. C 73, 153 (1996).

[58] F. Grassi, Braz. J. Phys. 35, 52 (2005).

[59] T. Osada, C.E. Aguiar, Y. Hama, and T. Kodama, nucl-th/0102011; C.E. Aguiar, Y. Hama, T. Kodama, and T. Osada, Nucl. Phys. A698, 639 (2002).

[60] C.E. Aguiar, Y. Hama, T. Kodama, and T. Osada, in AIP Conf. Proc. 631, 686 (2003); ''Event-by-event fluctuation of the initial conditions in hydrodynamic model'', Y. Hama, F. Grassi, O. Socolowski Jr., C.E. Aguiar, T. Kodama, L.L.S. Portugal, B.M. Tavares, and T. Osada, in Proc. of 32 ISMD, eds. A. Sissakian, G. Kozlov and E. Kolganova, (World Scientific - Singapore, 2003), 65-68.

[61] O. Socolowski Jr., F. Grassi, Y. Hama, and T. Kodama, Phys. Rev. Lett. 93, 182301 (2004).

[62] H. Appelshäuser et al., NA49 Collab., Phys. Rev. Lett. 82, 2471 (1999); A.M. Poskanzer, nucl-ex/0110013.

[63] PHOBOS Collaboration, B.B. Back et al., Phys. Rev. Lett. 87, 102303 (2001).

[64] STAR Collab. C. Adler et al., Phys. Rev. Lett. 87, 112303 (2001).

[65] S. Voloshin and Y. Zhang, Z. Phys. C70, 665 (1996).

[66] J.-Y. Ollitrault, nucl-ex/9711003.

[67] R. Hanbury-Brown and R.Q. Twiss, Phil. Mag. Ser. 7, Vol. 45, 663 (1954).

[68] Y. Hama and Sandra S. Padula, Phys. Rev. D 37, 3237 (1988); Bose-Einstein Correlation in Landau's Model, Y. Hama and Sandra S. Padula, in Hadronic Matter in Collision, eds. P. Carruthers and D. Strottman (World Scientific Pub. Co., Singapore, 1986) 63.

[69] S. Pratt, Phys. Rev. D 33, 1314 (1986); K. Kolehmainen and M. Gyulassy, Phys. Lett. B 180, 203 (1986); A.N. Makhlin and Yu.M. Sinyukov, Z. Phys. C 39, 69 (1988).

[70] Sandra S. Padula, Braz. J. Phys. 35, 70 (2005).

[71] S. Nickerson, T Csörgo, and D. Kiang, Phys. Rev. C57, 3251 (1998).

[72] STAR Collab., C. Adler et al., Phys. Rev. Lett. 87, 082301 (2001).

[73] PHENIX Collab., K. Adcox et al., Phys. Rev. Lett. 88, 192302 (2002).

[74] F. Grassi, Y. Hama, S.S. Padula, and O. Socolowski Jr., Phys. Rev. C62, 044904 (2000).

[75] Yu. M. Sinyukov, S. V. Akkelin, and Y. Hama, Phys. Rev. Lett. 89, 052301 (2002).

[76] C. Anderlik et al. Phys. Rev. C59, 3309 (1999); V. K. Magas et al., Acta Phys. Hung. New Ser. Heavy Ion Phys. 9, 193 (1999); V. K. Magas et al., Phys. Lett. B459, 33 (1999).

[77] C.E. Aguiar, E.S. Fraga, and T. Kodama, Braz. J. Phys. 34, 307 (2004).

Received on 15 November, 2004

Publication Dates

History