Three-fluid Hydrodynamics-based Event Simulator Extended by UrQMD final State interactions (THESEUS) for FAIR-NICA-SPS-BES/RHIC energies

We present a new event generator based on the three-fluid hydrodynamics (3FH) approach, followed by a particlization at the hydrodynamic decoupling surface and a subsequent UrQMD afterburner stage based on the microscopic UrQMD transport model that accounts for hadronic final state interactions. First results for Au+Au collisions are presented. The following topics are addressed: the directed flow, transverse-mass spectra, as well as rapidity distributions of protons, pions and kaons for two model equations of state, one with a first-order phase transition, the other with a crossover transition. Preliminary results on the femtoscopy are also discussed. We analyze the accuracy of reproduction of the 3FH results by the new event generator and the effect of the subsequent UrQMD afterburner stage.


Introduction
The onset of deconfinement in relativistic heavy-ion collisions and the search for a critical endpoint is now in the focus of theoretical and experimental studies of the equation of state (EoS) and the phase diagram of strongly interacting matter. This challenge is one of the main motivations for the currently running beam-energy scan (BES) at the Relativistic Heavy-Ion Collider (RHIC) at Brookhaven National Laboratory (BNL) [1] and at the CERN Super-Proton-Synchrotron (SPS) [2] as well as for a speaker, e-mail: Y.Ivanov@gsi.de arXiv:1711.07959v1 [nucl-th] 21 Nov 2017 EPJ Web of Conferences constructing the Facility for Antiproton and Ion Research (FAIR) in Darmstadt [3] and the Nuclotronbased Ion Collider fAcility (NICA) in Dubna [4].
Three-fluid hydrodynamics (3FH) [5] was designed to simulate heavy-ion collisions at moderately relativistic energies, i.e. precisely in the energy range of the expected onset of deconfinement. In recent years applications of 3FH demonstrated a strong preference of deconfinement scenarios for the explanation of available experimental data [6][7][8][9][10][11][12][13]. However, up to now 3FH has been facing certain problems. From the theoretical side, the model lacks an afterburner stage that can play an important role for some observables. From the practical point of view, the model was not well suited for data simulations in terms of experimental events, because the model output consisted of fluid characteristics rather than of a set of observable particles.
In this contribution, we present a new Three-fluid Hydrodynamics-based Event Simulator Extended by UrQMD final State interactions (THESEUS) [14] and its application to the description of heavy-ion collisions in the FAIR-NICA-SPS-BES/RHIC energy range. This simulator provides a solution to both the above-mentioned problems. It presents the 3FH output in terms of a set of observed particles and the afterburner can be run starting from this output by means of the UrQMD model [15,16]. Thus THESEUS as a new tool allows to discuss the multifaceted physics challenges at FAIR-NICA-SPS-BES/RHIC energies. The new simulation program has the unique feature that it can describe a hadron-to-quark matter transition of the first order which proceeds in the baryon stopping regime that is not accessible to previous simulation programs designed for higher energies. Besides this, with THESEUS one can address practical questions like the influence of hadronic final state interactions and of the detector acceptance, which are necessary to focus on important physics questions. These deal with a potential discovery and investigation of the first-order phase transition line, where during a heavy-ion collision the EoS reaches its softest point [17]. It remains an open question how this characteristic feature of the EoS manifests itself in observables such as flow, proton rapidity distributions and femtoscopic radii. Particular emphasis is on the robustness of the "wiggle" [18] in the energy scan of the midrapidity curvature of the proton rapidity distribution that has been suggested as a possible signal for a first order phase transition, expected just in the range of energies at NICA and FAIR experiments.
At present THESEUS is not an integrated approach. The simulation proceeds in two steps: first, a numerical solution of the 3-fluid hydrodynamics is computed with the corresponding code. Based on the output of the hydrodynamic part, a Monte Carlo procedure is used to sample the ensemble of hadron distributions and the UrQMD code is engaged to calculate final state hadronic rescatterings, as will be explained below. Another present limitation which we leave for future work is the absence of event-by-event hydrodynamic evolution. Therefore later by an event we mean a Monte Carlo sampled set of final hadrons, which correspond to the same (average) hydrodynamic evolution.

