HYDJET++: Hybrid Monte-Carlo Event Generator

• HYDJET++ is a hybrid Monte-Carlo event generator that models heavy-ion collisions by combining a soft hydro-inspired state with a hard multi-parton fragmentation component.
• The soft sector uses thermal freeze-out and explicit resonance decays via FAST MC, while the hard sector employs PYTHIA and PYQUEN to simulate jet production and energy loss.
• Designed for a wide energy range from FAIR to LHC, HYDJET++ provides insights into bulk observables, flow harmonics, and event fluctuations despite approximations in charge conservation and higher-order dynamics.

HYDJET++ is a Monte-Carlo event generator for simulation of relativistic heavy ion $AA$ collisions considered as a superposition of the soft, hydro-type state and the hard state resulting from multi-parton fragmentation. It was introduced as the development and continuation of HYDJET, with the main program written in object-oriented C++ under the ROOT environment, the hard part retained from the Fortran-written HYDJET, and the soft part implemented through an adaptation of FAST MC. In this formulation, HYDJET++ is intended as a fast, hybrid framework for studying multiparticle production over a wide energy range, from FAIR and NICA to RHIC and LHC, while retaining simultaneous access to bulk collective observables and hard-probe phenomena (0809.2708).

1. Event concept and computational architecture

The defining structural assumption of HYDJET++ is the independent generation and subsequent superposition of two components: a soft hydro-type hadronic state and a hard state from multi-parton production, in-medium energy loss, and fragmentation. This separation is not merely organizational. It is the central modeling device by which low-$p_T$ bulk production, freeze-out, flow, and resonance feed-down are treated together with high-$p_T$ jets, quenching, and shadowing in the same event record (0809.2708).

The soft sector is generated with FAST MC and represents a thermalized hadronic source emitted from chemical and thermal freeze-out hypersurfaces. The hard sector is built from PYTHIA-generated partons propagated through the medium with PYQUEN, followed by hadronization through the Lund string model. Later phenomenological applications preserve this same division: in particular, analyses of charm flow at RHIC and of Xe–Xe observables at the LHC continue to use HYDJET++ as a two-component or “hydrodynamics plus jets” framework, rather than as a full dynamical hydrodynamic simulation (Saxena et al., 14 Jun 2025, Pandey et al., 6 May 2025).

This architecture supports several practical switches and extensions. The generator can toggle jets, quenching, the soft sector, and nuclear shadowing; it stores event output in ROOT trees with particle kinematics, emission coordinates, genealogy, and origin labels. The soft hadron table includes more than 360 meson and baryon states, including charmed states, which is one reason the model has been repeatedly adapted to identified-particle, heavy-flavor, and fluctuation studies (Bravina et al., 2016).

2. Soft sector: thermal emission, freeze-out, and hadrochemistry

The soft component of HYDJET++ is a hydro-inspired freeze-out model rather than an event-by-event solution of relativistic viscous hydrodynamics. Hadrons are produced on chemical and thermal freeze-out hypersurfaces with longitudinal, radial, and elliptic flow, and resonance decays are included explicitly. The equilibrium distribution used for hadron species $i$ is written as

$$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$

with the corresponding equilibrium density $\rho_i^{\rm eq}(T,\mu_i)$, an effective thermal volume $V_{\rm eff}$, and Poissonian multiplicity fluctuations around the mean $\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$ (0809.2708).

A key design feature is the distinction between chemical and thermal freeze-out. HYDJET++ allows $T^{\rm ch} \ge T^{\rm th}$, and the chemically frozen evolution between them is implemented by preserving particle ratios through effective thermal chemical potentials. This structure has been used across disparate beam energies. In the Au–Au Beam Energy Scan analysis, the chemical freeze-out parameters were taken from the Cleymans–Redlich parameterization,

$$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$

with the paper emphasizing that lower beam energies produce larger $p_T$0 and consequently stronger baryon–antibaryon asymmetry (Nayak et al., 2024).

