SMASH-vHLLE Hybrid Approach

• SMASH-vHLLE hybrid approach is a modular framework that integrates microscopic hadronic transport with 3+1D viscous hydrodynamics to simulate the full evolution of heavy-ion collisions.
• It dynamically switches between transport and hydrodynamic modules to capture non-equilibrium early/late stages and the dense, approximately equilibrated mid-stage, reproducing key observables.
• The framework has been validated across collision energies from a few GeV to TeV scales and extended via Bayesian calibration and studies in small systems.

The SMASH-vHLLE hybrid approach is a modular transport-plus-hydrodynamics framework for relativistic heavy-ion collisions in which SMASH provides the microscopic hadronic transport description of the early pre-equilibrium stage and the late dilute afterburner, while vHLLE evolves the intermediate approximately locally equilibrated stage in full $3+1$ dimensions with viscous hydrodynamics. A dedicated hadron sampler implements Cooper–Frye particlization on a freeze-out hypersurface, with thermodynamic matching to a SMASH-consistent hadron resonance gas equation of state. The framework was introduced for collisions from $\sqrt{s_{NN}}=4.3$ GeV up to $5.02$ TeV, then extended to Bayesian calibration, dynamic fluidization at low beam energies, and exploratory applications to small systems such as O–O and Ne–Ne at $\sqrt{s_{NN}}=5.36$ TeV (Schäfer et al., 2021, Schäfer et al., 2021, Götz et al., 13 Mar 2025, Constantin et al., 6 Sep 2025).

1. Concept, scope, and modeling rationale

The central motivation of the SMASH-vHLLE program is to describe collision systems for which neither a pure hadronic cascade nor a purely macroscopic fluid picture is adequate over the full space-time history. In the high-baryon-density region of the QCD phase diagram, early and late stages are strongly non-equilibrium and naturally described by transport, whereas the dense intermediate stage may be close enough to local equilibrium to justify viscous hydrodynamics. The framework is therefore explicitly a dynamical hybrid: transport $\rightarrow$ hydrodynamics $\rightarrow$ transport (Schäfer et al., 2021).

This design was first validated for Au+Au and Pb+Pb collisions in the beam-energy-scan regime, where it reproduced rapidity and transverse-mass distributions of identified hadrons, excitation functions for $\mathrm dN/\mathrm dy|_{y=0}$ and $\langle p_T\rangle$, and a reasonable description of charged-particle $v_2$. In particular, it reproduced the transition from a Gaussian proton rapidity spectrum at lower energies to the double-hump structure at high energies, which is a direct signature of baryon stopping dynamics (Schäfer et al., 2021).

Subsequent work broadened the scope in two directions. First, Bayesian inference was applied to the $3+1$D hybrid with SMASH-generated initial conditions, using Au+Au collisions at $\sqrt{s_{NN}}=4.3$0, $\sqrt{s_{NN}}=4.3$1, and $\sqrt{s_{NN}}=4.3$2 GeV to constrain interface parameters and the temperature and baryochemical-potential dependence of shear and bulk viscosity (Götz et al., 13 Mar 2025). Second, the same hybrid logic was pushed to lower beam energies, $\sqrt{s_{NN}}=4.3$3 GeV, by replacing sudden iso-$\sqrt{s_{NN}}=4.3$4 fluidization with dynamic initial conditions sourced continuously from the hadronic transport stage (Góes-Hirayama et al., 25 Jul 2025).

At the opposite end of system size, the framework was used for O–O and Ne–Ne collisions at $\sqrt{s_{NN}}=4.3$5 TeV as a test of whether a heavy-ion-style hybrid description remains meaningful in intermediate small systems. In that setting it was compared on an equal basis to a pure SMASH hadronic cascade and to Angantyr as a baseline without collective effects (Constantin et al., 6 Sep 2025).

2. Modular architecture and switching structure

The hybrid consists of three main modules: SMASH, vHLLE, and a SMASH hadron sampler. SMASH is used twice: first as an initial-state and pre-equilibrium generator, later as the hadronic afterburner. vHLLE is the $\sqrt{s_{NN}}=4.3$6D viscous hydrodynamic solver. The hadron sampler converts fluid cells on the freeze-out hypersurface into particles consistent with the SMASH hadronic spectrum and equation of state (Schäfer et al., 2021).

