Particle production in a hybrid approach for a beam energy scan of Au+Au/Pb+Pb collisions between $\sqrt{s_\mathrm{NN}}$ = 4.3 GeV and $\sqrt{s_\mathrm{NN}}$ = 200.0 GeV

Heavy-ion collisions at varying collision energies provide access to different regions of the QCD phase diagram. In particular collisions at intermediate energies are promising candidates to experimentally identify the postulated first order phase transition and critical end point. While heavy-ion collisions at low and high collision energies are theoretically well described by transport approaches and hydrodynamics+transport hybrid approaches, respectively, intermediate energy collisions remain a challenge. In this work, a modular hybrid approach, the SMASH-vHLLE-hybrid coupling 3+1D viscous hydrodynamics (vHLLE) to hadronic transport (SMASH), is introduced. It is validated and subsequently applied in Au+Au/Pb+Pb collisions between $\sqrt{s_\mathrm{NN}}$ = 4.3 GeV and $\sqrt{s_\mathrm{NN}}$ = 200.0 GeV to study the rapidity and transverse mass distributions of identified particles as well as excitation functions for $\mathrm{dN}/\mathrm{d}y|_{y = 0}$ and $\langle p_\mathrm{T} \rangle$. A good agreement with experimental measurements is obtained, including the baryon stopping dynamics. The transition from a Gaussian rapidity spectrum of protons at lower energies to the double-hump structure at high energies is reproduced. The centrality and energy dependence of charged particle $v_2$ is also described reasonably well. This work serves as a basis for further studies, e.g. systematic investigations of different equations of state or transport coefficients.


Introduction
Collisions of atomic nuclei at relativistic velocities provide a unique opportunity to study and understand the fundamental properties of matter [1]. Such heavy-ion collisions are realized at different facilities, including for example the Relativistic Heavy-Ion Collider (RHIC) at BNL, the Large Hadron Collider (LHC) at Cern and SIS-18 at GSI. These different facilities provide, by virtue of distinct experiment geometries and accelerator dimensions (collision energies), access to different regions of the QCD phase diagram, ranging from high temperatures and vanishing baryon chemical potentials to moderate temperatures and high baryon chemical potentials. In recent years, heavy-ion collisions at intermediate energies have received particular attention in the quest for signatures of the postulated first order phase transition and critical end point [2,3]. In particular the BESII program at BNL [4], the NA61/SHINE experiment at Cern [5] as well as future FAIR [6], NICA [7] and J-PARC-HI [8] facilities have already or will soon provide experimental measurements in this designated region of interest [9]. There is a large variety of models to theoretically describe the underlying dynamics of heavy-ion collisions. At low collision energies, where the system is dominated by hadronic interactions, the dynamics are successfully captured by transport approaches relying on hadronic degrees of freedom [10][11][12][13][14]. At high collision energies, where a quark-gluon plasma state of matter becomes accessible, so-called hybrid models constitute the standard approach [15,[17][18][19]72]. Relativistic (viscous) hydrodynamic calculations are coupled to a hadronic afterburner; a setup successfully reproducing experimental observables [20]. For heavy-ion collisions at intermediate energies, where baryon densities are finite, there is however no such standard description yet. Previous works relying on different initial conditions, hydrodynamic evolutions and hadronic transport approaches include for example [21][22][23][24][69][70][71][72]. These have however either not been applied below √ s NN = 7.7 GeV, rely on relativistic hydrodynamics in 2+1D, or have issues properly reproducing baryon dynamics at low collision energies. The proper description of the baryon stopping dynamics is crucial for studies concerning fluctuations of conserved charges and to properly model the evolution at intermediate beam energies. Recent progress has been made in [25][26][27], where a framework towards more dynamical initialization of the hydrodynamic phase has been introduced. This is of importance for lower collision energies due to the longer-lived initial state.
In this work, the SMASH-vHLLE-hybrid is introduced. It is conceptually similar to the hybrid models presented in [21,22] but comprises of an improved description of the baryon stopping dynamics towards lower collision energies. The SMASH-vHLLE-hybrid relies on SMASH (Simulating Many Accelerated Strongly-Interacting Hadrons) [10,73] for the initial and the hadronic rescattering stage and on vHLLE [28,77] for the hydrodynamic evolution. It can be applied to model heavy-ion collisions ranging from √ s NN = 4.3 GeV to √ s NN = 5.02 TeV. It was already successfully employed to study the annihilation and regeneration of (anti-)protons and (anti-)baryons at a variety of collision energies and centralities [29]. Here, the details of the SMASH-vHLLE-hybrid are explained, the approach is validated and results for particle spectra and excitation functions confronted with experimental data, where a good agreement across a wide range of collision energies is observed. In particular, the proton dN/dy spectra as well as the proton p T excitation function agree well with experimental measurements.
This work is structured as follows: In Sec. 2 the SMASH-vHLLE-hybrid is explained with particular emphasis on the individual modules. The equation of state of the SMASH hadron resonance gas, which is of fundamental importance for the particlization process, is presented in Sec. 2.4. The approach is validated in Sec. 3 by demonstrating consistency at the interfaces, global conservation of quantum numbers and good agreement with other hybrid approaches consisting of SMASH+CLVisc [23] and vHLLE+UrQMD [21]. In Sec. 4, results obtained with the SMASH-vHLLE-hybrid are presented. First, dN/dy and dN/dm T spectra for π − , p and K − in Au+Au/Pb+Pb collisions ranging from √ s NN = 4.3 GeV to √ s NN = 200.0 GeV as well as excitation functions for dN/dy| y=0 and p T are confronted with experimental data. A good agreement is obtained for pions, kaons and also the usually challenging protons. Second, the centrality and collision energy dependence of the charged particle integrated v 2 is well reproduced.
To conclude, a brief summary and outlook are provided in Sec. 5.