The 3FH model
The 3FH model treats the collision process from the very beginning, i.e. from the stage of cold nuclei, up to the particlization from the fluid dynamics. This model is a straightforward extension of the two-fluid model with radiation of direct pions [19,20] and of the (2+1)-fluid model of Refs. [21,22]. The 3-fluid approximation is a minimal way to simulate the finite stopping power at the initial stage of the collision. Within the 3-fluid approximation a generally nonequilibrium distribution of baryonrich matter is simulated by counter-streaming baryon-rich fluids initially associated with constituent nucleons of the projectile (p) and target (t) nuclei. Therefore, the initial conditions for the fluid evolution are two Lorentz contracted spheres with radii of corresponding nuclei and zero diffuseness, ICNFP 2017 baryon density n B = 0.15 fm −3 and energy density m N n B 0.14 GeV/fm 3 . In addition, newly produced particles, populating the mid-rapidity region, are associated with a fireball (f) fluid. Each of these fluids is governed by conventional hydrodynamic equations. The continuity equations for the baryon charge read for α =p and t, where J µ α = n α u µ α is the baryon current defined in terms of proper (i.e. in the local rest frame) net-baryon density n α and hydrodynamic 4-velocity u µ α normalized as u αµ u µ α = 1. Eq. (1) implies that there is no baryon-charge exchange between p-, t-and f-fluids, as well as that the baryon current of the fireball fluid is identically zero, J µ f = 0, by construction. Equations of the energymomentum exchange between fluids are formulated in terms of energy-momentum tensors T µν α of the fluids where the F ν α are friction forces originating from inter-fluid interactions. F ν p and F ν t in Eqs.
(2)-(3) describe energy-momentum loss of the baryon-rich fluids due to their mutual friction. A part of this loss |F ν p − F ν t | is transformed into thermal excitation of these fluids, while another part (F ν p + F ν t ) gives rise to particle production into the fireball fluid (see Eq. (4)). F ν fp and F ν ft are associated with friction of the fireball fluid with the p-and t-fluids, respectively. Here τ f is the formation time, and is the 4-velocity of the free propagation of the produced fireball matter. In fact, this is the velocity of the fireball matter at the moment of its production. According to Eq. (4), this matter gets formed only after the time span U 0 F τ f upon the production, and in a different space point x − U F (x ) τ f , as compared to the production point x . The friction between fluids was fitted to reproduce the stopping power observed in proton rapidity distributions for each EoS, as it is described in Refs. [5,6] in detail.
Different equations of state (EoS) can be applied within the 3FH model. The recent series of simulations [6][7][8][9][10][11][12][13] was performed employing three different types of EoS: a purely hadronic EoS [23] (hadr. EoS) and two versions of the EoS involving deconfinement [24]. The latter two versions are an EoS with a first-order phase transition (2-phase EoS) and one with a smooth crossover transition (crossover EoS). Figure 1 illustrates the differences between the three considered EoS.
An application of the 3FH model is illustrated in Fig. 2 where the evolution of the proper (i.e. in the local rest frame) baryon density in the reaction plane is presented for a semi-central (impact parameter b = 6 fm) Au+Au collision at √ s NN = 6.4 GeV (E lab = 20 A GeV). The simulation was performed with the crossover EoS without freeze-out. As can be seen from that figure, very high baryon densities are reached in the central region of the colliding system. The freeze-out criterion used in the 3FH model is ε < ε frz , where ε is the total energy density of all three fluids in their common rest frame. More details can be found in Refs. [25,26] The freezeout energy density ε frz = 0.4 GeV/fm 3 was chosen mostly on the condition of the best reproduction of secondary particle yields for all considered scenarios, see [5]. An important feature of the 3FH freeze-out is an antibubble prescription, preventing the formation of bubbles of frozen-out matter inside the dense matter while it is still hydrodynamically evolving. The matter is allowed to be frozen EPJ Web of Conferences  out only if it is located near the border with the vacuum (this piece of matter gets locally frozen out). The thermodynamic quantities of the frozen-out matter are recalculated from the in-matter EoS, with which the hydrodynamic calculation runs, to the hadronic gas EoS 1 . This is done because a part of the energy is still accumulated in collective mean fields at the freeze-out instant. This mean-field energy needs to be released before entering the hadronic cascade in order to facilitate energy conservation.
The output of the model is recorded in terms of Lagrangian test particles (i.e. fluid droplets) for each fluid α (= p, t or f). Each particle contains information on space-time coordinates (t, x) of the frozen-out matter, proper volume of the test particle of the α fluid (V pr α ), hydrodynamic velocity (u µ α ) in the frame of computation, temperature (T α ), baryonic (µ Bα ) and strange (µ Sα ) chemical potentials.