Before the stage-by-stage description, the overall division of labor is compactly summarized as follows.

Stage	Code	Role
Pre-equilibrium	SMASH	Hadronic/string transport and coarse graining
Dense stage	vHLLE	$\sqrt{s_{NN}}=4.3$7D viscous hydrodynamic evolution
Particlization	SMASH hadron sampler	Cooper–Frye conversion on switching hypersurface
Dilute stage	SMASH	Hadronic rescattering and decays

In the standard initialization used across the original beam-energy-scan studies, SMASH is evolved until a constant proper time $\sqrt{s_{NN}}=4.3$8, identified with the passing time of the two nuclei up to a lower bound, and particles crossing that iso-$\sqrt{s_{NN}}=4.3$9 hypersurface are removed and mapped onto hydrodynamic fields. The passing-time estimate is

$5.02$0

with $5.02$1 in the implementations that impose a minimum hydrodynamization time (Schäfer et al., 2021, Götz et al., 2023, Constantin et al., 16 May 2026).

The coarse graining uses the standard transport-to-hydro mapping

$5.02$2

with analogous expressions for electric charge and, when needed, strangeness currents. In the longitudinally resolved SMASH initial conditions used for $5.02$3D hydrodynamics, the smearing kernel is Gaussian and Lorentz-contracted in the longitudinal direction (Schäfer et al., 2021, Constantin et al., 16 May 2026).

The hydro-to-transport switch is defined by a constant energy-density hypersurface. In the original conservation-law study this was $5.02$4, constructed with CORNELIUS; later Bayesian and low-energy studies treated $5.02$5 as a parameter or reused a calibrated value from earlier inference (Schäfer et al., 2021, Götz et al., 13 Mar 2025, Góes-Hirayama et al., 25 Jul 2025). On that hypersurface, hadrons are sampled according to the Cooper–Frye formula,

$5.02$6

after which SMASH propagates the sampled particles microscopically until interactions cease (Schäfer et al., 2021).

A practical consequence of this architecture is that the framework can be wrapped as a complete event-by-event pipeline. The 2025 SMASH review notes a public bash-based hybrid handler that automates SMASH $5.02$7 vHLLE $5.02$8 Cooper–Frye sampling $5.02$9 SMASH, and explicitly presents the SMASH-vHLLE hybrid as part of a broader ecosystem for Bayesian analysis and low-energy dynamical fluidization (Elfner et al., 29 Aug 2025).

3. Initial-state constructions and fluidization strategies

A distinguishing feature of SMASH-vHLLE is that the initial conditions can be generated by hadronic transport itself rather than imposed through a purely parametric entropy-deposition model. In the standard SMASH initial conditions, nuclei are initialized from Woods–Saxon distributions, evolved through nucleon–nucleon collisions, resonance dynamics, and, at high energies, string excitation and fragmentation via PYTHIA. Produced hadrons propagate with formation times and energy-dependent interaction rules, so the resulting pre-hydrodynamic state already contains nontrivial transverse flow, longitudinal structure, baryon stopping, and event-by-event fluctuations (Götz et al., 2023, Constantin et al., 16 May 2026).

This transport-based initialization is fully three-dimensional. In the longitudinally resolved formulation used for comparisons to saturation-based initial conditions, the deposited energy–momentum and conserved charges are functions of $\sqrt{s_{NN}}=5.36$0 at $\sqrt{s_{NN}}=5.36$1, generated by an anisotropic Gaussian kernel with transverse and longitudinal widths $\sqrt{s_{NN}}=5.36$2 and $\sqrt{s_{NN}}=5.36$3. That makes SMASH-compatible initialization directly usable by $\sqrt{s_{NN}}=5.36$4D hydrodynamics and avoids a boost-invariant ansatz at beam-energy-scan energies (Constantin et al., 16 May 2026).