Model Description
The SMASH-vHLLE-hybrid [76] presented in this work is a novel hybrid approach for the theoretical description of heavy-ion collisions between √ s NN = 4.3 GeV and √ s NN = 5.02 TeV. It is publicly available on Github 1 .
In general, hybrid approaches are a very successful combination of viscous fluid dynamics for the hot and dense region, while microscopic non-equilibrium transport approaches are applied for the initial and final stages where the system is far away from local equilibrium. In this work, we rely on the same transport approach (SMASH) for initial and final stages to allow for a smooth transition from the low energy regime, where hadronic transport approaches describe the full evolution to the higher energy region, which is where the phase transition to the quark-gluon plasma is expected. The SMASH-vHLLE-hybrid consists of multiple submodules: The initial state is provided by the hadronic transport approach SMASH [10,73], while the evolution of the hot and dense fireball is described by vHLLE, a 3+1D viscous hydrodynamics model [28,77]. The hydroynamical evolution is performed until the medium gets too dilute, estimated with the energy density e falling below a critical value of e crit = 0.5 GeV/fm 3 . Particlization is realized with the SMASH-hadron-sampler [75] and the subsequent hadronic afterburner evolution is again performed with the SMASH transport approach. These submodules are explained in more detail in the following.

SMASH
SMASH is a hadronic transport approach designed for the description of low-energy heavy-ion collisions and hadronic matter in and out of equilibrium [10,73]. It provides an effective solution of the relativistic Boltzmann equation by modeling the collision integral C(f ) through formations, decays and elastic scatterings of hadronic resonances. In Eq.
(1), f denotes the one-particle distribution function, p µ the particle's 4-momentum, m its mass and F µ an external force. The SMASH degrees of freedom consist of all hadrons listed by the PDG up to a mass of m ≈ 2.35 GeV [30]. Resonances are represented by vacuum Breit-Wigner spectral functions, but have mass-dependent widths following the Manley-Saleski ansatz [31]. The resonance properties are tuned to reproduce elementary cross sections. For high-energy hadronic interactions, a string model is employed where string fragmentation and the high energy hard scatterings are carried out within Pythia 8 [32,33]. In addition, SMASH has the perturbative production of dileptons [34] and photons [35] integrated. In this work, a geometric collision criterion is applied that is formulated in a covariant form [36], a stochastic collision criterion is however also implemented [37]. Inter alia, SMASH was already successfully applied to study particle production in a variety of collision setups [10,[37][38][39]. For a broader introduction into SMASH and its different features, the interested reader is referred to [10,38]. In this work, the nuclear mean fields are turned off and only the cascade mode of SMASH is applied, which is appropriate for the beam energies under consideration.

Initial conditions
The initial conditions for the hydrodynamical evolution are extracted from SMASH on a hypersurface of constant proper time τ 0 . By default, this proper time corresponds to the passing time of the two nuclei [21] and is determined from nuclear overlap via where R p and R t are the radii of projectile and target, respectively. √ s NN is the collision energy and m N is the nucleon mass. The passing time constitutes a reasonable assumption for the time where local equilibrium is reached [41] and for high beam energies a lower limit of τ 0 = 0.5 fm/c is implemented. To determine the isoτ hypersurface, the SMASH evolution, which is realized in Cartesian coordinates, is initialized and performed including all scatterings and interactions until the particles reach the iso-τ hypersurface. Upon crossing this hypersurface, they are removed from the evolution and stored such that they can be later utilized to initialize the subsequent hydrodynamical evolution. One might argue that the removal of particles introduces holes in the density distribution of the medium and could therefore modify the underlying dynamics. We have however verified that the final state observables remain unaffected if the particles are not removed from the SMASH evolution but can undergo further interactions.

Hadronic Rescattering
In the last stage of the SMASH-vHLLE-hybrid, SMASH is again applied to model the non-equilibrium hadronic afterburner evolution. For this, the hadrons obtained from particlization are read from an external file, backpropagated to the earliest time and successively appear once fully formed. They free stream before reaching their respective formation times. The afterburner stage is characterized by hadronic rescatterings as well as resonance decays. Both kinds of interaction take place until the medium is too dilute.