The same formal freeze-out machinery is routinely specialized to particular hadron species. For prompt $p_T$1-mesons in PbPb at $p_T$2 TeV, the default inclusive-hadron freeze-out description was sufficient, whereas $p_T$3 required early thermal freeze-out near chemical decoupling with reduced collective velocities (Eyyubova et al., 2015). In Beam Energy Scan studies of $p_T$4 and $p_T$5, earlier freeze-out hypersurfaces were introduced for multi-strange hadrons, motivated by their small hadronic interaction cross sections and, for the $p_T$6, OZI-suppressed coupling to non-strange hadrons (Devi et al., 21 Aug 2025). In Pb+Pb studies of pure multi-strange hadrons at the LHC, species-dependent $p_T$7 and centrality-dependent $p_T$8 were used instead of a universal freeze-out prescription (Devi et al., 2024).

3. Hard sector: jets, quenching, and the soft–hard separation scale

The hard component of HYDJET++ is inherited from HYDJET/PYQUEN and models hard $p_T$9 sub-collisions, in-medium rescattering, radiative and collisional energy loss, and hadronization. Initial hard events are generated with PYTHIA 6.4, hard production vertices are sampled according to nuclear geometry, and the final partons are fragmented with the Lund model after quenching (0809.2708).

The mean number of hard sub-collisions is determined from the hard $p_T$0 cross section above a minimum scale $p_T$1, combined with the nuclear overlap and shadowing factors. This threshold has a substantive physical meaning throughout later applications: partons produced below $p_T$2 are treated as thermalized and folded into the soft component. That convention is central, for example, to RHIC charm-flow studies where the model is explicitly described as combining hydrodynamic-like soft production with jet-quenching-driven hard production (Saxena et al., 14 Jun 2025).

Energy loss in the hard sector is modeled as the cumulative effect of repeated scatterings in an expanding QGP. HYDJET++ includes collisional loss in the high-momentum-transfer approximation and medium-induced radiative loss in the BDMS framework; for heavy quarks it implements dead-cone suppression of small-angle gluon radiation. Later charm applications preserve this same structure and emphasize that both collisional and radiative energy loss, realistic Glauber geometry, and impact-parameter-dependent shadowing are retained when moving from inclusive-hadron phenomenology to heavy flavor (Eyyubova et al., 2015).

The model’s hybrid logic becomes most transparent at intermediate transverse momentum. In multiple analyses, the fall-off of $p_T$3, the degradation of NCQ scaling at LHC energies, and the appearance of hump-like $p_T$4 dependences are all attributed to the changing relative weight of soft collective hadrons and jet fragments. This suggests that HYDJET++ should be read not as a purely hydrodynamic or purely transport model, but as a controlled superposition model in which the observable consequence of the soft–hard crossover is itself a central result (Bravina et al., 2016).

4. Flow generation, higher harmonics, and event-by-event fluctuations

HYDJET++ uses the standard Fourier decomposition

$p_T$5

with $p_T$6, $p_T$7, and $p_T$8 interpreted as elliptic, triangular, and quadrangular flow. In the model, anisotropic flow is not obtained from a full space-time hydrodynamic evolution. Instead, it is encoded directly in freeze-out geometry and flow-field parameterizations. Elliptic flow is controlled by a spatial anisotropy $p_T$9 and a momentum anisotropy $i$0, with the frequently used relation

$i$1

Triangular and higher structures are introduced through additional deformations such as $i$2, $i$3, and $i$4 (Bravina et al., 2015, Saxena et al., 14 Jun 2025).

A major line of HYDJET++ work concerns the emergence of higher harmonics through nonlinear mode coupling rather than independent primordial eccentricities. In PbPb at $i$5 TeV, the cross-talk of $i$6 and $i$7 was used to explain the $i$8- and centrality-dependence of $i$9, the basic trends of $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$0 and $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$1, and the ridge structure in long-range dihadron angular correlations. In that formulation, $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$2 is generated by the interference of $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$3 and $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$4, while $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$5 receives contributions from both sectors, with the approximate relation

$$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$6

explicitly emphasized (Bravina et al., 2013, Bravina et al., 2016).

