A Fully Integrated Transport Approach to Heavy Ion Reactions with an Intermediate Hydrodynamic Stage
Abstract
We present a coupled Boltzmann and hydrodynamics approach to relativistic heavy ion reactions. This hybrid approach is based on the Ultrarelativistic Quantum Molecular Dynamics (UrQMD) transport approach with an intermediate hydrodynamical evolution for the hot and dense stage of the collision. Eventbyevent fluctuations are directly taken into account via the nonequilibrium initial conditions generated by the initial collisions and string fragmentations in the microscopic UrQMD model. After a (3+1)dimensional ideal hydrodynamic evolution, the hydrodynamical fields are mapped to hadrons via the CooperFrye equation and the subsequent hadronic cascade calculation within UrQMD proceeds to incorporate the important final state effects for a realistic freezeout. This implementation allows to compare pure microscopic transport calculations with hydrodynamic calculations using exactly the same initial conditions and freezeout procedure. The effects of the change in the underlying dynamics  ideal fluid dynamics vs. nonequilibrium transport theory  will be explored. The freezeout and initial state parameter dependences are investigated for different observables. Furthermore, the time evolution of the baryon density and particle yields are discussed. We find that the final pion and proton multiplicities are lower in the hybrid model calculation due to the isentropic hydrodynamic expansion while the yields for strange particles are enhanced due to the local equilibrium in the hydrodynamic evolution. The results of the different calculations for the mean transverse mass excitation function, rapidity and transverse mass spectra for different particle species at three different beam energies are discussed in the context of the available data.
pacs:
25.75.q, 24.10.Lx, 24.10.Nz, 25.75.AgI Introduction
One of the main motivations to study high energy heavy ion collisions is the creation of a new deconfined phase of strongly interacting matter, the so called QuarkGluon Plasma (QGP) (1); (2). At the Relativistic Heavy Ion Collider (RHIC) many experimental observations like, e.g., jet quenching and high elliptic flow hint to the fact that a strongly coupled QGP (sQGP) might have been created (3); (4); (5); (6). At CERNSPS energies evidence for the creation of a new state of matter has been published, e.g., the enhanced K/ ratio (’horn’) and the step in the mean transverse mass excitation function for pions, kaons and protons (7). Especially the low energy (high ) program at SPS showed a culmination of exciting results. Therefore this energy regime will be the subject to further detailed studies at the CERNSPS, BNLRHIC, JINRNICA and at the FAIR facility.
Since the direct detection of free quarks and gluons is impossible due to the confining nature of QCD, it is important to model the dynamical evolution of heavy ion reactions to draw conclusions from the final state particle distributions about the interesting early stage of the reaction. One approach which aims at the description of heavy ion reactions consistently from the initial state to the final state is relativistic transport theory (8); (9); (10); (11); (12); (13). This microscopic description has been applied quite successfully to the partonic as well as to the hadronic stage of the collision. Unfortunately, most transport approaches are restricted to scattering processes. Thus, if the particle density increases it becomes questionable if a restriction to twoparticle interaction is still justified. While first attempts to include multiparticle interactions have been proposed (11); (14); (15); (16); (17), this extension of transport theory is still in its infancy. To explain hadronization and the phase transition between the hadronic and the partonic phase on a microscopic level is also one of the main open issues that still has to be resolved. It is therefore difficult to find an appropriate prescription of the phase transition in such a microscopic approach. First, however promising attempts to solve the microscopic hadronization problem can be found in the literature (18); (19); (20); (21); (22); (12).
Hydrodynamics, on the other hand, has been proposed many years ago as a tool for the description of the hot and dense stage of heavy ion reactions where the matter might behave like a locally thermalized ideal fluid (23); (24); (25). In this approach it is possible to model phase transitions explicitly because one of the major inputs to a hydrodynamic calculation is the equation of state (EoS). The hydrodynamic description has gained importance over the last few years because the high elliptic flow values that have been observed at RHIC seem compatible with some ideal hydrodynamic predictions (26); (27); (28). The initial conditions and freezeout prescription are the boundary conditions for a hydrodynamic calculation and therefore a further crucial input. Thus, the hydrodynamic results depend strongly on the initial and final state prescription that is applied in the specific calculation.
To get a more consistent picture of the whole dynamics of heavy ion reactions various so called microscopic plus macroscopic (micro+macro) hybrid approaches have been launched during the last decade. Most noteworthy in this respect are the pioneering studies related to a coupling between UrQMD (Ultrarelativistic Quantum Molecular Dynamics) and hydrodynamics (a detailed systematic investigation of this coupling procedure can be found in the following references (29); (30); (31); (32); (33); (34); (25); (35)).
Other approaches in the same spirit are, e.g., the NEXSpheRIO approach that uses initial conditions calculated in a nonequilibrium model (NEXUS) followed by an ideal hydrodynamic evolution (36); (37); (38); (39); (40); (41); (42) or a hybrid approach by Toneev et al. which uses QGSM initial conditions followed by a threedimensional hydrodynamic evolution (43). In this way eventbyevent fluctuations are taken into account and the calculation mimics more realistically the experimental case. For the freezeout NEXspheRIO employs a continuous emission scenario or a standard CooperFrye calculation. Other groups, e.g., Teaney et al. (44), Hirano et al. (45); (46), Bass/Nonaka (25), are using smooth Glauber or Color Glass Condensate initial conditions followed by a full two or threedimensional hydrodynamic evolution and calculate the freezeout by a subsequent hadronic cascade. The separation of chemical and kinetic freezeout and final state interactions like resonance decays and rescatterings are taken into account. There are two major conclusions from these previous studies: The treatment of the initial state fluctuations and the final decoupling is of major importance for a sound interpretation of the experimental data.
Unfortunately, all presently existing micro+macro approaches rely on a complete separation of the three main ingredients (initial conditions, hydrodynamic evolution, transport calculation). Thus, it is impossible to compare the evolution of the system between hydrodynamics and transport simulation directly and from the same initial conditions. This may provide essential new insights into the role of viscosity and local equilibration. In addition, the usual separation of the program code does not allow for a dynamical coupling between hydrodynamics and transport calculation, which would be desirable to consistently solve the freezeout puzzle (47); (48); (49); (50).
To overcome these restrictions, we go forward and build a transport approach with an embedded threedimensional ideal relativistic one fluid evolution for the hot and dense stage of the reaction. This allows to reduce the parameters for the initial conditions and the freezeout prescription. The aim is to compare calculations with different EoS within the same framework. It will be possible to extract the effect of changes in the EoS, e.g., a phase transition from hadronic matter to the QGP, on observables.
In this paper we describe the specific micro+macro hybrid approach that embeds a hydrodynamic phase in the UrQMD approach. First we explain the initial conditions, then introduce the basics of the hydrodynamic evolution including the hadron gas EoS and the transport calculation and illustrate how the freezeout is treated. In the second part, the sensitivity of the results on the parameters are tested and the time evolution of the baryon density and the particle numbers are compared. Results on particle multiplicities and rapidity as well as transverse mass spectra are presented in the third part.
At present we have calculated results imposing a hadron gas EoS to provide a baseline calculation to disentangle the effects of the different assumptions for the underlying dynamics in a transport vs. hydrodynamic calculation. The purely hadronic calculations can be compared in the broad energy regime from GeV where a vast amount of experimental data from BNLAGS and CERNSPS exists and which will be explored in more detailed energy scans by the FAIR project near GSI and the critRHIC program. Studies employing different EoS are delayed to future work to concentrate on effects of the underlying dynamics first.
Ii General aspects
The modelling of the dynamical evolution of heavy ion reactions is essential to gain further insights about the properties of the newly produced hot and dense QCD matter. Transport theory aims at the description of all stages of the collision on the basis of an effective solution of the relativistic Boltzmann equation (51)
(1) 
This equation describes the time evolution of the distribution functions for particle species and includes the full collision term on the right hand side. The interaction with external potentials leads to an additional term on the left hand side. The influence of potentials gets small at higher energies compared to the energy that is transferred by collisions. Therefore, they are dropped in Eqn. 1 and are not further discussed here. Usually, the collision kernel is truncated on the level of binary collisions and processes to keep the calculation numerically tractable. This microscopic approach has the advantage that it is applicable to nonequilibrium situations and the full phase space information is available at all stages of the heavy ion reaction. The restriction to binary collisions assumes large mean free paths of the particles. Between interactions the particle trajectories are given by straight line trajectories and particles are assumed to be in asymptotic states between the collisions (no “memory effect”).
This assumption might not be justified in the hot and very dense stage of heavy ion collisions anymore. In this regime the continuum limit in form of relativistic hydrodynamics might fit better to the characteristics of the system. The hydrodynamic evolution is governed by the energy and momentum conservation laws for given initial conditions, i.e. spatial distributions of energy and net baryon number densities. The coordinate space is divided into small cells in which the distribution functions correspond to equilibrium distributions (Fermi or Bose distribution). In this macroscopic approach the propagated quantities are net baryon number and energy densities which can be translated into information about the temperature and chemical potential via the specific equation of state (EoS). Since the evolution is driven by pressure gradients and the pressure is determined via the EoS, the EoS is the essential ingredient for the hydrodynamical evolution. Thus, hydrodynamics is a good tool to describe collective behaviour. Ideal hydrodynamics applies to systems with small mean free path, otherwise viscous effects have to be taken into account (52). A general advantage of hydrodynamics is the feature to explicitly incorporate phase transitions by changing the EoS.
However, in the late stage of the heavy ion reaction the system gets too dilute to apply ideal fluid dynamics. The hadronic rescatterings and decays of resonances have to be described, e.g., by using a transport description. Overall, there are two crucial points one has to take care of when building up a transport+hydrodynamics hybrid approach. The first is the initial switch from the microscopic to the macroscopic calculation where it has to be ensured that the local equilibrium assumption is fulfilled. The second one is the so called freezeout where the hydrodynamic fields are mapped to particles that are further propagated in a hadronic cascade. The freezeout transition must be placed in a region where both descriptions are valid at the same time,e.g., the phase transition region. In the following the specific implementation developed here will be discussed in more detail.
Iii Specific micro+macro approach
iii.1 UrQMD Approach
For our investigation, the UrQMD model (v2.3) (8); (9) is applied to heavy ion reactions from GeV. This nonequilibrium transport approach constitutes an effective solution of the relativistic Boltzmann equation (see Eqn. 1). The underlying degrees of freedom are hadrons and strings that are excited in high energetic binary collisions. Mean fields can in principle be taken into account in this framework, but the model is run in the so called cascade mode without interparticle potentials. To omit the potentials is reasonable because the inclusion of mean fields would not change the results in the energy range that we are considering here. Note that this is consistent with the calculation of the equation of state for the hydrodynamic evolution where no mean field has been taken into account as well.
The projectile and target nuclei are initialised according to a WoodsSaxon profile in coordinate space and Fermi momenta are assigned randomly for each nucleon in the rest frame of the corresponding nucleus. The hadrons are propagated on straight lines until the collision criterium is fulfilled. If the covariant relative distance between two particles gets smaller than a critical distance that is given by the corresponding total cross section a collision takes place,
(2) 
Each collision process is calculated in the rest frame of the binary collision. The reference frame that is used for the time ordering of the collisions and later on also for the switchings to and from the hydrodynamic phase is the equal speedsystem of the nucleusnucleus collision (for symmetric systems the equal speed system is identical to the center of mass system).
In UrQMD 55 baryon and 32 meson species, ground state particles and all resonances with masses up to GeV, are implemented with their specific properties and interaction cross sections. In addition, full particleantiparticle symmetry is applied. Isospin symmetry is assumed and only flavourSU(3) states are taken into account. The elementary cross sections are calculated by detailed balance or the additive quark model or are fitted and parametrized according to the available experimental data. For resonance excitations and decays the BreitWigner formalism, utilizing their vacuum properties is employed.
Towards higher energies, the treatment of subhadronic degrees of freedom is of major importance. In the present model, these degrees of freedom enter via the introduction of a formation time for hadrons produced in the fragmentation of strings (53); (54); (55). String excitation and fragmentation is treated according to the Lund model. For hard collisions with large momentum transfer ( GeV) PYTHIA is used for the calculation. A phase transition to a quarkgluon state is not incorporated explicitly into the model dynamics. However, a detailed analysis of the model in equilibrium yields an effective equation of state of Hagedorn type (56); (57); (58); (59); (60). The UrQMD transport model is successful in describing the yields and the spectra of various particles in protonproton, protonnucleus and nucleusnucleus collisions (61). A compilation of results of the actual version UrQMD2.3 compared to experimental data can be found in (62).
Apart from the success of transport simulations to describe spectra and yields certain problems remain:

Strangeness, especially multistrange baryons are not produced in sufficient amounts (67).
These observables that are sensitive to the early stage of the collision (pressure) or to the approach of thermal and chemical equilibrium during the collision history hint to the fact that a purely hadronic transport model may not be sufficient to describe the dynamics of the hot and dense stage of heavy ion reactions at higher energies (63); (68); (69); (64); (66). Therefore, these observations exemplify the need to embed a full threedimensional relativistic fluid dynamics description for these stages of the reaction.
For the results that are shown in this paper, the reference calculations are always performed employing the stateoftheart UrQMD2.3 model.
iii.2 Initial Conditions
The Ultrarelativistic Quantum Molecular Dynamics Model is used to calculate the initial state of a heavy ion collision for the hydrodynamical evolution (35). This is necessary to account for the nonequilibrium nature of the very early stage of the collision. Eventbyevent fluctuations of the initial state are naturally included by this setup. The coupling between the UrQMD initial state and the hydrodynamical evolution takes place when the two Lorentzcontracted nuclei have passed through each other. The initial time to begin with the hydrodynamical evolution is calculated via Eqn. 3 (and is assumed to be at least 1 fm/c):
(3) 
where R is the radius of the nucleus, is the nucleon mass and is the kinetic beam energy. This assures that (essentially) all initial baryonbaryon scatterings have proceeded and that the energy deposition has taken place. This is the earliest possible transition time where thermalization might be achieved (60). It is also convenient from the hydrodynamical point of view since at that time the two baryon currents that fly into oppposite directions have separated again.
In general, it is not wellestablished how and when chemical and kinetic equilibrium might have been reached in the early stage of the collision. One of the problems is, e.g., that the local equilibrium assumption might not apply equally well to all parts of the system at the same time in the computational frame which corresponds to the center of mass system of the two colliding nuclei. As a consequence, the faster particles have had less time in their local rest frame to equilibrate. For the bulk part and the high density region at midrapidity the difference between the two frames is small. These problems are present in all hydrodynamic/macroscopic approaches that rely on an equilibrium assumption and it is not our attempt to resolve these difficulties in this paper. One perspective might be the dynamical coupling between the initial transport calculation and the hydrodynamic evolution including source terms on both sides of the transition surface.
To allow for a consistent and numerically stable mapping of the ’point like’ particles from UrQMD to the 3dimensional spatialgrid with a cell size of , each hadron is represented by a Gaussian with a finite width. “Preformed” hadrons in the process of string fragmentation are also included in the transformation to the hydrodynamic quantities. I.e. each particle is described by a threedimensional Gaussian distribution of its total energy, momentum (in x, y, and zdirection) and baryon number density. The width of these Gaussians is chosen to be fm. A smaller Gaussian widths leads to numerical instabilities (e.g., entropy production) in the further hydrodynamical evolution, while a broader width would smear out the initial fluctuations to a large extent. To account for the Lorentzcontraction of the nuclei in the longitudinal direction, a Lorentzgammafactor is included. The resulting distribution function in the computational frame (cf), e.g., for the energy density, reads:
(4) 
where provides the proper normalisation, and are the energy density and total energy of the particle in the computational frame, while is the position vector of the particle. Summing over all single particle distribution functions leads to distributions of energy, momentum and baryon number densities in each cell.
To allow for calculations at finite impact parameter the spectators  nucleons that have not interacted at all before the start time of the hydrodynamic evolution  are propagated separately from the hydrodynamic evolution. The spectators are propagated on straight line trajectories in the usual cascade mode until the end of the hydrodynamic phase has been reached.
Instead of smearing out the initial distributions by describing the point like hadrons as Gaussian distributions, one could also obtain a smooth distribution by averaging over a large sample of UrQMD events. Our procedure of creating a new initial state for each event is motivated by the fact, that the experimental results all relate to observed (averaged) final, and not initial, states. Thus, eventbyevent fluctuations of the initial state can be observed, (e.g., in fluctuations) and have therefore been taken into account properly (for discussion of the importance of these fluctuations see, e.g., (36); (70)).
As an example, Figs. 1 and 2 show the energy and baryon number densities obtained in one single central ( fm) Pb+Pb collision at GeV after the initialisation of the hydrodynamic fields. The starting time is in this case fm and the densities in the figures correspond to the same time. All quantities are given in the local rest frame. The maximum values reach 6 GeV/fm for the energy density and around 8 times the nuclear ground state density for the baryon number density. The distributions are quite smooth which is necessary to provide possible initial conditions for the hydrodynamic evolution. One can see some peaks that correspond to local maxima of the distributions (“hot spots”). It is further clear that the single event distributions are not symmetric, neither in the transverse nor the longitudinal direction.
The remaining question is if the system is thermalized enough to assure that the local equilibrium assumption of ideal hydrodynamics is fulfilled. In our case, the hydrodynamic code transforms all the given quantities from the computational frame to the local rest frame of the energy momentum tensor which is also known as the Landau frame. This frame coincides in ideal hydrodynamics with the Eckart frame which is defined as the local rest frame of the baryon number current. The iterative calculation of the cell velocity succeeds if those two frames are close enough to each other. By this transformation the system is forced to local equilibrium.
iii.3 Hydrodynamic Evolution
Ideal relativistic one fluid dynamics is based on the conservation of energy, momentum and the net baryon number current. For the hydrodynamical evolution local equilibrium is assumed and zero viscosity which corresponds to zero mean free path. The two conservation equations that govern the evolution are (52); (71)
(5) 
where is the energymomentum tensor and is the baryon current. For an ideal fluid the energymomentum tensor and the net baryon number current take the simple form
(6) 
where and are the local rest frame energy density, pressure and net baryon density, respectively. is the four velocity of the cell and is the metric tensor. The local rest frame is defined as the frame where has diagonal form, (i.e. all offdiagonal elements vanish). The fourvelocity of the cells is calculated via the transformation into the local rest frame.
The equations of motion are solved in the following form by employing computational frame quantities and for the energy, momentum and net baryon number densities.
(7)  
(8)  
(9) 
In our case, the full (3+1) dimensional hydrodynamic evolution is performed using the SHASTA algorithm (72); (73). The partial differential equations are solved on a threedimensional spatial Eulerian grid with fixed position and size in the computational frame. The standard size of the grid is 200 cells in each direction while the cell size has been chosen to be which leads to timesteps of fm. Depending on the beam energy, the cell sizes may require adjustment to assure a stable solution of the differential equation.
The equation of state is needed as an additional input to calculate the pressure, temperature and chemical potential corresponding to the energy and the baryon number densities. Since the evolution of the system is driven by pressure gradients the EoS has the most important influence on the evolution.
iii.4 Equation of State
To solve the hydrodynamical equations, the EoS, the pressure as a function of energy and netbaryon number density, is needed as an input. Since the actual EoS of hot and dense QCD matter is still not precisely known, it may seem disadvantageous to have this additional uncertainty in the model. On the contrary it may prove to be an important trait of the model to be able to study changes on the dynamics of the bulk matter when changing the EoS thus finding observables for a phase transition in hot QCD matter. For recent discussions of different EoS and how to obtain EoS from lattice calculations, the reader is referred to (74); (75); (76).
The EoS used in the present calculations is a grand canonical description of a free, non interacting gas of hadrons. We will refer to it as the Hadron Gas (HG).
It follows from the hadronic chiral model presented in (77); (78). The chiral hadronic Lagrangian incorporates the complete set of
baryons from the lowest flavour octet, as well as
the entire multiplets of scalar, pseudoscalar, vector and axialvector
mesons (77). In meanfield approximation, the
expectation values of the scalar fields relevant for symmetric nuclear
matter correspond to the nonstrange and strange chiral quark
condensates, namely the and its counterpart
, respectively, and further the and vector
meson fields. Another scalar isoscalar field, the dilaton , is
introduced to model the QCD scale anomaly. However, if
does not couple strongly to baryonic degrees of freedom then it remains
essentially “frozen” below the chiral transition (77).
Consequently, we focus here on the role of the
quark condensates.
Interactions between baryons and scalar (BM) or vector (BV) mesons, respectively, are introduced as
(10)  
(11) 
Here, sums over the baryon octet (, , , ). A term with mass terms and quartic selfinteraction of the vector mesons is also added:
The scalar selfinteractions are
(12)  
Interactions between the scalar mesons induce the spontaneous breaking of chiral symmetry (first line) and the scale breaking via the dilaton field (last two terms).
Nonzero current quark masses break chiral symmetry explicitly in QCD. In the effective Lagrangian this corresponds to terms such as
(13) 
According to (10), the effective
masses of the baryons,
,
are generated through their coupling to the chiral condensates,
which attain nonzero vacuum expectation values due to their
selfinteractions (77) in (12).
The effective masses of the mesons are obtained as the second
derivatives of the mesonic potential
about its minimum.
All parameters of the chiral model discussed so far are fixed by either
symmetry relations, hadronic vacuum observables or nuclear matter
saturation properties (for details see (77)).
In addition, the model also provides a satisfactory description of realistic
(finitesize and isospin asymmetric) nuclei and of neutron stars
(77); (79); (80).
If the baryonic degrees of freedom are restricted to the members
of the lowest lying octet,
the model exhibits a smooth decrease of the chiral condensates
(crossover) for both high and high (77); (81).
However, additional baryonic degrees of freedom changes this into a firstorder phase transition
in certain regimes of the  plane, depending on the
couplings (82); (81); (83); (78); (35).
In what follows, the meson fields are replaced by their (classical)
expectation values, which corresponds to neglecting quantum and
thermal fluctuations. Fermions have to be integrated out
to oneloop. The grand canonical potential can then be written as
(14)  
where denote the baryonic and mesonic spinisospin degeneracy factors and are the corresponding single particle energies. The effective baryonchemical potentials are , with . Here is the quark, and the strange quark, chemical potential. The potentials of the mesons are given by the sum of the corresponding quark and antiquark chemical potentials. The vacuum energy (the potential at ) has been subtracted.
By extremizing one obtains selfconsistent gap equations
for the meson fields. Here, globally nonstrange matter is considered and
for any given and is adjusted to obtain a vanishing net strangeness.
Once the grand canonic potential is known as a function of and , all other thermodynamic quantities are derived straightforward. In its minima the grand canonic potential corresponds to , the pressure. The entropy density , number density and energy density then follow from the Euler relation,
(15) 
where the sum runs over all included hadron species.
Setting all hadron masses and chemical potentials to their vacuum values, and adding all reliably known heavy resonance states  with masses up to 2 GeV (84)  as free particles into (14), yields the above mentioned Hadron Gas EoS (85).
Hence, the hadronic degrees of freedom included in this EoS are consistent with the active degrees of freedom in the UrQMD model. This enables us to directly compare the dynamics of the hydrodynamic model with the transport simulation.
Using the chiral model and adding additional baryonic degrees of freedom as well as adjusting their scalar and vector coupling, an EoS with a phase structure including a first order phase transition and even a critical endpoint at finite can be obtained (78). This chiral EoS has already been successfully applied to a hydrodynamic calculation (35). Here the essentially different equation of state leads to distinguishable different results on the properties of bulk matter.
To emphasize these differences even more, a MITBag model EoS matched to an interacting hadron gas (73), generating a phase structure with a broad first order phase transition at all , can also be applied in our model.
Consequently these results will be compared to our present calculations, constituting observables for a phase transition in hot and dense QCD matter, even suggesting the order of this phase transition. Comparisons between the different EoS (i.e. free Hadron Gas, chiral EoS, and Bag model EoS) will be presented in a followup paper.
iii.5 Freezeout
Presently, the hydrodynamic evolution is stopped, if the energy density drops below five times the nuclear ground state energy density (i.e. MeV/) in all cells. This criterium corresponds to a Tconfiguration where the phase transition is expected (see dotted line in Fig. 4), i.e. a region where the hydrodynamic and the transport description are valid at the same time. The hydrodynamic fields are mapped to hadrons according to the CooperFrye equation (86)
(16) 
where are the boosted Fermi or Bose distributions corresponding to the respective particle species. Since we are dealing with an isochronous freezeout, the normal vector on the hypersurface is .
Let us note that it is of utmost importance to consider the same degrees of freedom on both sides of the hypersurface because otherwise energy and momentum conservation is violated. In our case, this is assured by the inclusion of the same particle species in the equation of state for the hydrodynamic calculation as in the transport calculation. In principle, it might also happen that particles are moving back into the hydrodynamic phase, however, the explosive character of heavy ion reactions, i.e. the rapid expansion flow suppresses the backstreaming effect. Therefore, this effect is negligible in our situation (87).
The assumption of an isochronous freezeout leads to fluctuations of the temperature and baryochemical potential distributions and not to single values for some of the thermodynamic quantities on the hypersurface. The quality of this approach and that the two (hydrodynamics and transport) prescriptions are valid through the applied range of switching temperatures is shown in the next section, where parameter tests and the time evolutions of the particle yields are explored in detail.
Fig. 3 shows the distribution of energy in the cells at freezeout with respect to temperature and baryochemical potential at GeV. We present here the energy distribution and not the number of cells because it is more interesting to know where the energy of the system sits than considering the many almost empty cells that do essentially not contribute to particle production. From Fig. 3 one obtains the mean values of the distributions that are in line with results from statistical model fits. The mean values of, e.g., the temperature can be calculated as
(17) 
where are the cell indices and the sum runs over all cells of one event. The net baryon number density has been used as a weighting factor.
To illuminate this finding in more detail, Fig. 4 shows the mean values of the temperature and baryochemical potential distributions at different energies in the plane for central ( fm) Au+Au/Pb+Pb collisions. Also, the widths of the distributions are depicted as “error” bars. Fig. 4 shows that the present freezeout distributions are similar to the parametrized curve for chemical freezeout as calculated by Cleymans et al. from statistical model fits to final particle multiplicities. The calculation by Dumitru et al. shows mean values as well as widths of temperature and baryochemical potential distributions that have been obtained by statistical model fits to final particle yields employing the assumptions of an inhomogeneous freezeout hypersurface. This calculation also leads to similar values as our calculation.
The effect that the mean temperature at the transition to the transport prescription saturates or even drops down a little at higher beam energies is related to the rapidity distribution of the temperature in the hydrodynamic cells at freezeout which is shown in Fig. 5 for three different beam energies. At low beam energies the midrapidity region coincides with the hottest region at freezeout. At higher SPS energies the situation changes. The hottest cells are at high rapidities while the midrapidity region has already cooled down well below the temperature of 170 MeV. This problem might be resolved by a different freezeout prescription on another hypersurface (e.g.,isotherm, iso) and is subject to future investigations. The best solution will be the dynamical coupling between hydrodynamics and transport which allows also for backstreaming contributions.
In the following the practical implementation will be explained in more detail. The implementation is based on a Monte Carlo sampling of Eqn. 16 and follows the general steps:

The particle numbers are calculated according to the following formula,
(18) where the index runs over the different particle species like, e.g., , p, or . is the boost factor between the computational frame and the cell. is the volume of the cell in the computational frame and is the particle number density. All cells with temperatures that are lower than 3 MeV are discarded from the following procedure because of numerical reasons. The local rest frame equilibrium distribution function is denoted by . To simplify the calculation, a relativistic Boltzmann distribution is used for all particles, except pions. It has been checked, that the Boltzmann approximation is sufficient to describe all particle species it is applied to. For the Boltzmann distribution the momentum integration leads to the following result for the particle number density
(19) where is the degeneracy factor for the respective particle species, is the mass of the particle to be produced, the temperature of the cell and is the modified Bessel function. The chemical potential includes the baryochemical potential and the strangeness chemical potential in the following way
(20) where is the quantum number for strangeness and is the baryon number.
For pions the Bose distribution has to be taken into account because the pion mass is on the order of the temperature of the system. In this case, the momentum integration involves an infinite sum over modified Bessel functions
(21) To calculate the number of particles in the computational frame the particle number density has to be multiplied with the Lorentzstretched volume of the cell ().

The average total number of particles in the cell, , is the sum over all particle numbers
(22) 
The total number of particles emitted from a cell, , is obtained from a Poisson distribution according to .
In the limit of small mean values, the Poisson distribution becomes . Thus it can be decided by one random number between 0 and 1 if a particle is produced in the respective cell. If the random number is smaller than one particle is produced and there is no particle production otherwise. The full Poisson distribution is used, if the particle number is larger than (). This assures an accuracy better than 1 %.