vHLLE
The evolution of the hot and dense fireball is modeled by means of the 3+1D viscous hydrodynamics code vHLLE [28]. It provides a solution of the hydrodynamic equations with the energy-momentum tensor decomposed as follows: where is the local rest frame energy density, p and Π are the equilibrium and bulk pressure, and π µν is the shear stress tensor. Eqs. (3) and (4) encapsulate the conservation of energy and momentum via the energy momentum tensor T µν as well as the conservation of net-baryon, net-charge and net-strangeness numbers N B , N Q , N S , respectively. The code solves relativistic viscous hydrodynamic equations in the second-order Israel-Stewart framework in the 14-momentum approximation, including shear-bulk coupling terms [66,67]. That yields the following relaxation type equations for the viscous corrections: Dπ µν = 2ησ µν − π µν τ π − δ ππ τ π π µν θ + φ 7 τ π π µ α π ν α − τ ππ τ π π µ α σ ν α + λ πΠ τ π Πσ µν ., where θ is expansion scalar, σ µν is the stress tensor, τ π and τ Π are the relaxation times for the shear and bulk corrections and λ Ππ , λ πΠ , δ ππ , τ ππ and φ 7 are higher-order couplings, whose ansatzes are taken from [66]. In this work, we only apply fixed values of effective shear viscosity over entropy density for the hydrodynamic evolution. vHLLE is a Godunov-type algorithm in conservative form. In its original form, in Cartesian coordinates, the algorithm preserves the total energy, momentum and conserved charges in the system by construction. However, hydrodynamic equations in Milne coordinates contain non-vanishing source terms, which are integrated numerically with finite accuracy. Therefore, in Milne coordinates the total energy and momentum of the system varies a little, typically increasing from the beginning of the hydrodynamic evolution towards its end by few percents [44].
The hydrodynamical evolution is initialized with the energy-momentum tensor calculated from the iso-τ particle list created with SMASH as detailed in Sec. 2.1.1. To prevent shock waves, these particles are smeared relying on a Gaussian smearing kernel as described in detail in [21]. The magnitude of the smearing is controlled by R η and R ⊥ , the longitudinal and transversal smearing parameters, respectively. The specific values used in this work are provided in Sec. 2.5. Once the particles are transformed into fluid elements, the medium is evolved relying on a chiral model equation of state [42] which is matched to a hadron resonance gas at low energy densities. The equation of state [42] is constructed for isospin symmetric matter, i.e. zero electric charge chemical potential µ Q = 0. Nevertheless, the electric charge density is included in the simulation. As such, the density is computed in the initial state and is evolved in the hydrodynamic stage. At particlization, a hadron gas equation of state is used, which includes both charges and corresponding chemical potentials and as such is not limited to isospin symmetric matter. In contrast to the baryon and electric charge densities, the strangeness density is set to zero in the hydrodynamic stage, as explained in Sec. 2.4. The viscous evolution is performed until the medium drops below the predefined critical energy density of e crit = 0.5 GeV/fm 3 and the corresponding freezeout hypersurface is created. The latter is constructed dynamically during the hydrodynamical evolution with the CORNELIUS subroutine [43]. For the sake of quantum number conservation in the subsequent particlization process, it is important that the thermodynamical properties of the hypersurface elements match those of the SMASH hadron resonance gas [44]. The SMASH equation of state is thus used to re-calculate the temperature and chemical potentials at each element. It is introduced in Sec. 2.4. The hydrodynamical evolution terminates once all cells have dropped below the critical energy density e crit = 0.5 GeV/fm 3 . An energy density criterion results in a reasonable transition region in the QCD phase diagram. For a more detailed introduction into vHLLE as well as its initialization from an external particle list, the reader is referred to [21, 28].

SMASH-hadron-sampler
To perform the non-equilium hadronic afterburner evolution it is necessary to transform the fluid elements on the freezeout hypersurface into particles, a process usually referred to as particlization. In the presented model, the SMASH-hadron-sampler [75] is applied to achieve particlization according to the properties of the SMASH hadron resonance gas. Grand-canonical ensemble is employed, which allows to particlize each surface element independently. Technically this is realized in two steps: first, the sum of thermal multiplicities N tot i of all hadron species 2 passing through each surface element i is computed. Next, an actual number of hadrons to generate at each surface element N gen i is sampled according to Poisson distribution with a mean N tot i . If, for a given surface element N gen i > 0 then the generation proceeds to the second step: the type of a hadron is sampled based on the relative abundances from the thermal average, and the Cooper-Frye formula [45] is applied to sample the momentum of each of the N gen i hadrons: In Eq. (8), N denotes the species' abundance, d its degeneracy factor, p its 3-momentum, p µ its 4momentum, and E p its energy. Furthermore, dΣ µ is the normal vector to the hypersurface patch, f 0 is the vacuum one-particle distribution function and δf shear and δf bulk are the viscous and bulk corrections to f 0 , respectively. For the presented study, only shear viscous corrections are considered though. The integration is performed over the full freezeout surface Σ to account for all contributions from the individual hypersurface 2 Note though, that we exclude the leptons e ± , µ ± , τ ± , the photon γ and the σ meson from the SMASH particle list for the sampling process. These are not among the usual hadronic degrees of freedom of SMASH , but form part of the particle list to allow for additional features. They are also excluded for the determination of the SMASH equation of state in Sec. 2.4.

