NuclearConfectionery: Multi-stage Simulation Framework for Modeling Relativistic Heavy-ion Collisions

TL;DR

NuclearConfectionery delivers a GPU-accelerated, open-source, multi-stage framework for relativistic heavy-ion collisions, unifying 3+1D SPH hydrodynamics with a 4D BSQ equation of state, dynamical jet coupling, and exact BSQ-conserving particlization followed by a hadronic afterburner. The core innovation is CCAKE 2.0, which supports Israel–Stewart, DNMR, or ADNH hydrodynamics in Cartesian or hyperbolic coordinates, with offline 4D EoS inversion and BSQ diffusion, enabling consistent simulations across RHIC BES to LHC energies. The authors validate the code against semi-analytical solutions (e.g., Gubser and Landau–Khalatnikov) and perform comprehensive convergence studies, demonstrating robustness and scalability on CPUs and GPUs. This framework provides a flexible, high-performance platform for precision QCD matter studies across finite and vanishing densities, with extensive options for observables, initial-state models, and jet interactions, paving the way for systematic Bayesian parameterizations and new explorables such as critical-point behavior.

Abstract

We present the NuclearConfectionery, a modular framework for simulating the full dynamical evolution of relativistic heavy-ion collisions. Its core hydrodynamic module, CCAKE 2.0, represents a major advance over previous SPH-based relativistic hydrodynamic codes. CCAKE 2.0 simultaneously evolves energy-momentum and multiple conserved charges (B, S, Q) with a four-dimensional equation of state, and can be run in either Cartesian or hyperbolic coordinates, enabling consistent simulations from the RHIC Beam Energy Scan to LHC energies. We have implemented a particlization module that supports global BSQ charge conservation on the freeze-out surface; the resulting hadron ensemble is then propagated through a hadronic transport afterburner. A source term is included in the equations of motion to couple jets to the fluid, allowing simultaneous bulk and hard-probe evolution or, alternatively, for stopped baryons at low beam energies. The framework offers flexible choices of equations of motion (Israel-Stewart, DNMR, ADNH) and transport coefficients, along with GPU-ready performance via Kokkos/Cabana, offline equation of state inversion for 4D tables, and containerized portability. We validate the code with semi-analytical benchmarks (including BSQ Gubser and Landau-Khalatnikov solutions) and extensive convergence studies. The NuclearConfectionery provides a user-friendly, high-performance, open-source tool for event-by-event simulations across collision energies, offering flexibility to study QCD matter at both vanishing and finite densities.

Figures (34)

• Figure 1: Heavy-ion collisions dynamical simulation framework. The boxes in white are modules that are always run whereas the boxes in gray are optional. Each initial state is an event that is indexed using $\left\{\mathrm{ev}\right\}$ and run through the entire simulation chain, where the final output is a set of observables $\vec{\mathcal{O}}_{\left\{\mathrm{ev}\right\}}$ for that particular event. The time scales, energy density, or temperatures that are relevant characteristic scales used to switch between modules are indicated above and explained in the text. Note that the particlization and observables modules have no internal dynamics such that they do not have a characteristic scale. Since this is only for a single event, the Analysis Script module is not shown.
• Figure 2: Hybrid simulation framework: Initial Conditions, Pre-hydrodynamics (optional in pink), hydrodynamics (requires EoS and EoM), Particlization, Afterburner, and Event analysis. One option is to run Jets (in pink) coupled to the hydrodynamics phase that is then fed into Event Analysis afterwards. Example YAML files and the simulation chain can be found at https://github.com/the-nuclear-confectionery/confectionery_chain.
• Figure 3: Passing time for two heavy nuclei with an estimated total radius of $R_0=6$ fm across different center-of-mass beam energies $\sqrt{s_{NN}}$.
• Figure 4: Illustration comparing c c a k e runtime for the same a m p t initial condition. The a m p t event is initialized in a similar way using either a hyperbolic (orange) or a Cartesian coordinate system (blue).
• Figure 5: Comparing the observables for the aforementioned a m p t event, with the same initialization methods described before. Each initial energy density was normalized to match experimental data from star STAR:2017sal.