The framework is also modular enough to accept alternative initial-condition models. A direct comparison among SMASH IC, Trento, and IP-Glasma within the same SMASH-vHLLE evolution showed that average eccentricities can be similar while the distributions of eccentricities, the correlations among initial-state properties, and the correlations between initial-state and final-state properties differ substantially. A central result of that study is that initial momentum anisotropy itself does not significantly determine the final $\sqrt{s_{NN}}=5.36$5 or $\sqrt{s_{NN}}=5.36$6, whereas initial radial flow or initial transverse momentum $\sqrt{s_{NN}}=5.36$7 systematically improves the prediction of final-state flow beyond $\sqrt{s_{NN}}=5.36$8 alone (Götz et al., 2023).

At low beam energies, the assumption of a sudden global initialization at one proper time becomes increasingly problematic. The dynamic-fluidization extension therefore replaces the iso-$\sqrt{s_{NN}}=5.36$9 switch by a local core–corona separation based on the local rest-frame energy density during the pre-equilibrium SMASH evolution. The coarse-grained tensor and currents are evaluated continuously on a $\rightarrow$0D lattice, and a particle is fluidized only if the local energy density, after subtracting its own self-contribution, exceeds a threshold $\rightarrow$1. Additional constraints are imposed for string products, leading nucleons, and core–corona interactions, and source terms are deposited continuously into vHLLE as the collision evolves (Góes-Hirayama et al., 25 Jul 2025).

This low-energy construction is not merely a technical refinement. It reflects the fact that at $\rightarrow$2 GeV some parts of the system may form a dense core while others remain corona-like throughout the event. The resulting dynamic initial conditions allowed the hybrid to achieve good agreement with measured bulk observables between $\rightarrow$3 and $\rightarrow$4 GeV, whereas the older iso-$\rightarrow$5 scheme becomes unreliable toward the lower end of that range (Góes-Hirayama et al., 25 Jul 2025).

At high energies, however, the validity of SMASH as an initial-state model becomes more restricted. A direct comparison to the saturation-based McDipper model found that SMASH initial conditions are best constrained for low to intermediate collision energies, particularly up to around $\rightarrow$6 GeV, with caution needed up to $\rightarrow$7 GeV. Above that, SMASH underestimates transverse-energy deposition and exhibits an unphysical plateau in midrapidity net-baryon density, so saturation-based models are argued to be more appropriate for top RHIC and LHC initial conditions (Constantin et al., 16 May 2026).

4. Hydrodynamic sector, equation of state, and transport coefficients

In all applications the intermediate stage is evolved with vHLLE, a $\rightarrow$8D viscous hydrodynamics code that solves

$\rightarrow$9

with an Israel–Stewart-type second-order viscous formulation in the full version of the code (Schäfer et al., 2021, Götz et al., 2023). Depending on the study, the hydrodynamic equation of state is a chiral model fitted to lattice QCD data, a parity-doublet mean-field equation of state matched to a hadron resonance gas, or a tabulated holographic equation of state introduced into the vHLLE-SMASH workflow (Schäfer et al., 2021, Götz et al., 13 Mar 2025, Anufriev et al., 3 Oct 2025).

A recurrent theme in the SMASH-vHLLE literature is that transport-coefficient inference is highly sensitive to how one parameterizes viscosity across the QCD phase diagram. One line of work therefore replaced constant $\rightarrow$0 by an energy-density and net-baryon-density dependent parametrization, rather than a purely temperature-dependent one, so that the effective shear viscosity is expressed directly in terms of hydrodynamic variables and less tied to a specific equation of state. In that formulation, the generalized $\rightarrow$1 substantially reduces the sensitivity of final $\rightarrow$2 to the hydro–transport switching point, and the explicit $\rightarrow$3 dependence of elliptic flow was found to be negligible and relevant only in the early stages of the collision (Götz et al., 2022, Götz et al., 2023).

The same studies also showed that anisotropic flow in the hybrid is generated by all three sectors: initial conditions, viscous hydrodynamics, and hadronic transport. A pure eccentricity-response picture is therefore incomplete. In particular, inclusion of initial radial flow in a linear response fit improves the prediction of final-state flow from initial-state properties, while initial momentum anisotropy itself adds little explanatory power once eccentricity and initial radial flow are accounted for (Götz et al., 2023).