patches.
For further information about the specific sampling procedure, the reader is referred to Section 2 C of [21]. Note that, in the SMASH-hadron-sampler, because of the grand-canonical nature of the sampling procedure, quantum numbers are not conserved eventby-event, but only on average. Since we are only interested in averaged single-particle observables, this approximation is sufficient at this point. The resulting particle list is subsequently used to initialize the nonequilibrium hadronic afterburner evolution as described in Section 2.1.2 to model the late and dilute hadronic rescattering stage.

SMASH hadron resonance gas equation of state
As described in Sec. 2.2, the equation of state of the hadron gas corresponding to the degrees of freedom in the afterburner evolution is an essential ingredient for the creation of the freezeout hypersurface of the hydrodynamic stage. Otherwise, quantum number conservation in the particlization process is impossible. The equation of state of the hadron gas made up of the SMASH degrees of freedom is thus briefly described in the following. The equation of state (EoS) generally provides a mapping of the thermodynamic quantities energy density e, net baryon density n B , net charge density n Q and net strangeness density n S to the temperature T , the pressure p, the baryon chemical potential µ B , the charge chemical potential µ Q , and the strange chemical potential µ S . For heavy-ion collisions however, the net strangeness density can be approximated as n S = 0 fm −3 . We thus neglect the explicit dependence on the net strangeness density. Whereas the hydrodynamic code is capable of evolving the strangeness current, in the present study n S = 0 is set in the initial state for the hydrodynamical stage. As such, local net strangeness is zero in the hydrodynamical as well as the particlization stage, which is consistent with the hadrnoic equation of state employed. Hence, the SMASH equation of state presented herein provides the mapping: This mapping follows from the properties of the underlying gas of hadrons and can be determined by solving the set of coupled equations assuming an ideal Boltzmann gas within the grandcanonical ensemble [46]. The solutions of Eqs. (10) can numerically be obtained with a root solver algorithm. Unfortunately however, this solver is highly-sensitive to the choice of the initial approximation and might thus fail to converge. This is particularly problematic at low energy densities, where the gas is loosely populated, and close to the kinematic thresholds 3 . The SMASH equation of state presented herein is thus accurate for high energy densities as well as low baryon and charge densities, else it constitutes an approximation. The latter is obtained from a combination of (i) variation of the grid size within the solver algorithm, (ii) variation of the initial approximations, and (iii) interpolations and extrapolations of the previously obtained solver results, where suitable. The resulting equation of state, ranging from e = 0.01 GeV/fm 3 up to e = 1.0 GeV/fm 3 , and accounting for all hadronic degrees of the SMASH hadron resonance gas, is available in tabularized format under [78]. For validation purposes, the SMASH equation of state at vanishing net baryon, net electric charge and net strangeness densities is compared to lattice QCD results obtained in 2+1-flavour QCD [47] in Fig. 1

Model Validation
Before applying the SMASH-vHLLE-hybrid to study different observables in heavy-ion collisions, we first want to systematically validate the approach. This is achieved with four different checks: First, the interfaces at the initial and final state of vHLLE are investigated with particular emphasis on the smeared initial state and the properties of the cells on the freezeout hypersurface, respectively. Second, the approximate global conservation of quantum numbers throughout the different stages of the evolution is demonstrated. Third, an apples-to-apples comparison to a hybrid model consisting of SMASH and CLVisc is performed, and fourth, results of the SMASH-vHLLE-hybrid are compared to results from a UrQMD+vHLLE hybrid model, as used in [21]. These validations and cross checks are elaborated on in more detail in the following.

Interface with vHLLE
It is well known that the interfaces between the non-equilibrium initial and final evolution and hydrodynamics are crucial in any hybrid approach. During the initial transition from SMASH to vHLLE the assumption of being close to local equilibrium leads to discontinuities. On the Cooper-Frye hypersurface, there is again a somewhat sharp transition between the strongly coupled fluid evolution and a hadronic transport approach with vacuum properties. Those fundamental problems are beyond the current work and are existent in all other hybrid approaches applied to heavy-ion collisions as well. Here, we are going to concentrate on ensuring minimal consistency in terms of matching quantum numbers at the initial interface.
As described in Sec. 2.2, it is necessary to smear the initial SMASH particle list upon initialization of vHLLE. The effect of this smearing is demonstrated in Fig. 2, containing the dE/dη (upper), dB/dη (center), dQ/dη (lower) distribution as a function of space-time rapidity η. The solid lines correspond to one single underlying SMASH event in a Pb+Pb collision at √ s NN = 8.8 GeV and the dashed line to its smeared counterpart, indicating the initial condition for the hydrodynamical evolution. It can be observed that the spiky structure of the underlying SMASH event is smoothed successfully, without significant loss of quantum numbers, thus validating the handling of the SMASH → vHLLE interface. Fig. 2 gives also an idea on how large the event-by-event fluctuations in our initial state are.  GeV. The SMASH state (solid lines) is compared to the smeared state used to initialize the hydrodynamical evolution in vHLLE (dashed line). Ideal hydrodynamics is applied and spectators are excluded.