The particle type is chosen according to the probabilities .

The component of the isospin is distributed randomly because UrQMD assumes full isospin symmetry. To conserve the overall charge of the system and the initial isospinasymmetry the probability to generate the isospin component that leads to the right value of the charge that should be obtained in the end is favoured. The other isospin components are exponentially suppressed. The power of the exponential is proportional to the difference of the total charge generated by this produced particle and the required value.

The 4momenta of the particles are generated according to the CooperFrye equation (see Eqn. 16). For baryons and strange mesons the chemical potentials for baryon number and strangeness are taken into account.

The particle vector information is transferred back into the UrQMD model. The subsequent hadronic cascade calculation incorporates important final state effects as, e.g., rescatterings of the particles and resonance decays.
The above mentioned steps are pursued on random cells until the initial net baryon number is reached. Strangeness and charge are also conserved in each event separately, energy conservation is fulfilled for the mean values averaged over several events. Aiming at a realistic description of heavy ion reactions we perform the freezeout for each event separately and do not average over many freezeouts for one hydrodynamical evolution.
Iv Parameter tests
iv.1 Starttime and FreezeOut Criterium
In this section we investigate the dependences of observables on parameters of the implementation. Two important parameters have to be determined. The first one is the starting time which defines the initial switch from UrQMD to the hydrodynamic evolution. The second parameter is the freezeout criterium which is parametrized as an energy density criterium. While varying one parameter we have fixed the other one to the default value ( or ).
Fig. 6 (top) shows calculations of the total, i.e. pion and kaon () multiplicities, for four different starting times at two beam energies. The open symbols depict always the result at the highest AGS energy (GeV) and the filled symbols are the results at the SPS energy (GeV). The starting time is varied in factors of the default value that has been calculated via Eqn. 3. Displayed are results from halved to doubled initial time. One observes a higher pion production for earlier starting times compared to the pion production in the standard setup (1 ). This may be explained by the fact that the system is forced more strongly to equilibrium and the cascade evolution lasts longer. If the hydrodynamic evolution is started at later times (1.5 or 2 ) the resulting pion multiplicities are not affected anymore. The kaon yield is essentially not sensitive to the switching time. To summarize, varying the starting time by a factor of 4 results in a change in the pion and kaon production of less than compared to the pion and kaon production in the default configuration (1 ). In Fig. 6 (bottom) the mean transverse mass of pions and kaons at midrapidity is shown. The mean transverse mass values are calculated for the same four different starting times at the two exemplary beam energies as before. Here the results do also not change more than for a large spectrum of starting times. Therefore, our choice of the starting time as the geometrical criterium when the nuclei have passed through each other is sensible and stable. It is the earliest possible time where thermalization may have been achieved and the baryon currents have disconnected.
Fig. 7 (top) is the equivalent picture to Fig. 6 (top), however displaying the dependence of the total pion and kaon multiplicities on the freezeout criterium. The default value for the transition energy density is while we have varied it from . The higher the freezeout energy density the earlier the hydrodynamic evolution is stopped because the cells reach this critical energy density value earlier. As a consequence the kaon yields raises with an increase of the energy density criterium, while the pions remain virtually unchanged for all investigated transition criteria. Fig. 7 (bottom) shows the results for the pions and kaons mean transverse mass for two beam energies and different freezeout criteria. Again one observes only a very weak dependence on the freezeout criterion. Furthermore, the mean transverse mass values for the two different meson species are very similar since they acquire the same transverse flow. These findings confirm that our choice for the freezeout criterium as is robust.
iv.2 Timescales
In this Section the time scales that are important will be explored. Let us start with a table that summarizes the mean durations of the hydrodynamic and the hadronic phase of the collision for different starting times and freezeout parameters at GeV.
[fm]  [fm]  