HYDJET++ has also been used to study event-by-event flow fluctuations and factorization breaking. In the hydro-inspired freeze-out picture, the model already generates intrinsic fluctuations from momenta, multiplicities, resonance decays, mini-jets, and impact-parameter smearing. With unfolding procedures matched to ATLAS analyses, the default model was found to produce $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$7 and $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$8 distributions that were too narrow, which led to a minimal extension in which $$f_{i}^{\rm eq}(p^{*0};T^{\rm ch},\mu_i,\gamma_s) = \frac{g_i}{\gamma_s^{-n_i^s}\exp([p^{*0}-\mu_i]/T^{\rm ch}) \pm 1},$$9 and $\rho_i^{\rm eq}(T,\mu_i)$0 are smeared event by event. That modification allowed the model to reproduce unfolded $\rho_i^{\rm eq}(T,\mu_i)$1 and $\rho_i^{\rm eq}(T,\mu_i)$2 distributions and the corresponding eccentricity fluctuations much more accurately (Bravina et al., 2015).

Sub-leading flow modes have been analyzed with principal component analysis. In PbPb at $\rho_i^{\rm eq}(T,\mu_i)$3 TeV, HYDJET++ reproduced the leading modes $\rho_i^{\rm eq}(T,\mu_i)$4 and $\rho_i^{\rm eq}(T,\mu_i)$5 rather well, while the sub-leading elliptic mode $\rho_i^{\rm eq}(T,\mu_i)$6 was found to be small but nonzero and increasing toward peripheral collisions, and the sub-leading triangular mode $\rho_i^{\rm eq}(T,\mu_i)$7 even smaller and nearly centrality independent. These results supported the interpretation of nonzero $\rho_i^{\rm eq}(T,\mu_i)$8 PCA modes as fluctuation-driven factorization breaking rather than artifacts of a purely averaged flow picture (Cirkovic et al., 2016).

5. Phenomenological applications across systems and energies

HYDJET++ has been used across a broad range of colliding systems, deformation scenarios, and beam energies. In PbPb collisions at LHC energies, it has been applied to bulk flow, higher harmonics, dihadron ridge formation, event-by-event flow fluctuations, charmed hadrons, multi-strange hadrons, charge balance functions, and net-charge fluctuations (Bravina et al., 2016, Eyyubova et al., 2015, Devi et al., 2024).

In Xe–Xe collisions at $\rho_i^{\rm eq}(T,\mu_i)$9 TeV, the model has been extended to deformed nuclei and studied in body-body and tip-tip configurations. For charged hadrons, HYDJET++ was used to calculate $V_{\rm eff}$0, $V_{\rm eff}$1 spectra, and $V_{\rm eff}$2, with minimum-bias results at midrapidity matching ALICE and CMS data reasonably well and generally tracking the data more closely than AMPT string melting. In the harmonic-flow analysis up to $V_{\rm eff}$3, HYDJET++ reproduced the overall centrality dependence of $V_{\rm eff}$4, $V_{\rm eff}$5, and $V_{\rm eff}$6, observed mass ordering for pions, kaons, and protons at low $V_{\rm eff}$7, and captured the qualitative structure of $V_{\rm eff}$8–$V_{\rm eff}$9, $\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$0–$\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$1, and $\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$2–$\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$3 correlations, while also showing that body-body results are generally higher than tip-tip results (Pandey et al., 6 May 2025, Pandey et al., 2021).

At RHIC energies, HYDJET++ has been used both for Beam Energy Scan systematics and for heavy flavor. In Au–Au collisions from $\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$4 to $\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$5 GeV, the model described identified-hadron ratios, low-$\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$6 spectra, and elliptic flow using freeze-out parameters tied to the Cleymans–Redlich parameterization and a scaling relation $\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$7 between initial eccentricity and freeze-out anisotropy. The same study reported that the invariant-yield ratio of central and peripheral collisions is independent of beam energy and that strange hadrons require smaller $\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$8 than non-strange hadrons (Nayak et al., 2024). For open-charm hadrons in Au+Au at $\bar N_i = \rho_i^{\rm eq} V_{\rm eff}$9 GeV, HYDJET++ reproduced STAR $T^{\rm ch} \ge T^{\rm th}$0 $T^{\rm ch} \ge T^{\rm th}$1 and $T^{\rm ch} \ge T^{\rm th}$2 data up to $T^{\rm ch} \ge T^{\rm th}$3 GeV/$T^{\rm ch} \ge T^{\rm th}$4, as well as NCQ scaling trends, mass ordering, and baryon–meson grouping, thereby supporting substantial charm-medium coupling and partial thermalization (Saxena et al., 14 Jun 2025).