Global Conservation Laws
The conservation of quantum numbers throughout all stages of the hybrid approach is essential in order to properly describe the evolution of the collision system. In this section we thus study the conserved quantities (energy E, net baryon number B and net electric charge Q) as they evolve through the different modules and stages of the SMASH-vHLLE-hybrid. Their evolution is presented in Fig. 3, where the left panel refers to the total energy E, the center panel to the net baryon number B, and the right panel to the net electric charge Q. These plots are to be interpreted as follows: The evolution of the conserved quantity from the initial to the final stage of the SMASH-vHLLE-hybrid is presented from left to right. The value displayed for E, B, or Q is normalized to its respective initial value to identify deviations more easily. The segments of curves connect the normalized values at the beginning and the end of each stage in the modelling. As quantum number conservation is enforced in SMASH there are no violations up to the SMASH → vHLLE interface. Deviations in the vHLLE stage represent a cumulative change to the quantum numbers from two sources. The first source is the finite accuracy of the hydrodynamic algorithm, due to the used Milne coordinate frame. In the Milne frame, geometrical source terms in the fluid-dynamical equations need to be time-integrated numerically, and the integration has a finite accuracy. The second source stems from the SMASH → vHLLE interface, however for visual purposes we combine it with the first source and display the result as an overall energy deviation during the vHLLE stage. At the initialization of the hydrodynmical evolution there is a small fraction of cells, typically located at the very borders of the fireball, whose energy density is already below e ≤ 0.01 GeV/fm 3 and thus fall outside the applicability range of the SMASH equation of state (see Sect. 2.4). Lacking possibilities to assign thermodynamic quantities (T, p, µ B , µ Q , µ S ) to these cells, we have decided to neglect their contribution, for the sake of not introducing additional uncertainties. Naturally, this treatment results in loss of E, B, and Q, which becomes more severe the lower the collision energy, as relatively more cells are characterized by small energy densities. We have however verified, that at most 1.57 % of the total energy, 0.74 % of the total baryon number, and 0.29 % of the total electric charge are lost in Au+Au collisions at √ s NN = 4.3 GeV. Combining the afore-explained effects, the resulting violations in the SMASH-vHLLE-hybrid can be quantified to be of the order ≤ 5 % up to completion of the vHLLE stage. In the subsequent sampling process, quantum numbers are approximately conserved. Only small violations are found, which are most pronounced at low collision energies. These violations are related to the fact that, as described in Sec. 2.4, the hadron gas equation of state in this work is not perfectly accurate across the entire (e, n B , n Q )-range, but constitutes an approximation in the problematic regions. Finally, in the SMASH afterburner stage, quantum numbers are exactly conserved again, owing to its enforcement within the SMASH approach. Further discussions can be found in [44].
Regarding the entire evolution of the system within the SMASH-vHLLE-hybrid, we find violations of quantum number conservation by at most 7 % for collisions between √ s NN = 4.3 GeV and √ s NN = 200.0 GeV, which we consider an overall good validation.

Comparison to CLVisc
The evolution up to the Cooper-Frye hypersurface is now validated by comparing its outcome to a hybrid model in which the hydrodynamical evolution is performed by CLVisc Results obtained with vHLLE are denoted with solid lines, those obtained with CLVisc with filled circles. We find a perfect agreement for the dN/dy spectra as well as for the anisotropic flow coefficients between the vHLLE and CLVisc hydrodynamical evolution, relying on an identical initial state. This is considered another validation of the presented model, in particular regarding the hydrodynamical evolution. While both fluid dynamics codes have been very well tested, it is very nice to see that they lead to exactly the same results starting from the same SMASH initial condition.

Comparison to UrQMD+vHLLE hybrid approach
To finally put the results obtained with the SMASH-vHLLE-hybrid into context with other state-ofthe-art approaches, particle spectra are compared to those obtained within a hybrid model consisting of UrQMD and vHLLE, which was successfully applied to 4 The corona region denotes the collection of cells whose energy density is below the critical energy density already at the initialization of the hydrodynamical evolution. These cells are directly written to the freezeout hypersurface and are not further propagated in the hydrodynamical evolution.  The conceptual ideas of the SMASH-vHLLE-hybrid and the UrQMD+vHLLE hybrid approach are identical. The only difference is the transport approach applied to provide the initial state and to perform the afterburner evolution. Although UrQMD and SMASH are conceptually similar, there are small differences regarding  the degrees of freedom, cross sections, or in the application of Pythia to perform the string excitations. In particular the latter is important to properly model the baryon dynamics at intermediate and high collision energies [38]. Fig. 5 shows the rapidity distributions dN/dy (upper) and m T spectra (lower) of π − and K − in Pb+Pb collisions at √ s NN = 8.8 GeV and 17.3 GeV as obtained within the SMASH-vHLLE-hybrid (solid lines) and within the UrQMD+vHLLE hybrid (dashed lines, taken from [21]). For consistency, identical smearing parameters and viscosities are applied to model the hydrodynamical evolution in both cases. Note that these differ from the ones listed in Table 1, the applied values are indicated directly in the figures. It can be observed that, at both collision energies, the dN/dy distribution of pions from the SMASH-vHLLE-hybrid is wider than the one from the UrQMD+vHLLE hybrid, but peaks at lower midrapidity yields. The shape of the kaon dN/dy distribution is similar, but the UrQMD+vHLLE hybrid produces smaller kaon yields than the SMASH-vHLLE-hybrid. The m T spectra of both approaches are in good agreement, only at low transverse masses the UrQMD+vHLLE hybrid produces more pions but less kaons than the SMASH-vHLLE-hybrid, which is not surprising given the differences in the dN/dy spectra. It is reassuring to find that both hybrid approaches provide similar results in terms of dN/dy and m T spectra. This puts the novel hybrid approach introduced in this work into context with other approaches existent in the fied and gives some hints towards differences between the two transport approaches SMASH and UrQMD.