A more global statistical-learning program was carried out in the 2025 Bayesian analysis with SMASH-generated initial conditions. That study inferred a posterior distribution preferring a near-vanishing specific shear viscosity in the high-temperature QGP phase, combined with moderate-to-large bulk viscosity around the phase-transition region, while the constraints on the baryochemical-potential dependence were weak. It also found that hadronic initial conditions constrain the evolution more strictly at intermediate energies, making parameters such as the hydrodynamic onset time highly sensitive (Götz et al., 13 Mar 2025). This result differs substantially from previous Bayesian analyses based on more parametric initial-condition models; a plausible implication is that extracted transport coefficients remain strongly model-dependent when the initial-state sector is changed.

In small systems, the same hydrodynamic machinery was used to test the onset of collectivity in O–O and Ne–Ne at $\rightarrow$4 TeV. There the baseline calculation employed temperature-dependent $\rightarrow$5 and $\rightarrow$6 from a previous Bayesian tuning, alongside a comparison run with constant $\rightarrow$7. The hybrid response was shown to be highly sensitive to the viscosity choice: the constant-$\rightarrow$8 run produced systematically smaller $\rightarrow$9 than the baseline with very low shear viscosity in the hottest cells (Constantin et al., 6 Sep 2025).

5. Particlization, hadronic microphysics, and conservation constraints

The hydro-to-particle conversion is one of the technically most constrained parts of the framework because it is where macroscopic fields must be mapped back to a discrete hadronic ensemble without introducing systematic violations of energy and conserved charges. The SMASH-vHLLE program therefore developed a dedicated SMASH hadron resonance gas equation of state for particlization, built from the same hadronic degrees of freedom used in SMASH transport, and tabulated as the map

$\mathrm dN/\mathrm dy|_{y=0}$0

under strangeness neutrality (Schäfer et al., 2021).

This consistency turns out to be essential. In the original conservation-law validation, using a raw, unmodified solver-based particlization equation of state led to severe violations of energy, baryon number, and electric charge conservation up to approximately $\mathrm dN/\mathrm dy|_{y=0}$1 after sampling, especially at lower collision energies where freeze-out cells populate low-energy-density and high-$\mathrm dN/\mathrm dy|_{y=0}$2 regions. Replacing it with an improved SMASH-consistent equation of state removed those large jumps and reduced total deviations over the whole hybrid evolution to values not exceeding about $\mathrm dN/\mathrm dy|_{y=0}$3, with the remaining mismatch attributable to known numerical limitations such as grid resolution, Milne-coordinate source terms, and Monte Carlo sampling fluctuations (Schäfer et al., 2021).

The actual particlization algorithm is grand canonical. Surface elements are treated independently, species multiplicities are drawn from Poisson distributions, and momenta are sampled from the Cooper–Frye distribution with viscous corrections. Exact event-by-event global conservation is therefore not imposed by the sampler itself, although SMASH again enforces exact conservation once the particles enter the hadronic afterburner (Schäfer et al., 2021, Götz et al., 13 Mar 2025).

The afterburner stage solves the relativistic Boltzmann equation microscopically with a broad hadron resonance spectrum. SMASH includes mesons, baryons, and resonances up to $\mathrm dN/\mathrm dy|_{y=0}$4 GeV, with elastic and inelastic scatterings, resonance formation and decay, and string dynamics at higher energies. In the hybrid context this stage is responsible for hadronic rescattering, kinetic decoupling, resonance feeddown, and late-stage modifications of spectra and correlations (Schäfer et al., 2021, Götz et al., 2023).

Because the afterburner retains the full transport microphysics, the hybrid can also be extended beyond single-hadron bulk observables. The 2025 SMASH review highlights light-cluster production via coalescence or stochastic multiparticle reactions, the propagation of critical fluctuations through hadronic rescattering, and electromagnetic probes such as dileptons with collisional broadening. It also notes that hadronic rescattering dynamics are crucial for reproducing the transverse-momentum spectra of light nuclei in hybrid calculations (Elfner et al., 29 Aug 2025).

6. Phenomenology, validity domain, and debated regimes