1  4  7.68  15.63 
1  5  7.72  16.07 
1  6  6.84  16.49 
1  10  4.60  17.29 
0.5  5  7.03  14.59 
1.5  5  6.22  17.50 
2  5  4.91  17.61 
The duration of the hydrodynamic evolution is a welldefined period for each event because of the isochronous freezeout. The average is therefore an average over 100 events. The hadronic stage starts when the hydrodynamic evolution is stopped and it ends when the particles have undergone their last interaction. An interaction can be an inelastic or an elastic collision or a decay.
The averaged time duration of the stages of the reaction (given in Table 1) reflect the expectations. The lower the freezeout energy density the later the hydrodynamic freezeout proceeds and therefore the hydrodynamic evolution lasts longer while the hadronic stage is shortened. The later the hydrodynamic evolution starts, the bigger is, the shorter the hydrodynamic evolution lasts. The hadronic phase does not show a clear trend. To first approximation the final UrQMD stage lasts for fm independent of the parameters.
Fig. 8 shows the beam energy dependences of the timescales for the chosen values of and the freezeout energy density, from GeV. The starting time decreases as a function of beam energy from more than 10 fm at low energies to less than 2 fm at higher SPS energies. The mean duration of the hydrodynamic as well as the hadronic phase of the reaction grow with raising energy.
iv.3 Time Evolutions
In the following Section we investigate the time evolution of different quantities and compare the results of the hybrid model calculation to the UrQMD calculation without an hydrodynamic stage. Since the net baryon density is directly propagated on the hydrodynamic grid, it serves as a good example to compare to the default UrQMD calculation. In the microscopic approach the local rest frame density is calculated as the zero component of the net baryon number current in the frame where the corresponding local velocity vanishes (91).
Fig. 9 shows the time evolution of the net baryon number density at the center of a central Pb+Pb collision at GeV. The blue full line depicts the default UrQMD calculation while the red full line depicts the result of the hybrid model calculation. The result has some spikes because here we compare single events. There are two important observations: The absolute values of the net baryon number densities are very similar in both cases and there are no obvious discontinuities at the switching points to and from the hydrodynamic model calculation. The smoothness of the curve confirms our choice of parameters.
In Fig. 10 (top and middle) the time evolution of the particle yields in the two different models for the most central Pb+Pb collisions at GeV is compared. The multiplicities at different timesteps are extracted from the hydrodynamic evolution by converting the number densities to particle numbers via the freezeout procedure (see Section III.5). Fig. 10 (top) depicts the total particle multiplicity (red circles and full line) and the midrapidity multiplicity (blue squares and full line). The full lines indicate the default UrQMD calculation, while the symbols show the results of the hybrid model. The multiplicities increase rapidly in the initial 3 fm/c and then decrease a little, followed by a slower constant raise until the final multiplicity is reached. This qualitative behaviour is very similar in both approaches. The decrease of the multiplicity can be associated with the thermalization because absorption and production processes are on the same order. Note again that there are no discontinuities at the switching times in the hybrid calculation.
Next, we explore the time evolution for two particle species in more detail. In Fig 10 (middle), the pions (red circles and full line) represent the newly produced particles while the nucleons (blue squares and full line) are already there in the beginning as the constituents of the two incoming nuclei. The qualitative behaviour of the temporal evolution of the pion yield is similar to that discussed above for the total multiplicity. The decrease at the starting time of the hydrodynamic evolution fm/c is much stronger than in the model without hydrodynamic phase, because of the instant thermalization at the transition time. The default UrQMD transport calculation results in a similar, but much smoother, shape of the curve. This similarity hints to the fact that the microscopic calculation also equilibrates the hot and dense matter to a rather large degree. The number of nucleons decreases in the beginning due to the production of resonances and string excitations. At the thermalization the minimum is reached and the number of nucleons increases slowly until the final value is reached. In this case, not only the qualitative behaviour is independent of the underlying dynamics but also the absolute values are very close to each other.
To check quantum number and energy conservation the temporal evolution of all the important quantities in both approaches for central Pb+Pb collisions at GeV is depicted in Fig. 10 (bottom). The default UrQMD calculations are indicated by full lines while the hybrid model calculations with integrated hydrodynamic evolution are depicted by symbols. The net baryon number (orange triangles and full line), the charge (blue squares and full line) and the net strangeness (green diamonds and full line) are exactly conserved in both approaches. The total energy (red circles and full line) is only on average over several events conserved in the hybrid model calculation due to the freezeout prescription. But the fluctuations are on a level. Note however, that the total strangeness in the system ( quarks) is very different in both approaches. In the default transport calculation (black line) the total strangeness increases in the early stage of the collision and remains constant. This is contrasted by the hybrid calculation (dots). Due to the local thermal equilibration and the thermal production of strange particles in the hybrid calculation the yield of strange quarks jumps to a higher value at the switching time (). The total strangeness then decreases as the system cools down, but the final value remains higher than in the default transport calculation.
iv.4 Final State Interactions
The last step is the analysis of the freezeout process, i.e. how much hadronic interactions happen after the hydrodynamic evolution. For this purpose, we have calculated the number of collisions during the hadronic cascade calculation in dependence of time and of the elementary collisions. The corresponding distributions are shown in Figs. 11 and 12, respectively, for three different beam energies (11, 40 and 160 GeV, from top to bottom ).
Fig. 11 shows the collision rates for mesonmeson, mesonbaryon and baryonbaryon interactions. The grey area indicates the average time span of the hydrodynamic phase. One observes that substantial collision rates are present directly after the transition to UrQMD. The collision rates stay high for 5 fm/c, and after 3040 fm/c the system is completely frozen out. Only some resonance decays proceed for longer. According to the composition of the system the baryonbaryon and baryonmeson interactions dominate the lower beam energy result, while at the highest SPS energy the mesonmeson together with the mesonbaryon interactions are the most abundant type of collision indicating the transition from baryon dominated to meson dominated systems. Note that the overlap of the hadronic interaction phase with the hydrodynamic evolution results from the fact that the duration of the different stages fluctuates in the present approach. Shown here is the average duration of the hydrodynamic evolution.
Fig. 12 shows the distribution of the elementary collision in the freezeout process. One nicely observes all the resonance peaks in the corresponding channels. This figure suggests that the most abundant mesonbaryon collisions are excitations of the resonance (i.e. interactions) since there is a sharp peak at the mass ( GeV). For mesonmeson reactions the and the peaks are clearly present. This result indicates that there is still resonance regeneration even at this late stage of the system’s evolution.
V Results
v.1 Multiplicities and Particle Spectra
We start with a comparison of the multiplicities and particle spectra in the two frameworks. Calculations with the embedded hydrodynamic evolution employing a hadron gas equation of state for the high density stage of the collisions are compared to the reference results of the default transport calculation (UrQMD2.3). Since both calculations use the same initial conditions and freezeout prescription it allows to extract, which observables are sensitive to the change in the underlying dynamics, thus allowing to explore the effect of local equilibration, of viscosities and heat conductivity. The hybrid model calculation is in the following always depicted as full lines while the default UrQMD2.3 calculations are depicted as dotted lines. Note that we do not tune any parameters for different energies or centralities for the hybrid model calculation. The starting time for the hydrodynamic expansion is always calculated using Eqn. 3 and the fixed energy density criterium () for the freezeout (as explained in Section III.5) is always employed.
In Fig. 13 the excitation functions of the total multiplicities are shown for central Au+Au/Pb+Pb collisions for GeV. The present hybrid approach simulations have been restricted to this energy range because for calculations at higher energies some numerical subtleties have to be resolved, e.g. a dynamical grid size for the hydrodynamical evolution. Compared to the default simulation, the pion and proton multiplicities are decreased over the whole energy range in the hybrid model calculation due to the conservation of entropy in the ideal hydrodynamic evolution. The nonequilibrium transport calculation produces entropy and therefore the yields of nonstrange particles are higher. The production of strange particles however, is enhanced due to the establishment of full local thermal equilibrium in the hybrid model calculation. Since the abundance of strange particles is relatively small they survive the interactions in the UrQMD evolution that follows the hydrodynamic freezeout almost without rethermalization.
Fig. 14 shows the midrapidity yields of protons, pions, ’s and kaons as a function of the beam energy for central Au+Au/Pb+Pb collisions at GeV. For pions, kaons and ’s the same trend as for the multiplicities is observed. There are less pions produced in the hybrid model calculation due to entropy conservation, but more strange particles because of the production according to the local thermal equilibrium distributions. The proton yield at midrapidity is very similar in both calculations while there are less antiprotons produced in the hybrid calculation. Also for the strange antiparticles a reduction of the midrapidity yield in the hybrid calculation at higher SPS energies can be seen.
To explore the kinetics of the system in more detail, Fig. 15 shows the rapidity distribution for at three different energies ( and GeV). The general shape of the distribution is very similar in both approaches and in line with the experimental data. At higher energies even the absolute yields become very close to each other in both approaches.
Fig. 16 shows the rapidity distributions. In this case, the yield is higher in the hybrid calculation as already discussed above and also the shape of the distribution fits very nicely to the experimental data at SPS energies. Overall the rapidity distributions seem not to be too sensitive to the details of the dynamics for the hot and dense stage, but strangeness yields are influenced by the local equilibrium assumption. It seems that the local equilibrium assumption provides similar strangeness enhancement as previous calculations including additional strong color fields (67); (96). It is remarkable how well the hybrid calculation matches the rapidity spectra at lower energies (GeV), even though the transport calculation provides a slightly better description to the experimental data at this energy. One might still conclude from the rapidity spectra that the local equilibrium is not a good assumption for AGS energies.
v.2 Transverse Dynamics
After the longitudinal dynamics which reflects more the stopping power in the initial state we turn now to the transverse dynamics of the system. Transverse spectra are a promising candidate to be sensitive to the change in the underlying dynamics because they emerge from the transverse expansion which is mostly dominated by the evolution in the hot and dense stage of the reaction.
Figs. 17, 18 and 19 display the transverse mass spectra for pions, protons and kaons at midrapidity for central Au+Au/Pb+Pb reactions at three different beam energies. At GeV (Fig. 17) the differential transverse mass spectra are very similar for both calculations and are in line with the experimental data.
At GeV (Fig. 18) first differences become visible. Most notably is the strong flow of protons in the hybrid approach, that results in an overestimate of protons at high transverse momenta.
At the highest SPS energy (GeV, Fig. 19) all the transverse mass spectra are flatter in the hybrid approach. The initial pressure gradients are higher in the hydrodynamic calculation due to the hadronic equation of state without phase transition. Therefore, the matter expands faster in the transverse plane and higher transverse masses are reached. At this energy either the introduction of a mixed phase (first order phase transition) or nonequilibrium effects are necessary to explain the experimental data.
Fig. 20 shows the mean transverse mass excitation function for pions. It confirms the observations from the differential spectra. Up to GeV beam energy the hybrid model calculation leads to similar results as the default UrQMD calculation and is in line with the experimental data. The mean value of the transverse mass of pions is proportional to the temperature of the system and very different in the two calculations at higher energies. The UrQMD approach shows a softening of the equation of state in the region where the phase transition is expected because of nonequilibrium effects, while the hadron gas hydrodynamic calculation continuously rises as a function of the energy. This behaviour is well known, see, e.g., (113).
Finally Fig. 21 shows the mean transverse momenta as a function of particle mass for and particles. Here we observe the behaviour known from previous hybrid studies, that with increased strangeness, baryons accumulate less flow than in a complete hydrodynamic approach. This effect can be traced back to the small cross sections of multi strange baryons in the hadronic cascade, thus showing, that the freezeout/decoupling process proceeds gradually.
Vi Summary
We have presented the first fully integrated Boltzmann+hydrodynamics approach to relativistic heavy ion reactions. This hybrid approach is based on the Ultrarelativistic Quantum Molecular Dynamics (UrQMD) transport approach with an intermediate hydrodynamical evolution for the hot and dense stage of the collision. The specific coupling procedure including the initial conditions and the freezeout prescription have been explained. The eventbyevent character of the hybrid approach has been emphasized. The present implementation allows to compare pure microscopic transport calculations with hydrodynamic calculations using exactly the same initial conditions and freezeout procedure.
The parameter dependences of the model have been investigated and the time evolution of different quantities were explored. These tests led to the conclusion that the choice of the starting time and the freezeout criterium does generally alter the multiplicities and transverse mass spectra only on a level. The time evolution has shown that there are no discontinuities at the switching times in the hybrid model calculation. The importance of the final state interactions has been emphasized by demonstrating that there is still resonance regeneration after the hydrodynamic evolution.
The effects of the change in the underlying dynamics  ideal fluid dynamics vs. nonequilibrium transport theory  have been explored. The final pion and proton multiplicities are lower in the hybrid model calculation due to the isentropic hydrodynamic expansion while the yields for strange particles are enhanced due to the local equilibrium in the hydrodynamic evolution. The results of the different calculations for the mean transverse mass excitation function, rapidity and transverse mass spectra for different particle species at three different beam energies have been discussed in the context of the available data. The transverse expansion of the system is much faster in the hybrid model calculation, especially at higher energies which leads to differences in the observables that are sensitive to the transverse dynamics. This finding indicates qualitatively that “new” physical effects like, e.g., nonequilibrium effects or a phase transition have to be taken into account.
Forthcoming work will be devoted to the study of different equations of state and the effect of changes in the equation of state on observables. Also in progress are calculations within the hybrid model at higher beam energies (RHIC and LHC), but therefore specific numerical subtleties have to be resolved like, e.g., a dynamical grid size for the hydrodynamical evolution.
Acknowledgements.
We are grateful to the Center for Scientific Computing (CSC) at Frankfurt for providing the computing resources. The authors thank Dirk Rischke for providing the 1fluid hydrodynamics code. H. Petersen gratefully acknowledges financial support by the Deutsche TelekomStiftung and support from the Helmholtz Research School on Quark Matter Studies. This work was supported by GSI and BMBF.References
 J. W. Harris and B. Muller, Ann. Rev. Nucl. Part. Sci. 46, 71 (1996)
 S. A. Bass, M. Gyulassy, H. Stoecker and W. Greiner, J. Phys. G 25, R1 (1999)
 J. Adams et al. [STAR Collaboration], Nucl. Phys. A 757, 102 (2005)
 B. B. Back et al., Nucl. Phys. A 757, 28 (2005)
 I. Arsene et al. [BRAHMS Collaboration], Nucl. Phys. A 757, 1 (2005)
 K. Adcox et al. [PHENIX Collaboration], Nucl. Phys. A 757, 184 (2005)
 C. Alt et al. [NA49 Collaboration], Phys. Rev. C 77, 024903 (2008)
 S. A. Bass et al., Prog. Part. Nucl. Phys. 41, 255 (1998) [Prog. Part. Nucl. Phys. 41, 225 (1998)]
 M. Bleicher et al., J. Phys. G 25, 1859 (1999)
 D. Molnar and P. Huovinen, Phys. Rev. Lett. 94, 012302 (2005)
 Z. Xu and C. Greiner, Phys. Rev. C 71, 064901 (2005)
 Z. W. Lin, C. M. Ko, B. A. Li, B. Zhang and S. Pal, Phys. Rev. C 72, 064901 (2005)
 G. Burau, J. Bleibel, C. Fuchs, A. Faessler, L. V. Bravina and E. E. Zabrodin, Phys. Rev. C 71, 054905 (2005)
 H. W. Barz and B. Kampfer, Nucl. Phys. A 683, 594 (2001)
 A. B. Larionov, O. Buss, K. Gallmeister and U. Mosel, Phys. Rev. C 76, 044909 (2007)
 J. Bleibel, G. Burau, A. Faessler and C. Fuchs, Phys. Rev. C 76, 024912 (2007)
 J. Bleibel, G. Burau and C. Fuchs, Phys. Lett. B 659, 520 (2008)
 B. Andersson, G. Gustafson and C. Peterson, Nucl. Phys. B 135, 273 (1978).
 J. R. Ellis and K. Geiger, Phys. Rev. D 52, 1500 (1995)
 T. S. Biro, P. Levai and J. Zimanyi, Phys. Rev. C 59, 1574 (1999)
 C. T. Traxler, U. Mosel and T. S. Biro, Phys. Rev. C 59, 1620 (1999)
 M. Hofmann, M. Bleicher, S. Scherer, L. Neise, H. Stoecker and W. Greiner, Phys. Lett. B 478, 161 (2000)
 T. Hirano, Phys. Rev. C 65, 011901 (2002)
 P. F. Kolb and U. W. Heinz, arXiv:nuclth/0305084.
 C. Nonaka and S. A. Bass, Phys. Rev. C 75, 014902 (2007)
 P. F. Kolb, P. Huovinen, U. W. Heinz and H. Heiselberg, Phys. Lett. B 500, 232 (2001)
 P. F. Kolb, U. W. Heinz, P. Huovinen, K. J. Eskola and K. Tuominen, Nucl. Phys. A 696, 197 (2001)
 P. Huovinen, P. F. Kolb, U. W. Heinz, P. V. Ruuskanen and S. A. Voloshin, Phys. Lett. B 503, 58 (2001)
 A. Dumitru, S. A. Bass, M. Bleicher, H. Stoecker and W. Greiner, Phys. Lett. B 460, 411 (1999)
 S. A. Bass, A. Dumitru, M. Bleicher, L. Bravina, E. Zabrodin, H. Stoecker and W. Greiner, Phys. Rev. C 60, 021902 (1999)
 S. A. Bass and A. Dumitru, Phys. Rev. C 61, 064909 (2000)
 S. Soff, S. A. Bass and A. Dumitru, Phys. Rev. Lett. 86, 3981 (2001)
 S. Soff, S. A. Bass, D. H. Hardtke and S. Y. Panitkin, Phys. Rev. Lett. 88, 072301 (2002)
 C. Nonaka and S. A. Bass, Nucl. Phys. A 774, 873 (2006)
 J. Steinheimer, M. Bleicher, H. Petersen, S. Schramm, H. Stocker and D. Zschiesche, Phys. Rev. C 77, 034901 (2008)
 S. Paiva, Y. Hama and T. Kodama, Phys. Rev. C 55, 1455 (1997).
 C. E. Aguiar, Y. Hama, T. Kodama and T. Osada, Nucl. Phys. A 698, 639 (2002)
 O. . J. Socolowski, F. Grassi, Y. Hama and T. Kodama, Phys. Rev. Lett. 93, 182301 (2004)
 F. Grassi, Y. Hama, O. Socolowski and T. Kodama, J. Phys. G 31, S1041 (2005).
 R. Andrade, F. Grassi, Y. Hama, T. Kodama, O. . J. Socolowski and B. Tavares, Eur. Phys. J. A 29, 23 (2006)
 R. Andrade, F. Grassi, Y. Hama, T. Kodama and O. . J. Socolowski, Phys. Rev. Lett. 97, 202302 (2006)
 C. E. Aguiar, T. Kodama, T. Koide and Y. Hama, Braz. J. Phys. 37, 95 (2007).
 V. V. Skokov and V. D. Toneev, Phys. Rev. C 73, 021902 (2006)
 D. Teaney, J. Lauret and E. V. Shuryak, arXiv:nuclth/0110037.
 T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and Y. Nara, Phys. Lett. B 636, 299 (2006)
 T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and Y. Nara, Phys. Rev. C 77, 044909 (2008)
 C. Anderlik et al., Phys. Rev. C 59, 3309 (1999)
 V. K. Magas et al., Heavy Ion Phys. 9, 193 (1999)
 K. A. Bugaev and M. I. Gorenstein, arXiv:nuclth/9903072.
 K. A. Bugaev, Phys. Rev. Lett. 90, 252301 (2003)
 S. R. De Groot, W. A. Van Leeuwen and C. G. Van Weert, Amsterdam, Netherlands: Northholland ( 1980) 417p
 L. D. Landau, Izv. Akad. Nauk Ser. Fiz. 17, 51 (1953).
 B. Andersson, G. Gustafson and B. NilssonAlmqvist, Nucl. Phys. B 281, 289 (1987).
 B. NilssonAlmqvist and E. Stenlund, Comput. Phys. Commun. 43, 387 (1987).
 T. Sjostrand, Comput. Phys. Commun. 82, 74 (1994).
 S. A. Bass et al., Phys. Rev. Lett. 81, 4092 (1998)
 L. V. Bravina et al., J. Phys. G 25, 351 (1999)
 M. Belkacem et al., Phys. Rev. C 58, 1727 (1998)
 L. V. Bravina et al., Phys. Rev. C 60, 024904 (1999)
 L. V. Bravina et al., arXiv:0804.1484 [hepph].
 E. L. Bratkovskaya et al., Phys. Rev. C 69, 054907 (2004)
 H. Petersen, M. Bleicher, S. A. Bass and H. Stocker, arXiv:0805.0567 [hepph].
 M. Bleicher and H. Stoecker, Phys. Lett. B 526, 309 (2002)
 H. Petersen, Q. Li, X. Zhu and M. Bleicher, Phys. Rev. C 74, 064908 (2006)
 S. Kniege et al. [NA49 Collaboration], AIP Conf. Proc. 828, 473 (2006)
 Q. Li, M. Bleicher and H. Stocker, Phys. Lett. B 659, 525 (2008)
 S. Soff et al., Phys. Lett. B 471, 89 (1999)
 D. Molnar and M. Gyulassy, Phys. Rev. Lett. 92, 052301 (2004)
 U. W. Heinz and P. F. Kolb, arXiv:hepph/0204061.
 H. J. Drescher and Y. Nara, Phys. Rev. C 75, 034905 (2007)
 R. B. Clare and D. Strottman, Phys. Rept. 141, 177 (1986).
 D. H. Rischke, S. Bernard and J. A. Maruhn, Nucl. Phys. A 595, 346 (1995)
 D. H. Rischke, Y. Pursun and J. A. Maruhn, Nucl. Phys. A 595, 383 (1995) [Erratumibid. A 596, 717 (1996)]
 M. Bluhm, B. Kampfer, R. Schulze, D. Seipt and U. Heinz, Phys. Rev. C 76, 034901 (2007)
 K. Redlich, B. Friman and C. Sasaki, J. Phys. G 35, 044013 (2008)
 T. S. Biro, P. Levai, P. Van and J. Zimanyi, J. Phys. G 32, S205 (2006)
 P. Papazoglou, D. Zschiesche, S. Schramm, J. SchaffnerBielich, H. Stoecker and W. Greiner, Phys. Rev. C 59, 411 (1999)
 D. Zschiesche, G. Zeeb and S. Schramm, J. Phys. G 34, 1665 (2007)
 S. Schramm, Phys. Lett. B 560, 164 (2003)
 S. Schramm, Phys. Rev. C 66, 064310 (2002)
 D. Zschiesche, S. Schramm, H. Stoecker and W. Greiner, Phys. Rev. C 65, 064902 (2002)
 J. Theis, G. Graebner, G. Buchwald, J. A. Maruhn, W. Greiner, H. Stoecker and J. Polonyi, Phys. Rev. D 28, 2286 (1983).
 D. Zschiesche, G. Zeeb, S. Schramm and H. Stoecker, J. Phys. G 31, 935 (2005)
 S. Eidelman et al. [Particle Data Group], Phys. Lett. B 592, 1 (2004).
 D. Zschiesche, S. Schramm, J. SchaffnerBielich, H. Stoecker and W. Greiner, Phys. Lett. B 547, 7 (2002)
 F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).
 K. A. Bugaev, Phys. Rev. C 70, 034903 (2004)
 J. Cleymans and K. Redlich, Phys. Rev. Lett. 81, 5284 (1998)
 J. Cleymans, H. Oeschler, K. Redlich and S. Wheaton, J. Phys. G 32, S165 (2006)
 A. Dumitru, L. Portugal and D. Zschiesche, Phys. Rev. C 73, 024902 (2006)
 S. Vogel, H. Petersen, J. Aichelin and M. Bleicher, arXiv:0710.4463 [hepph].
 J. L. Klay et al. [E0895 Collaboration], Phys. Rev. C 68, 054905 (2003)
 C. Pinkenburg et al. [E895 Collaboration], Nucl. Phys. A 698, 495 (2002)
 P. Chung et al. [E895 collaboration], Phys. Rev. Lett. 91, 202301 (2003)
 S. V. Afanasiev et al. [The NA49 Collaboration], Phys. Rev. C 66, 054902 (2002)
 M. Bleicher, W. Greiner, H. Stoecker and N. Xu, Phys. Rev. C 62, 061901 (2000)
 T. Anticic et al. [NA49 Collaboration], Phys. Rev. Lett. 93, 022302 (2004)
 A. Richard [NA49 Collaboration], J. Phys. G 31, S155 (2005).
 M. K. Mitrovski et al. [NA49 Collaboration], J. Phys. G 32, S43 (2006)
 C. Alt et al. [NA49 Collaboration], arXiv:0804.3770 [nuclex].
 C. Blume [NA49 Collaboration], J. Phys. G 31, S685 (2005)
 S. V. Afanasiev et al. [NA49 Collaboration], Phys. Lett. B 538, 275 (2002)
 L. Ahle et al. [E866 Collaboration], Phys. Lett. B 476, 1 (2000)
 K. Adcox et al. [PHENIX Collaboration], Phys. Rev. C 69, 024904 (2004)
 J. H. Lee et al. [BRAHMS Collaboration], J. Phys. G 30, S85 (2004).
 L. Ahle et al. [E866 Collaboration], Phys. Lett. B 490, 53 (2000)
 D. Ouerdane [BRAHMS Collaboration], Nucl. Phys. A 715, 478 (2003)
 S. Ahmad et al., Phys. Lett. B 382, 35 (1996).
 A. Mischke et al. [NA49 Collaboration], J. Phys. G 28, 1761 (2002)
 J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 92, 112301 (2004)
 Y. Akiba et al. [E802 Collaboration], Nucl. Phys. A 610, 139C (1996).
 C. Alt et al. [NA49 Collaboration], Phys. Rev. C 73, 044910 (2006).
 H. Sorge, Phys. Lett. B 402, 251 (1997)