Particlization
In the multi-fluid approach one simulates the heavy ion collision from its very first moment using fluid dynamics. However, once the system becomes too dilute, the fluid approximation loses its applicability and individual particles are the relevant degrees of freedom. The process of changing from a fluid to a particle description is called "particlization" [27]. Since we supplement the 3FH with a hadronic cascade, the particlization is not freeze-out anymore. By definition there are only resonance decays after freeze-out, whereas in the present generator final state hadronic rescattering processes are simulated as well using the UrQMD code.

ICNFP 2017
The particlization criterion is chosen to be the same as freeze-out criterion in [5], e.g., where ε tot is defined as: ε tot = T * 00 p + T * 00 t + T * 00 f and the asterisk denotes a reference frame where the nondiagonal components of the total energy momentum tensor are zero. This choice allows to study the influence of hadronic rescatterings to the observables by comparing them with the ones calculated in previous 3-fluid hydrodynamic models.
For the details of fluid to particle conversion the reader is referred to [5], whereas here we repeat the details important for the construction of the Monte Carlo sampling procedure. Both the baryonrich projectile and target fluids as well as the fireball fluid are being frozen out in small portions, therefore the output of the particlization procedure is a set of droplets (or surface elements). Each droplet is characterized by its proper volume V pr , temperature T , baryon, µ B , strange chemical potentials µ S , and the collective flow velocity u µ .
The thermodynamic parameters of the droplets correspond to a free hadron resonance gas. Therefore, we proceed with sampling the hadrons according to their phase space distributions (see Eq. (33) in [5]), which are expressed in the rest frame of the fluid element (FRF) as where the asterisk denotes momentum in the fluid rest frame, where u * µ α = (1, 0, 0, 0), µ αi = B i · µ αB + S i · µ αS is the chemical potential of hadron i with baryon number B i , strangeness S i , degeneracy factor g i , and the α summation runs over droplets from all (p, t and f) fluids.
The use of temperature and chemical potentials implies a grand canonical ensemble for each surface element. The sampling is therefore organized as a loop over all droplets, every iteration of which consists of the following steps [28,29] • average multiplicities of all hadron species are calculated according to together with their sum ∆N tot,α = i ∆N i,α ; • total (integer) number of hadrons from each surface element is sampled according to Poisson distribution with mean ∆N tot,α . If the number is greater than zero, sort of hadron is randomly chosen based on probabilities ∆N i,α /∆N tot,α ; • hadron's momentum in FRF p * is sampled according to (6), which is isotropic in momentum space; • momentum vector is Lorentz boosted to the global frame of the collision.
In the present version of the generator, also from the arguments of consistency with preceding hydrodynamic evolution, we do not apply any corrections over the grand canonical procedure to account for effects of charge or energy conservation. Therefore, particle multiplicities fluctuate from event to event according to the composition of grand canonical ensembles given by the individual droplets.

UrQMD simulation of final state interactions
The Ultra-relativistic Quantum Molecular Dynamics (UrQMD) approach [15,16] treats hadrons and resonances up to a mass of ∼ 2.2 GeV. All binary interactions are treated via the excitation and decay of resonances or string excitation and decay and elastic scatterings. It is crucial for a state-of-the-art event generator to treat the interactions during the late non-equilibrium hadronic stage of heavy ion reactions properly. At RHIC and LHC notable differences in the proton yields have been observed and the identified particle spectra and flow observables show an effect of the hadronic rescattering (for a review of hybrid approaches see [30]). At lower beam energies as they are investigated in this work, the hadronic stage of the reaction is of utmost importance. In [31] it has been shown, that the excitation function of elliptic and triangular flow can only be understood within a combined hydro-dynamics+transport approach. UrQMD constitutes an effective solution of the relativistic Boltzmann equation and therefore provides access to the full phase-space distribution of all individual particles at all times. In this work the effect of hadronic rescattering in the final state on the identified particle spectra and the rapidity dependent directed flow is demonstrated in detail.

Results
In this section we present a selection of first results from THESEUS for the energy scan ( √ s NN = 4 − 11 GeV) planned at the NICA-MPD collider experiment, which has overlap with the FAIR-SPS-BES/RHIC energy range.