The model has also been adapted to deformed U+U collisions at $T^{\rm ch} \ge T^{\rm th}$5 GeV. In that setting, the default Woods–Saxon density was replaced by a deformed profile, and the freeze-out conditions were made centrality dependent. The resulting simulations supported only a small correlation between charged-hadron observables and the initial geometric orientation, consistent with the experimental observation that tip-tip and body-body separation in U+U is weaker than some earlier models predicted (Singh et al., 2017).

6. Known limitations, model tensions, and major modifications

A recurring point in the HYDJET++ literature is that the model is hydro-inspired, not a full time-dependent viscous hydrodynamic simulation. This is a modeling choice rather than a defect, but it constrains what can be expected from the generator. Several studies explicitly note that its Bjorken boost-invariant hydrodynamic implementation works best near midrapidity and becomes less reliable at larger rapidity or away from the kinematic domain where the soft component dominates (Pandey et al., 6 May 2025, Singh et al., 2017).

Some discrepancies are structural. In the PbPb harmonic-correlation study against ATLAS, HYDJET++ reproduced the slope of the $T^{\rm ch} \ge T^{\rm th}$6–$T^{\rm ch} \ge T^{\rm th}$7 correlation but failed for $T^{\rm ch} \ge T^{\rm th}$8–$T^{\rm ch} \ge T^{\rm th}$9 and $$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$0–$$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$1, predicting slopes that were too steep and missing the boomerang-like structure observed in data. The paper explicitly connected this to the fact that fourth-order anisotropy is not explicitly implemented in the configuration used there, which limits the model’s ability to describe $$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$2-related observables (Dordevic et al., 2019).

Charge-sensitive observables exposed another important limitation of the default implementation. In the soft sector, hadron yields are generated in a grand canonical ensemble, so electric charge is conserved only on average. This is adequate for many single-particle observables, but it fails for the low-$$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$3 charge balance function and for net-charge fluctuations. Two later studies therefore introduced a modified version of HYDJET++ with explicit event-by-event charge conservation in the soft sector: half of the charged hadrons are discarded, each remaining charged hadron is paired with an opposite-charge counterpart, and the partner’s $$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$4 and $$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$5 are sampled from Gaussian distributions around the original particle. This canonical-like local balancing substantially improved the description of balance-function widths and of the strongly intensive fluctuation measures $$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$6 and $$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$7 in comparison with ALICE and CMS data (Chernyshov et al., 2022, Ambaryan et al., 2024).

Other tensions arise in species-dependent phenomenology. For pure multi-strange hadrons in Pb+Pb at the LHC, HYDJET++ provides a reasonable description in central and semi-central collisions, but deviations in peripheral bins and at higher $$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$8 were interpreted as evidence for missing physics such as coalescence and more detailed jet–medium coupling. In RHIC Beam Energy Scan studies of $$T(\mu_B)=a-b\mu_B^2-c\mu_B^4,\qquad \mu_B(\sqrt{s_{NN}})=\frac{d}{1+e\sqrt{s_{NN}}},$$9, the model overpredicts $p_T$00 across energies and centralities, suggesting that $p_T$01 baryons may not achieve the same degree of thermal equilibration as $p_T$02 mesons and that additional mechanisms such as baryon coalescence are not fully captured (Devi et al., 2024, Devi et al., 21 Aug 2025).

These limitations clarify the model’s domain of validity. HYDJET++ is strongest when the problem requires a fast, unified Monte Carlo baseline that preserves the interplay of soft collective emission, resonance decays, jet quenching, and event geometry. It is less suited when the observable depends sensitively on exact local conservation laws, detailed higher-order initial-state dynamics, or hadronization channels absent from the default soft-plus-fragmentation picture. Within those bounds, it remains one of the most widely repurposed hybrid generators for heavy-ion phenomenology from the RHIC Beam Energy Scan to LHC PbPb and Xe–Xe collisions (0809.2708, Bravina et al., 2016).