Results
In this section, bulk hadronic observables are presented as obtained from the SMASH-vHLLE-hybrid in heavyion collisions at different collision energies. This includes dN/dy and m T spectra of pions, protons and kaons as well as their excitation functions for the midrapidity yield, p T , v 2 , and v 3 . The results from the hybrid approach are compared to a pure transport calculation to indicate the differences due to the intermediate hydrodynamic evolution. Before the final results for observables are discussed, we show results on the location of our switching transition between vHLLE and on a hypersurface of constant energy density in the phase diagram.
Here, A is a placeholder for the quantity of interest, T or µ B , N is the total number of hypersurface elements, and e n is the energy density of the respective hypersurface element. In Fig. 6, the results for central collisions are presented, for a centrality-dependent version of this figure, the reader is referred to Appendix A.
The solid line in Fig. 6 corresponds to a parametrization of the chemical freezeout line, deduced from experimentally measured hadron abundances within the statistical hadronization model [48,49]. Regarding the locations of the Cooper-Frye transition surface within the SMASH-vHLLE-hybrid it can be observed that, as expected, higher-energy collisions are characterized by higher temperatures and lower baryon-chemical potentials, while lower collision energies by relatively lower temperatures and higher baryon-chemical potentials. This is generally in line with the shape of the parametrization. Let us also point out here that the spread of the system is significant even in single events, so a large part of the phase diagram can be covered with a limited range of beam energies (see also [50]). As expected the transition from hydrodynamics to hadronic transport happens slightly above the chemical freezeout line, since the system should be mainly hadronic at that stage. The actual chemical freezeout happens naturally within the rescattering dynamics, when the inelastic collisions cease depending on the hadronic species.    (solid lines) are compared to those obtained when running only SMASH (dashed lines). The E866 data is taken from [56,57], the STAR data from [58], and the NA49 data from [52][53][54][55].

Hadron production: dN/dy and dN/dm
energies. While at low energies the single-peak structure in proton dN/dy spectra persist, baryon transparency has a greater impact for rising collision energies [59]. The SMASH-vHLLE-hybrid is successful at describing the single-peak structure observed at √ s NN = 4.3 GeV as well as the double-peak structure at √ s NN = 7.7 GeV and √ s NN = 17.3 GeV. The new hybrid approach is thus able to reproduce the magnitude of baryon stopping at intermediate collision energies. The agreement with experimental data for proton dN/dy spectra is improved significantly with the application of the SMASH-vHLLE-hybrid instead of a pure transport evolution, yet the final agreement is not perfect. Let us note in that respect that finite baryon diffusion will also affect the proton rapidity distribution, which is not yet included in the SMASH-vHLLE-hybrid [26,27]. Nevertheless, the qualitative (and almost quantitative) agreement with the proton rapidity distribution over a large range of beam energies is considered an important pre-requisite for further studies of phase transition signals.
In the case of the dN/dm T spectra in Fig. 7, a significant hardening of the π − , p, and K − spectra is observed at all collision energies if the SMASH-vHLLE-hybrid is employed instead of a pure transport evolution utilizing SMASH. These results are again in line with observations made in other hybrid approaches [21,22]. The observed hardening of the spectra in the hybrid setup noticeably improve the agreement with experimental measurements as compared to a pure transport evolution at all collision energies probed in Fig. 7, in particular regarding the slopes of the m T spectra of all particle species. Summarizing Fig. 7 it can be stated that the SMASH-vHLLE-hybrid approach is reasonably describing dN/dy and dN/dm T spectra of π − , p, and K − in Au+Au/Pb+Pb collisions at √ s NN = 4.3, 7.7 and 17.3 GeV. A decent agreement with experimental measurements at these energies is obtained. Also, the baryon stopping power is well reproduced at the above listed energies, thus demonstrating that longitudinal baryon dynamics are captured properly.