Tests of the particlization routine: Spectra of pions, kaons and protons
We start by showing the transverse momentum distributions of pions [(π + + π 0 + π − )/3] and kaons [(K + + K 0 )/2] in Fig. 3. They are calculated from a sample of 30000 events generated according to the Monte Carlo procedure described above, and are compared in the plot to direct integration according to Eq. (6) and inclusion of resonance decay contributions performed in the basic 3FH part of THESEUS (the afterburner is turned off for this comparison). The 3FH evolution simulates Au+Au collisions at E lab = 30 A GeV with the two-phase EoS. We observe excellent agreement up to p T = 2.2 GeV, which is limited by the generated event statistics. In Fig. 4 we show the rapidity distributions for the same setup. The rapidity distributions reveal a small difference in kaon yields, and an even smaller one for pions, which is attributed to differences in the large mass sector of the resonance tables and branching ratios between 3FH and THESEUS. Nevertheless, the shapes of rapidity distributions agree beautifully. In Figs which are included in THESEUS. They lead to a slight steepening of the p T spectrum for pions and to a reduction of the double-peak structure in the kaon rapidity spectrum. Both are sufficiently gentle effects to not spoil our conclusion. The tests demonstrate that both the procedure of particle sampling at particlization and the resonance decay kinematics are implemented correctly.

Directed flow of protons and pions
Next we test whether more subtle features of particle distributions are preserved by the particlization procedure and how they are affected by the hadronic cascade. First we calculate the directed flow coefficient v 1 for pions and protons as a function of rapidity using the reaction plane method where Ψ RP = 0 in the model, since the impact parameter is always directed along x-axis. Although the generator makes it possible to apply different methods of flow analysis over generated events, we use the reaction plane method in order to perform a one-to-one comparison between results from THESEUS with and without UrQMD and the corresponding ones from the basic 3FH model. We present the results of THESEUS with and without UrQMD afterburner for the directed flow v 1 of protons and pions at E lab = 8 A GeV (Fig. 5) and 30 A GeV (Fig. 6), comparing the case of the 2-phase EoS (first order phase transition, upper panels) with that of the crossover EoS (lower panels) at central (left panels), semicentral (middle panels) and peripheral (right panels) Au+Au collisions. Dashed lines show the results from THESEUS without hadronic cascade, where we quantitatively reproduce the results from basic 3FH model. This figure shows the influence of hadronic final state interactions on the patterns of directed flow of protons and pions in the the broad rapidity range −1.5 < y < 1.5 and how it evolves from low energies to high energies. At E lab = 8 A GeV in Fig. 5 we observe that hadronic rescattering causes the transition from flow to antiflow for pions due to the shadowing by a dense baryonic medium. The flow of protons is not affected by the hadronic rescattering, which remains so for all energies. The shadowing effect on the pion directed flow becomes less important at higher energies. At 30 A GeV hadronic rescattering has no effect on the directed flow of pions in the central rapidity region.  This behaviour can be understood as follows. If there is only hydrodynamics, the pions are emitted along the fluid flow, while when there is rescattering they are blocked by the baryonic matter in the projectile and target region, therefore the anti-flow appears. This effect was first demonstrated in Ref. [15]. This effect of the pion shadowing is more spectacular in Fig. 5 where the directed flow of protons and pions at E lab = 8 A GeV is presented. As seen, the proton v 1 is practically insensitive to the UrQMD afterburner, while the pion v 1 is strongly affected by this afterburner. The afterburner even changes the pion v 1 flow to an antiflow. The effect of the pion shadowing becomes weaker with the collision energy rise, as it is seen at E lab = 30 A GeV, because the midrapidity region becomes ICNFP 2017 less baryon abundant. Though, at larger collision energies and peripheral rapidities, this shadowing is still noticeable. This shadowing results in better agreement with the STAR data on pion v 1 [33].
There is hardly any difference to be noticed in the pion directed flow patterns between the case of a 2-phase EoS and a crossover EoS. Figure 7. Energy scan for the curvature C y of the net proton rapidity distribution at midrapidity for central Au+Au collisions with impact parameter b = 2 fm. We compare the 3FH model result (black solid lines) with THESEUS (blue short-dashed lines) and THESEUS without UrQMD (red long-dashed lines). The results for the two-phase EoS (upper row) are compared to those for the crossover EoS (lower row). The "wiggle" as a characteristic feature for the EoS with a first order phase transition is rather robust against hadronic final state interactions. Data from AGS experiments are shown by filled squares, data from NA49 by filled triangles.