Across the beam-energy-scan range, the hybrid established itself primarily through its description of identified-particle production and baryon stopping. For Au+Au/Pb+Pb collisions between $\mathrm dN/\mathrm dy|_{y=0}$5 and $\mathrm dN/\mathrm dy|_{y=0}$6 GeV, it reproduced rapidity distributions, transverse-mass spectra, excitation functions of $\mathrm dN/\mathrm dy|_{y=0}$7 and $\mathrm dN/\mathrm dy|_{y=0}$8, and a reasonable centrality and energy dependence of charged-particle $\mathrm dN/\mathrm dy|_{y=0}$9. The proton rapidity evolution from a single central peak at low energy to a double-hump structure at higher energy is one of its clearest successes (Schäfer et al., 2021).

At still lower beam energies, dynamic initial conditions from hadronic transport extended the framework to $\langle p_T\rangle$0 GeV and yielded good agreement with measured bulk observables. At the same time, that low-energy program clarified an important limitation: because the hydrodynamic stage does not include spectator shadowing, the hybrid produces positive elliptic flow where the data exhibit negative squeeze-out at the lowest energies. This suggests that some observables in the very low-energy regime remain strongly sensitive to mean fields and spectator dynamics beyond the present hydro-plus-cascade baseline (Góes-Hirayama et al., 25 Jul 2025).

A related case study comes from light nuclei at $\langle p_T\rangle$1 GeV. There the SMASH+vHLLE hybrid improved deuteron yields and described femtoscopic correlation functions well when combined with coalescence and CRAB analyses, but it failed to reproduce the negative elliptic flow of light nuclei. A pure SMASH calculation with a hard Skyrme equation of state reproduced the sign and magnitude of the measured $\langle p_T\rangle$2 much better, indicating that at such low beam energy the mean-field transport sector can remain more appropriate than a viscous-fluid intermediate stage for some anisotropic observables (Bailung et al., 2024).

In small systems the central question is different: not whether transport or hydrodynamics dominates globally, but whether a hydrodynamic stage is justified at all. In O–O and Ne–Ne at $\langle p_T\rangle$3 TeV, the hybrid predicts significant radial flow, $\langle p_T\rangle$4 for $\langle p_T\rangle$5 GeV with baryon enhancement over mesons, nonzero $\langle p_T\rangle$6 and $\langle p_T\rangle$7, negative four-particle cumulants, and sensitivity of flow to $\langle p_T\rangle$8-clustered nuclear structure. By contrast, pure SMASH and Angantyr show essentially vanishing four-particle cumulants and non-flow-dominated $\langle p_T\rangle$9. The same study, however, also evaluates an opacity measure and finds that hydrodynamics is reasonably justified mainly for the top $v_2$0 most central O–O and Ne–Ne events, with applicability deteriorating toward peripheral collisions (Constantin et al., 6 Sep 2025).

At high collider energies another debate concerns the appropriateness of SMASH-based initial conditions themselves. The 2026 comparison to McDipper concludes that SMASH and saturation-based initial conditions are comparable around $v_2$1 GeV, but that SMASH underestimates transverse-energy deposition and overestimates baryon stopping as one approaches and exceeds $v_2$2 GeV. This suggests a practical domain split: SMASH-initialized vHLLE is best constrained for low to intermediate energies, whereas top RHIC and LHC applications should preferentially use saturation-based initial-state models if baryon stopping and longitudinal energy deposition are central observables (Constantin et al., 16 May 2026).

The framework is also being generalized rather than frozen in its original form. Holographic equations of state calibrated to lattice QCD by machine-learning procedures have been implemented in vHLLE-SMASH workflows, with SMASH-generated initial conditions and Hadron Sampler+SMASH freeze-out chains, and they yield hadronic spectra comparable to those from standard reference equations of state in first tests (Anufriev et al., 3 Oct 2025, Anufriev et al., 31 Oct 2025). This suggests that the SMASH-vHLLE architecture is best viewed not as a single fixed model, but as a modular baseline in which the initial state, the dense-matter equation of state, and the switching strategy can all be varied while preserving a common transport–hydrodynamics–transport logic.