Excitation functions
To investigate the particle production over a broader range of collision energies, the midrapidity yield and mean p T are calculated for π − , p and K − and compared to the pure transport calculation. Reproducing the yield and the mean transverse momentum is almost equivalent to reproducing the full transverse mass spectra. Au+Au/Pb+Pb collisions are calculated between √ s NN = 4.3 GeV and √ s NN = 200.0 GeV, again relying on the viscosities and smearing parameters listed in Table 1. As above, central collisions are modeled where the impact parameter range for 0-5% central collisions is approximated by b ∈ [0, 3.3] fm and b ∈ [0, 3.1] in Au+Au and Pb+Pb collisions, respectively. For centrality-dependent excitation functions, the interested reader is referred to Appendix A.
In Fig. 8   the SMASH-vHLLE-hybrid shows a good agreement with experimental measurements. There is a nearly-perfect agreement with the measured π − and p yields, the K − yields are systematically overestimated though. The latter could be improved by more dynamical initial conditions, as recently realized in [25]. It can further be observed that the excitation function from the SMASH-vHLLE-hybrid shows a smooth behaviour for rising collision energies, while a kink is observed between √ s NN = 6 GeV and √ s NN = 8 GeV in the pure SMASH curve. This kink is unphysical and stems from the non-trivial transition from resonance dynamics to string dynamics in SMASH. The SMASH-vHLLE-hybrid is thus better suited to smoothly and consistently describe the dynamics at intermediate collision energies. Let us note here, that it is well known that the application of a pure transport model is not sufficient towards higher collision energies, as too little radial flow is produced. At the same time, pure hydrodynamics is not sufficient either for lack of non-equilibrium dynamics in the initial and final stages. In hydrodynamics+transport models, transport theory and hydrodynamics are employed in their respective regions of applicability to provide as good a description of heavy-ion collisions as possible. Ideally, unphysical properties arise only if these models are employed far from their region applicability. For the p T | y=0 excitation function displayed in the lower panel of Fig. 8 the agreement with experimental measurements is also improved once a hybrid model is applied instead of a pure SMASH evolution. Most importantly, the p T of protons rises with rising collision energies, which is in line with experimental observations, while in the pure SMASH case it shows a decreasing trend for collisions above √ s NN ≈ 30 GeV. The application of a hybrid model is thus essential to properly capture the transversal baryon dynamics, especially towards higher collision energies.

Collective flow
Third, the integrated elliptic flow v 2 and integrated triangular flow v 3 of charged particles are determined for Au+Au/Pb+Pb collisions between √ s NN = 4.3 GeV and √ s NN = 200.0 GeV as modeled with the SMASH-vHLLE-hybrid in 0-5%, 5-10%, 10-20%, 20-30%, 30-40%, and 40-50% most central collisions 5 . The event plane method is employed and the resulting v 2 and v 3 is compared to experimental data from the STAR collaboration [64]. Note though, that the STAR results are obtained from two-particle correlations, hence the comparison is not perfectly equivalent. It has been shown however that both methods yield similar, yet not identical, results [65]. The integrated v 2 (upper panel) and integrated v 3 (lower panel) are presented as a function of collision energy √ s NN in Fig. 9. More central collisions are marked by darker colours and more peripheral collision by lighter colours.
Regarding v 2 , it can be observed that for central collisions the v 2 extracted from the SMASH-vHLLE-hybrid is in decent agreement with measurements from STAR, while for more peripheral collisions the agreement is good at high collision energies, but the measured v 2 is significantly underestimated towards lower collision energies. It can further be seen that the v 2 determined from the SMASH-vHLLE-hybrid in 20-30%, 30 Fig. 9 Integrated v 2 (upper) and integrated v 3 (lower) excitation functions of charged particles at midrapidity as a function of collision energy, obtained with the SMASH-vHLLE-hybrid in Au+Au/Pb+Pb collisions. More central collisions are marked by darker colours, more peripheral collisions by lighter colours. The STAR data is taken from [64].
much more than in more central collisions when moving from higher to lower collision energies. This can be explained with a too short lifetime of the hot and dense fireball in the hydrodynamic stage. In Fig. 10 the mean lifetime of the hydrodynamic stage, τ hydro , is presented as a function of collision energy for different centralities. Again, more central collisions are marked by darker colours and more peripheral collision by lighter colours. The hydrodynamical lifetime is determined via where τ 0 is the proper time at initialization of hydrodynamics and τ end the proper time at which the last cell falls below the critical energy density. events denotes the average over all events. It can be seen that the lifetime of the fireball decreases continuously for increasing centralities. This is expected from the geometrically smaller overlap region resulting in smaller volumes and thus shorter-lived fireballs. One can further observe that the hydrodynamical lifetime is nearly independent of the collision energy above √ s NN ≈ 50 GeV.
Below √ s NN ≈ 50 GeV however, it is reduced substantially with decreasing collision energies. This implies, that especially in peripheral collisions at low and intermediate collision energies, the hydrodynamical evolution lasts only a few fm/c, which is seemingly too short for collective dynamics to develop. The resulting v 2 is thus significantly underestimated. This might be alleviated by implementing more dynamical initial conditions that are better suited to capture the underlying dynamics at low collision energies, as recently done in [25].
Regarding v 3 on the other hand (c.f. lower panel of Fig. 9), the experimentally measured v 3 is underestimated across the entire energy range and for all centralities by the SMASH-vHLLE-hybrid. This might be related to the Gaussian smearing employed at initialization of hydrodynamics to provide smooth density profiles. Initial state fluctuations, which are responsible for the generation of v 3 , might thus be smeared too much and the resulting v 3 is hence underestimated. Additionally, a too short hydrodynamical evolution, as detailed above, could also be responsible for an underestimation of the integrated v 3 .
In general, collective flow observables are very sensitive to the presence of a fluid dynamic evolution and to the employed transport coefficients. It will be interesting to see, how these coefficients behave when more realistic temperature and net baryon density dependent transport coefficients are implemented in the hybrid approach.