Baryon stopping signal for a first-order phase transition
In Fig. 7 we show the reduced curvature of the net proton rapidity distribution C y = y 2 cm (d 3 N net−p /dy 3 )/(dN net−p /dy), where y cm is the rapidity of the center of mass of the colliding system in the frame of the target [18,35,36]. Because of a narrower collision energy range chosen here, we observe only the peak-dip part of the so-called "peak-dip-peak-dip" structure reported in [18,35,36]. The reduced curvature is calculated by fitting the rapidity distribution with a 2 nd order polynomial of the form P 2 (y) = ay 2 + by + c for which then C y = y 2 beam 2a/c results. Contrary to the basic 3FH model which can calculate C y with any given precision, in the Monte Carlo procedure the accuracy depends on the event statistics and binning. A reliable determination of C y requires not less than 10 4 events for central and semi-central collisions and 10 5 events for peripheral collisions. Larger required statistics for peripheral events is a consequence of the lower average event multiplicity.
The robustness of the baryon stopping signal for a first-order phase transition against experimental cuts in the p T acceptance has been discussed in [18]. It is found that the baryon stopping signal is robust against hadronic rescattering.

Preliminary results on femtoscopy
Correlation femtoscopy allows one to measure the space-time characteristics of particle production in relativistic heavy-ion collisions due to the effects of quantum statistics and final state interactions. EPJ Web of Conferences  The femtoscopy was intensively studied at AGS and SPS accelerators and is being studied now at the BES/RHIC in the context of exploration of the QCD phase diagram. In this contribution we present preliminary results on femtoscopy observables for central Au-Au collisions at √ s NN = 7.7 GeV calculated by means of the event generator THESEUS. The femtoscopic analysis within THESEUS was performed very similarly to that described in Ref. [37]. A two-pion correlation function is fitted by the conventional Gaussian form where q = p 1 − p 2 with p 1 and p 2 being three-momenta of two considered particles, N is the normalization factor and λ is the correlation strength parameter, which can differ from unity due to the contribution of long-lived emitters and a non-Gaussian shape of the correlation function; R out , R side and R long are the Gaussian femtoscopy radii in the in the longitudinally co-moving system (LCMS), where the longitudinal pair momentum vanishes.
Preliminary results on R out , R side and R long as functions of the transverse mass of the particles, m T = (p 1 + p 2 ) 2 /4 + m π , for central Au-Au collisions at √ s NN = 7.7 GeV simulated by means of THESEUS without the UrQMD afterburner are presented in Fig. 8. The corresponding experimental data [38] are also displayed in Fig. 8. As seen, the femtoscopy radii are strongly underestimated without the afterburner stage. Figure 9 demonstrates these radii after the afterburner stage. Now these radii are already in a reasonable agreement with the STAR data. The R out /R side ratio is even in a very good agreement with the data. An unexpected result is that the femtoscopy radii turn out to be very close within the first-order-transition and crossover scenarios. The analysis of the femtoscopy observables is still in progress.

Conclusions
We have assembled the new event generator THESEUS that is based on a three-fluid hydrodynamics description of the early and dense stages of the collision, followed by a particlization as input to the UrQMD "afterburner" accounting for hadronic final state interactions. The new simulation program has the unique feature that it can describe a hadron-to-quark matter transition which proceeds in the baryon stopping regime that is not accessible to previous simulation programs that are designed for higher energies.
We presented first results from THESEUS for the FAIR-NICA-SPS-BES/RHIC energy scan addressing the directed flow of protons and pions as well as the rapidity distribution of protons, pions and kaons for two model EoS, one with a first order phase transition, the other with a crossover type softening at high densities. Preliminary results on the femtoscopy are also discussed. Another important application of THESEUS, not discussed in the present contribution, is prediction of the directed flow of deuterons in semicentral Au+Au collisions in the NICA energy range [39]. In Ref. [39] it is argued that light clusters are unique rare probes of in-medium characteristics such as phase space occupation and early flow.
We have found that the hadronic cascade which is switched on after the particlization has little effect on the proton flow observables and on the predicted baryon stopping signal for a first-order phase transition in heavy-ion collisions at NICA/FAIR energies. However, for pions in non-central collisions at lower energies the hadronic cascade leads to a qualitative change of the emission pattern (from flow to antiflow). The femtoscopy observables are strongly affected by the afterburner stage. The preliminary analysis of the femtoscopy radii manifested their reasonable agreement with the STAR data.