Conclusions and Outlook
In this work a novel modular hybrid approach based on the hadronic transport approach SMASH coupled to the 3+1D with experimental data improved for π − , p and K − as compared to employing a pure SMASH evolution, especially in view of baryon dynamics. The integrated v 2 and v 3 of charged particles as a function of collision energy were further confronted with STAR data for different centralities. A good agreement was found for v 2 in central collisions as well as at high collision energies. Otherwise, v 2 as well as v 3 were underestimated, which can be explained with a too short lifetime of the hydrodynamical fireball as well as with the smearing kernel employed at the initialization of hydrodynamics.
The SMASH-vHLLE-hybrid was successfully applied to study a broad range of hadronic observables relying on initial conditions extracted on a hypersurface of constant proper time. This assumption is however questionable at low collision energies, where the dynamics tend to be comparably slow. The extension of the SMASH-vHLLE-hybrid by more dynamical initial conditions, similar to efforts made in [25], thus constitute the next step, to improve the applicability at FAIR/NICA collision energies. In the future, the SMASH-vHLLE-hybrid can be applied with different equations of state to evolve the hot and dense fireballwith or without a critical end point or first order phase transition -as well as with varying transport coefficients to systematically study the impact on particle spectra, flow, or other observables.
In addition to the results presented in Sec. 4.3, the p T excitation functions of π − , p and K − were extracted from the SMASH-vHLLE-hybrid for different centralities between 0-5% and 40-50%. These results are presented in Fig. 11, where π − are displayed by dashed-dotted lines, p with solid lines and K − with dashed lines. Darker colors denote more central collisions, lighter colours more peripheral collisions. It can be observed that the p T of all three particle species generally depend little on the centrality of the collision, yet more peripheral collisions result in smaller mean transverse momenta of pions, protons and kaons. The centrality dependence increases with decreasing collision energies and is most pronounced in the case of protons. The properties of the freezeout hypersurface of the hydrodynamical stage, as presented for central collisions in Fig. 6, can further be analyzed as a function of collision centrality. This is displayed in Fig. 12, where central Au+Au/Pb+Pb collisions are again represented by darker coloured markers, more peripheral collisions by lighter coloured markers. The collision energies rise from √ s NN = 4.3 GeV in the lower right corner to √ s NN = 200.0 GeV in the upper left corner. It can be observed that the general dependence on the collision centrality is moderate. Especially for high-energy collisions, the freeze-out properties are nearly unaffected by centrality variation. Towards lower collision energies a stronger centrality dependence is observed, but it is small.
Appendix B: Parametrizing the SMASH hadron resonance gas equation of state As mentioned in Sec. 2.4 we have made different attempts to parametrize the SMASH equation of state for it to be applied more easily within the SMASH-vHLLE-hybrid and potentially also in other works. Unfortunately, these efforts remain unsuccessful up to now. We still want to briefly document the different attempts we made in the following. The purpose of the equation of state is to perform the mapping (e, n B , n Q ) → (T, p, µ B , µ Q , µ S ). Currently, the coupled equations (c.f. Eqs. 10) are solved once for different combinations of e, n B and n Q and stored in a table, where the resulting thermodynamic quantities need to be looked up for a given set of (e, n B , n Q ) and interpolated between the grid points, if necessary. This procedure is time-consuming and introduces uncertainties related to finite grid effects when interpolating between the predefined grid points. An alternative solution would be to directly solve the coupled thermodynamic equations directly within the hydrodynamical simulation for each cell individually, relying on its specific densities. This procedure is however even less effective and requires large computational resources. A parametrization from which one could directly calculate the thermodynamic quantities and that relies on a limited number of parameters only would thus be the most effective solution.
We have tried to use the following polynomial-inspired ansatzes in order to find a parametrization of the SMASH hadron resonance gas equation of state: 1. f (e, n B , n Q ) = Here, a ijk , b ijk , or c ijk denote the sets of parameters we want to determine. The major problem we faced was the rapidly growing number of parameters N Param for larger N Dim (for example in the second case N Dim = 3 results in N Param = 55). Large values for N Dim are however necessary as the equation of state is a non-trivial function and small values for N Dim resulted in unsatisfactory parametrizations. In addition, a good initial approximation of a ijk , b ijk , and c ijk is essential for the solver to succeed. Besides a uniform sampling in parameter space, we have also explored a Latin Hypercube Sampling (LHS) procedure [68], where the parameter space is covered in a more efficient way to provide better initial guesses. Both ways of determining initial guesses and later checking for a best fit did however not result in any satisfying parametrization of the underlying equation of state. On top of this we utilized Chebyshev nodes as grid points to minimize strong oscillations at the edges of the equation of state, which did not alleviate our issues either. In summary, wie explored different possibilities to find a suitable parametrization for the SMASH hadron resonance gas equation of state but did unfortunately not succeed. We thus fall back to our look-up