Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions

The microscopic description of heavy-ion reactions at low beam energies is achieved within hadronic transport approaches. In this article a new approach SMASH (Simulating Many Accelerated Strongly-interacting Hadrons) is introduced and applied to study the production of non-strange particles in heavy-ion reactions at $E_{\rm kin}=0.4-2A$ GeV. First, the model is described including details about the collision criterion, the initial conditions and the resonance formation and decays. To validate the approach, equilibrium properties such as detailed balance are presented and the results are compared to experimental data for elementary cross sections. Finally results for pion and proton production in C+C and Au+Au collisions is confronted with HADES and FOPI data. Predictions for particle production in $\pi+A$ collisions are made.


I. INTRODUCTION
Heavy-ion collisions offer the opportunity to study hot and dense strongly interacting matter under extreme conditions. High energy programs at the Large Hadron Collider (LHC) and the Relativistic Heavy Ion Collider (RHIC) are delivering a lot of detailed experimental data [1][2][3] relevant for the high temperature and low net baryochemical potential part of the phase diagram which corresponds to the situation shortly after the Big Bang. Scanning the beam energies to lower values as currently done at the CERN-Super Proton Synchrotron (SPS) [4] and the RHIC beam energy scan [5][6][7] program or in the future at FAIR and NICA provides access to regions in the phase diagram where a first order transition to the quarkgluon plasma is expected to take place. One of the goals of these programs is to search for a critical endpoint in the QCD phase diagram [8].
Since there is no first principle solution of the manybody problem in quantum chromodynamics including a non-equilibrium evolution through a phase transition up to date, effective theoretical approaches are necessary to describe the full dynamical evolution of heavy-ion reactions from the early to the late stages. By comparison of the output of these calculations with experimental data on particle distributions and their correlations in the final state, it is possible to draw conclusions about the properties of the hot and dense strongly interacting matter that was created for a very short time and in a very small volume.
Following the realization that the quark-gluon plasma behaves like an almost perfect fluid in contrast to the ideal gas expectation, within recent years the community has converged towards a standard model for the description of the evolution of heavy-ion reactions at high beam energies. The early stage of the collision is described by a non-equilibrium evolution likely based on fluctuating color fields/strings until approximate local equilibrium is reached [9,10]. The hot and dense stage of the evolution is governed by relativistic dissipative hydrodynamics [11][12][13][14] incorporating the QCD equation of state provided by lattice calculations [15][16][17][18]. The later dilute stages are described by a hadron transport approach [19]. Even though most of the dynamical features are captured within the hydrodynamic calculation, the hadronic rescattering stage becomes necessary as soon as one wants to address identified particle spectra or correlation and fluctuation observables that are affected by resonance decays and baryon annihilation [20,21].
The other limit where the description of the dynamical evolution of heavy-ion reactions is to some degree under control is at very low beam energies that are dominated by hadronic reactions and not yet affected by quark-gluon plasma formation. The region of intermediate beam energies that is of great interest with respect to the discovery of features in the QCD phase diagram poses a challenge to the current dynamical approaches. There are attempts to adapt the above described hybrid approaches and extend them to finite baryo-chemical potential [22]. The other option is to start from a vacuum hadronic transport approach that is extensible by including effects of the hot arXiv:1606.06642v2 [nucl-th] 17 Jan 2017 and dense medium such as many-body interactions. This second approach is the motivation for the development of a new hadronic transport approach, SMASH.
Hadronic transport approaches have been developed for 20-30 years and some models are still under active development [23][24][25][26]. The new experimental data that is available to constrain the resonance properties at low beam energies [27] and profiting from the experience of the existing transport approaches is the reason for developing a modern flexible open source code that can be adapted as a standard reference for a purely hadronic system with vacuum properties. To summarize, we have gained a lot of new experimental and theoretical insights over the past two decades that make the development of a new transport approach a timely endeavor. This new transport approach will also be highly relevant to provide a better understanding of the late stage evolution of hadronic rescattering at RHIC and LHC energies.
In this paper the newly developed approach is described in detail. In Section II the ingredients of the approach are explained including the general setup, the collision criterion, the initial conditions and treatment of potentials, Pauli blocking and resonance formation and decay. In Section III basic checks of detailed balance and comparisons with elementary cross sections are shown. In Section IV we present calculations of observables in comparison with experimental data from HADES and FOPI at E kin = 0.4 − 2A GeV and predictions for π − A collisions.

A. General Setup
The main advantage of a microscopic transport approach is that the full phase-space information of all particles is available at all times. SMASH constitutes a solution of the non-equilibrium dynamics of hadrons in the regime where the inelastic interactions are treated by resonance excitations and decays with vacuum properties. The underlying equation is the relativistic Boltzmann equation where C i coll is the collision term, F α is the force experienced by individual particles and m i is the particle mass. For high beam energy collisions, F α = 0, while for low beam energy collisions, F α = −∂ α U (x) where U (x) is the mean-field potential. The relativistic Boltzmann equation is an integro-differential equation in 6+1 dimensions. f i (x, p) is the single particle distribution for each species i that is represented by test particles. Along the lines of quantum molecular dynamics each particle is in principle represented by a Gaussian wave packet. In practice, all particles are treated as point particles and the finite spatial extent is only invoked to calculate thermodynamic properties like the particle density. In our case, per default each real particle is represented by one test particle, but more test particles can be created if necessary.

Collision Criterion
One of the major challenges for solving the Boltzmann equation in a relativistic situation is to define an appropriate collision criterion. The Kodama criterion [28] is a fully covariant collision criterion, but since it involves boosts of several four vectors it is rather inefficient. In the current approach we have chosen to use the geometrical criterion employed in the UrQMD (Ultra-relativistic Quantum Molecular Dynamics) approach [23], that is defined as follows: where r and p are the coordinates and momenta of the two particles a and b in the center of mass frame of the binary collision. The time of the collision is determined as the time of the closest approach in the computational frame: where now all coordinate and momentum vectors have to be taken in the computational frame. The computational frame is usually chosen to be the equal velocity frame of the two nuclei which is the same as the center of mass frame in case of symmetric systems. The computational system is the one that carries the clock that is relevant for ordering of the collisions, therefore it is crucial to transform the collision times to the same frame to decide which collision happens first. This geometrical criterion effectively encodes an instantaneous interaction over a finite distance and gives rise to causality violations [29]. We have compared the UrQMD criterion to the covariant Kodama criterion and found no significant differences. Since the above explained criterion is numerically more efficient, we stick to this definition in the following.
A different option to include all relevant scatterings at high density is to implement the solution of the Boltzmann equation by stochastic rates [30][31][32]. This approach has the advantage that multi-particle scatterings can be taken into account in a straightforward way. On the hadronic level there are of course a lot of different possibilities that one would need to take into account in such an approach, therefore this is left for future work. Also, the stochastic rates approach is relying on having a large number of test particles in each cell, therefore it is not clear how to model event-by-event fluctuations properly.

Test particles
Another method to circumvent the locality issues is the test particle method: all cross sections are scaled by a factor N −1 test , while the number of initially sampled particles is increased by the same factor N test .
N test is referred to as "test particle number". After substitution (Eqs. (5) and (6)) the scattering rate (number of collisions per unit time per particle) remains unchanged, but the cross sections become smaller and collisions are "more local". Locality is restored in the limit N test → ∞. As shown in [29], experimental observables such as particle spectra and flow obtained using transport models depend on N test and saturate when N test is sufficiently large (in case of [29] N test = 16 was large enough for saturation). Another important application of test particles is to provide statistics for density or phase-space density estimates, which are required for evaluating potentials accurately.

Time Steps/Propagation
To solve the Boltzmann equation numerically, time and space need to be divided into cells. The granularity of the time steps are crucial, since the time steps need to be small enough to catch all collisions (when assuming a maximum of one collision per particle in a timestep) and as large as possible to ensure a fast evaluation of the evolution. Therefore, SMASH has two different options for the propagation. Either fixed time steps are chosen or the time steps are dynamically determined from the collision times. The first setup has the advantage that calculations with nuclear potentials are feasible while the second one adapts nicely to high and low density regions and is more efficient.
In an algorithm with fixed time step size the actions are only determined for a short time dt in advance. In addition to the time step size dt a start t start and an end time t end have to be chosen. The search for collisions needs to extrapolate the movement of the particles. The assumption in SMASH is that the time step size ∆t is small enough that the effect of potentials on the trajectory of the particles can be neglected during this time interval. Therefore, the movement is extrapolated without taking potentials into account.
Any interaction of two particles is called an action. If the criterion Eq. (2) is satisfied, then the collision is added to the list of collisions and decays with a time stamp t coll . After all actions are found, they are sorted according to their associated time. Iterating over the sorted list, all actions are first tested whether they are still valid. We rely on the assumption that each particle only interacts once during one time step. Valid actions are performed which involves replacing the incoming particles with the outgoing particles. The actions are performed before the propagation, which means that the global time of the particles is still the time of the beginning of the time step. When all valid actions have been performed, all particles are propagated taking potentials into account (if they are present).
When propagating without potentials or for the later dilute stages of the collision, it is useful to abandon any fixed time steps and rather switch to an algorithm that takes the actions themselves to determine the next propagation step. The general idea is that the program keeps a list of actions which is constantly updated. Actions are removed from the list as they are performed and added as they are newly discovered. At the beginning of the simulation, the end time of the simulation t end is specified. This is used to find all possible actions for all particles until t end as described before. Next, there is a loop over all actions that starts with the first action according to the time of execution of the actions and checks if this action is still valid. The check consists of verifying that the incoming particles were not part of another action since this action was found. If they were, the action is discarded.
If the action is valid, all particles are propagated to the point in time where the action is supposed to happen. Then the action is performed as described before. As a result, all actions that involved the incoming particles are implicitly rendered invalid, since the state of these particles has changed. In the last step, all possible actions of the outgoing particles are added to the list of actions. This algorithm realizes all actions that are supposed to happen, assuming the time ordering of the actions is correct (which depends on the collision criterion). Fig. 1 shows the time difference between consecutive interactions ∆t for a Cu-Cu collision, averaged over 100 events. Each data point corresponds to one action. The average of this time between actions becomes as small as 0.01 fm/c but the variance of ∆t is high and outliers can reach 0.0001 fm/c.

Mean-field potentials
To create a more realistic simulation at low beam energies, a minimal version of mean-field potentials between nucleons is included. The equations of motions have to be adjusted according to the modified one-particle Hamiltonian H i where m eff is the mass for stable hadrons and the effective mass for resonances in accordance with their mass distribution (e.g. Breit-Wigner). At this point, the potential depends only on the coordinates, but not on the momentum of the particles. The corresponding equations of motion are then This formulation leads to the fact that momentum conservation is fulfilled only on average. Event by event momentum conservation requires that d p i /dt = −∂H tot /∂ r i , where H tot = i H i . The potential is calculated as a function of the local density Here ρ is the Eckart rest frame baryon density and ρ I3 is the Eckart rest frame baryon isospin density of the relative isospin projection I 3 /I. ρ 0 = 0.1681/fm 3 is the nuclear ground state density. Parameters for the Skyrme potential are by default set to a = −209.2 MeV, b = 156.4 MeV and τ = 1.35, while S pot = 18 MeV is the default value for the symmetry potential. These parameters were agreed on for a recent transport code comparison [33] and correspond to a rather soft potential with an incompressibility of K = 240 MeV. For the equations of motion one does not need the potential itself, but its gradient, ∂U/∂ r. In the symmetry term the positive sign is applied for the potential acting on neutrons and the minus sign is applied for the potential acting on protons. Currently, the potential acts only on baryons. The potentials are always calculated after the actions are performed, right when the propagation happens. We note that electromagnetic potentials (Coulomb and Lorentz force) are currently being neglected in the model, since they are typically much weaker than the hadronic mean fields (even if they are more long-ranged). The Coulomb potential can only play a role for collisions of large nuclei at very low energies and is completely negligible at higher energies (FAIR/RHIC/LHC).

Nearest neighbor search
To determine if two particles will scatter, their distance needs to be compared to the total cross section. In principle, every particle has to be paired with every other particle in the system and the complexity of the search will scale with N 2 where N is the number of particles, which is computationally intensive. In order to reduce the combinatorics of this search, the space can be divided into cells whose sizes are chosen such that, accounting for the time step size ∆t and the maximal possible cross section σ max tot , all collisions will happen within one cell or among neighboring cells, but not beyond that. SMASH uses such a grid structure with a minimal cell size of (2.5 fm) 3 . The value of d = 2.5 fm corresponds to a maximum cross section of σ max tot = πd 2 ≈ 200 mb that can be handled. This maximum is reached in the ∆ peak of the π + p cross section, see Fig. 13. The only exception of physical cross sections going above 200 mb are the elastic NN cross sections, which diverge at the threshold. Those are effectively being cut at 200 mb with our minimal cell size (at least without test particles). When larger numbers of test particles are used, the cross sections are scaled down accordingly and this limitation is lifted.
In the actual algorithm for iterating over the cells, a distinction is made between in-cell search and neighbor search. The in-cell search is used to find decays and collisions within a given cell. When searching for collisions, particle 1 is checked with particles 2 and 3, but not with particle 4.
The neighbor search looks for actions between the particles in a given cell and its neighbors. To avoid finding duplicate actions, not all neighboring cells are used in the neighbor search. Consider the case depicted in Fig. 2: When starting the neighbor search from the dark gray cell, the actions between particle 1 and particles 2 and 3 will be found. Afterwards, when starting the search from a light gray cell, the dark gray cell can be omitted from the search, because there would be no new actions. After some analysis one comes to the conclusion that each cell needs to check only half of its neighbors (except for cells at the border which need to check even fewer).
The grid is also used to realize periodic boundary conditions. With periodic boundaries, particles that are on opposite sides of a fixed-sized box can interact. When this feature is activated, the neighbor search for a cell at the border does not only check the actual neighboring cells. Instead, so-called ghost cells are added that contain the mirrored particles from the opposite side of the grid. Note that in this case it is important that no grid cell size at the boundaries is smaller than the minimal cell size. Therefore, the cell size is scaled up to fit the total volume with the minimal number of complete cells.

Elastic Box Test
To test the collision finding algorithm, we employ a simple setup, which we further call "elastic box". A box with periodic boundary conditions is uniformly filled with N pions. The momenta are distributed according to a Boltzmann distribution where the temperature T is taken to be 0.13 GeV. The pions are only allowed to scatter elastically with a constant isotropic cross section σ. In this simple setup the scattering rate s should be s = nσ, where it is taken into account that the relative velocity between particles is close to the speed of light (a calculation using Eq. (52) of [34] gives v rel ≈ 0.98 for our setup), and n is the particle density. From Fig. 3 one can see that s = nσ is fulfilled for SMASH regardless of cross section σ or test particle number N test , but only if the time step is sufficiently small. It is also interesting to observe that the scattering rate s for all combinations of different σ, N test and ∆t lies on one universal line For small time steps nσ∆t = ∆t/λ 1 one retrieves the expected ideal gas behavior, while in the limit of large time steps there is one collision per particle per time step. This is expected, because more than one collision per particle per time step is prohibited by the SMASH algorithm. If one wants to have the full collision rate, the propagation from collision to collision without timesteps n σ 22    can be used. The universality of the curve can be explained in terms of dimensions. One can construct only three dimensionless quantities from s, n, σ and ∆t. Let us choose s∆t, nσ∆t = ∆t/λ and nσ 3/2 , then s∆t = f (nσ∆t, nσ 3/2 ). Assuming independence on the second argument one obtains our universal curve. We have additionally checked that these results depend only on the density n = N/V , but not on N or V separately. In other words, one can vary number of particles N , box volume V or both, but the results are identical if n = N/V is the same.
B. Initial Conditions

Nuclear Collisions
a. Nucleon Distribution in coordinate space To generate initial conditions for heavy-ion collisions the whole phase-space distribution of the initial nucleons needs to be sampled. In coordinate space, 'round' nuclei like gold or lead can be described by Woods-Saxon distributions where d is the diffusiveness of the nucleus which controls the quick fall-off of the distribution. For d → 0, the nucleus is a hard sphere. ρ 0 = 0.168 fm −3 and r 0 are, in this limit, the nuclear ground state density and the nuclear radius. The default value for the diffusiveness is d = 0.545 fm, where more specific values are used for Au, Pb, Cu and U (see Table I).
Other values can be provided via the corresponding parameters in the configuration input file if necessary. Within the sampling procedure the finite size of the nucleons and nucleon-nucleon correlations are neglected for simplicity [35]. In Fig. 4 it is shown that the sampling in coordinate space for a lead nucleus works as expected.
The initial positions of nuclei and the time of initialization are chosen as shown on Fig. 5. We use Cartesian coordinates, where the z-direction corresponds to the beam direction and x is the impact parameter direction. At the initialization the projectile center is at xzcoordinates (b/2, −∆z − γ −1 P (R P + d P )) and the target center is at (−b/2, v T v P ∆z +γ −1 T (R T +d T )). Here R P,T are the projectile/target radii and d P,T are the corresponding diffusiveness parameters from the Woods-Saxon distribution. By v T,P we denote absolute values of the velocities, while γ P,T = (1 − v 2 P,T ) −1/2 are the associated gammafactors. The separation of the centers of the nuclei in x-direction equals the impact parameter b. For deformed (R P +d P )/γ P nuclei an additional rotation along all three angles is applied. In this way, the simulation is started at such an initial separation that the potential of one nucleus does not influence the other one yet, otherwise initialization in the ground state would not be justified. The initial coordinates and time are chosen in such a way that Lorentzcontracted spheres of radii (R + d) P,T will touch at t = 0 in a central collision. An alternative definition would be that t = 0 fm corresponds to the maximal overlap of the two nuclei. The additional distance ∆z = 2 fm is added to avoid missing any nucleon-nucleon collisions. Since the nucleons are distributed according to Woods-Saxon distributions, there is a small, but non-zero probability to position a nucleon at a large distance from the nucleus center. The initial separation distance ∆z is chosen such that all collisions are taken into account. The initial time is t 0 = ∆z/v P , which implies that the projectile is always moving, v P > 0, while the target can be at rest depending on the reference frame for the calculation.
b. Fermi Motion In momentum space nucleons optionally get Fermi momenta, then target and projectile are boosted in z direction according to the chosen energy of the reaction and computational frame. The gammafactor of the boost is γ = E A /M A , where E A is the energy of the nucleus and M A is its mass. The velocity of the boost is β = p A /E A . Note that in E A and M A one has to account for the binding energy of the nucleus. For this we adopt an approximation used in the JAM transport code [24], which assumes that all nucleons are equally bound. Thus, the energy of each nucleon in the rest frame of the nucleus is E i = M A /A, where A is the number of nucleons. With this assumption the boost of the longitudinal momenta p iz to the computational frame becomes where p beam is the beam momentum per nucleon and p iz are the momenta of nucleons in the rest frame of the nucleus. In our implementation p beam and γ themselves are computed without accounting for binding energy. We note that there is no well-established procedure of boosting nuclei with the account of their binding energy. Codes like UrQMD [23], JAM [24] and GiBUU [26] apply different methods. Though the typical binding energy per nucleon is much smaller than the nucleon mass ( 8 MeV/938 MeV ≈ 1%), we found that the different methods of accounting for the binding energy produce small but noticeable differences in pion multiplicities and mean transverse momentum at low collision energies of The momentum distribution of nucleons in the ground state nucleus is a uniformly filled sphere in momentum space, known as the Fermi sphere. The radius of the Fermi sphere is given by the formula where ρ( r) is the density of nucleons at the point r. A more detailed description of the density calculation is given in Section III D. A typical value of p F ≈ 300 MeV corresponds to an energy of p 2 F /(2m N ) ≈ 45 MeV. To make sure that Fermi momenta are generated correctly, in Fig. 6 we plot the momentum distribution of the neutrons in a lead nucleus from SMASH and compare it to theoretical expectation computed as follows. In the central part of the nucleus, where the density is uniform, the expected normalized distribution , while the analogous momentum distribution integrated over the whole nucleus with ρ(r) where x = p/p 0 F and α = d/R. Including the Fermi motion is only sensible if potentials are turned on simultaneously. Otherwise, the nucleus will fly apart due to the finite transverse momenta of the nucleons that need to be compensated by the attractive mean field interaction. Alternatively, one may 82 Pb nucleus neutrons r < 3.5 fm, SMASH full nucleus, SMASH r < 3.5 fm, expected full nucleus, expected N 1 d(p/pF) 3  employ the so-called frozen Fermi approximation: Fermi momenta are used for collisions, but not for propagation. This option is currently not implemented, but will be considered in the future. Fig. 7 shows the nuclear stability over a large time range, much larger than what is actually relevant for a nucleus-nucleus collision. The nucleons fly apart as expected, if only Fermi motion without potentials to stabilize the nucleus are included. With potentials there is the expected oscillatory behavior: The nucleons drift apart due to Fermi motion and the potentials counteract and push them closer together again.
Computations with potentials require that time step is small enough -the energy change per timestep should be much smaller than the energy of the particle: As an estimate for the maximal |∂U/∂r| max let us take two nucleons at the same point and consider |∂U/∂r| max = 2m N /σ, where σ is the width of the Gaussian smearing (as defined in Eq. (68)) ∆t σ/2 .
Assuming the default value of σ = 1 fm, a time-step size of ∆t = 0.1 fm is reasonable for physically relevant cases. Since the potential becomes smoother with higher number of test particles N test , the estimate becomes in this case c. Deformed Nuclei Despite the rather symmetric nuclei that are most often used for heavy-ion collisions, sometimes it is of interest to study deformed nuclei as well. For example, uranium has a prolate shape according to its nuclear many-body wave function. At RHIC U+U collisions at √ s NN = 193 GeV have been studied to evaluate the multiplicities and anisotropic flow as a  function of geometry. Especially the case of tip-tip collisions, where the multiplicity is very high but the elliptic flow is close to zero, and the case of body-body collisions, where the elliptic flow is maximal, are of great interest. To differentiate between the different geometries, Monte Carlo event generators that yield the correct trends for the observables are needed that take into account the more involved geometry of deformed nuclei [36].
In SMASH, the Woods-Saxon distribution is enhanced with an angular dependent radius r(θ, ϕ) The deformation dependent nuclear radius r(θ, ϕ) can be described using the β parameterization [37] r(θ, ϕ) Here r 0 is the initial nuclear radius and the coefficients in front of the spherical harmonics Y m l are called the β shape parameters (or deformation parameters). Note that our deformed nuclei are azimuthally symmetric and hence all terms with nonzero magnetic quantum number will vanish.
In SMASH, the deformed Woods-Saxon has been implemented by a rejection sampling routine. For the deformation we use the β shape parameters up to angular momentum quantum number l = 4 from [37]. For the initial nuclear radius r 0 , we use values from [38] (see Two-Parameter Fermi Model, abbreviated 2pF). We also have a default initial radius that uses the empirical relation 1.2A 1/3 . The diffusiveness parameter d is on average given by 0.54 fm [39]. Adjustments for specific nuclei come from [38]. We sample a polar angle from the uniform solid angle distribution, and for the radius we set our maximum sampled value to be r max = r 0 /d + r 0 d. Finally, the saturation density ρ 0 represents our normalization condition: A deformed nucleus is no longer invariant to rotation. We therefore need to rotate the nucleus during the initialization phase. To do so, we treat the system of nucleons like a rigid body. A set of Euler angles is uniformly sampled that determines the rotation of the specific nucleus, for which we use the notation convention (ϕ, θ, ψ).
To visualize the differences between collisions of symmetric and deformed nuclei in Fig. 8, the ratio between total transverse energy of charged particles and the beam energy scaled with the number of nucleons (which is used to determine the centrality classes in low energy collisions), is compared for Au+Au and U+U collisions at √ s NN = 3 GeV. It can be seen that there are less high multiplicity events and more intermediate multiplicity events in U+U than Au+Au collisions. The difference comes from the fact that many tip-body and nonoverlapped body-body collisions (with impact parameter b = 0 fm) do not produce as many new particles as in Au+Au most central collisions. Those events are selected as semi-central collisions at experiments but provide much smaller elliptic flow. For most peripheral collisions, there are more non-empty events in U+U than Au+Au collisions, since the uranium can touch each other with much larger impact parameter along their long axes. d. Frame Invariance To define the kinematics of the heavy-ion collision, different variables are commonly used. At lower beam energies often the kinetic energy per nucleon E kin or the momentum per nucleon p lab is given. Up to moderate beam energies of around 160 AGeV per nucleon, most experiments are fixed target experiments to increase the luminosity. Only if the whole available energy out of the accelerated bunches is needed to reach higher energies for the collision, the center-of-velocity frame of the two nuclei equals the laboratory frame in a collider setup. In this case, the center of mass energy for binary nucleon-nucleon collisions is usually specified √ s NN to characterize the collision energy. This frame is the standard computational frame for SMASH calculations, which is equal to the center-of-mass frame for symmetric systems.
To give an estimate on how much the Lorentz invariance is violated by the non-local collision criterion, the number of interactions in one physically identical heavyion reaction is counted for calculations in different reference frames (see Fig. 9). The calculations have been performed as a function of beam energy. All interactions above a cut-off of √ s = 2 GeV per binary collision are counted until the particles freeze out, therefore no Lorentz transformation of the runtime is necessary. The cut-off is necessary to exclude collisions within the nuclei at low momenta that are not relevant for the actual heavy-ion collision. These calculations do not assume specific time-steps, but all particles are propagated to the next interaction. Apart from the general trend that there are more collisions at higher beam energies, the relative difference between the reference frames is very small at all beam energies. The two calculations coincide within the statistical error bars at all energies.

Infinite matter calculations
To simulate infinite hadronic matter or other simple systems like an ideal massless or massive gas and investi-gate its thermodynamic properties, box calculations are performed. This section describes the initialization of N particles of species i in such a box. In general, every particle j is characterized by coordinates (x j , y j , z j ), the four-momentum (E j , p j ) and a spectral function, therefore the particle mass is given as m = E 2 − p 2 and is not necessarily equal to its pole mass. The coordinates of the N particles (x j , y j , z j ) are sampled uniformly in the box: where U denotes the uniform distribution and L is the length of the box. The momenta of the particles are sampled using the thermal Boltzmann distribution with temperature T : (24) where w(p) is a probability to generate momentum p, θ and ϕ are angles in spherical coordinates and N is a normalization factor. In other words, momentum directions are sampled uniformly in the solid angle dΩ = sin θdθdϕ. Let us denote the total momentum of N particles sampled from this distribution p tot . One can see that the ensemble average of p tot is zero, because it involves an integral over an odd function. However, in each single event p tot = 0, which is corrected by changing the momentum of every particle p j → p j − p tot /N . After this procedure the thermal distribution is slightly spoiled, the total energy is changed and angle uniformity is disturbed. This is a small effect for large numbers of particles N 1. After letting the system thermalize, the temperature differs by 1-2% from the initialization temperature. One also has to note that the total energy is not the same from event to event, it is fluctuating, even without this momentum shift. Fixed are the volume V , the number of particles N and the temperature T . This picture corresponds to the canonical ensemble (CE) with the temperature and the particle number as independent parameters.
After initialization particles propagate along straight lines with velocities v i = p i /E i and collide with each other. The simulation is time-step-based and uses a grid to increase the performance of the collision finder as described above in Section II A 5. The box has per default periodic boundary conditions: at the end of each time step, particles outside of the box with coordinates r are returned to the coordinate r mod (L, L, L).

Expanding sphere
A simplified scenario including expansion can be initialized using a three-dimensional sphere. For this purpose, N particles of different species i are uniformly distributed in a sphere with radius R. The momenta are sampled from a thermal distribution analogously to the box initialization. Then, the system expands freely. This setup provides the opportunity to analyze the numerical stability of the code by comparing to analytic solutions like [40], which is left for future work.

Afterburner for hydrodynamic simulations
The fourth method for initializing SMASH is to provide an externally generated particle list based on which the calculation is started. In heavy-ion collisions at high beam energies the hydrodynamic hot and dense stage is followed by a dilute phase that is dominated by hadronic rescattering and resonance decays. To include this late stage dynamically, a hadron transport approach like SMASH needs to be run for each particle configuration that is provided by sampling on the Cooper-Frye hypersurface. The hadronic transport calculation can be coupled in a similar way to other approaches than hydrodynamics, if necessary. Note that input particles in the list are not required to be at the same time t in the computational frame, successive appearance of particles are implemented by setting non-zero formation times. These particles are propagated back to the earliest time in the list and free stream before their formation time.
C. Particle Properties

Particle species
We implement the most well-established hadronic states from the Review of Particle Properties [41] with their corresponding decays and cross sections as detailed below. These particles and their properties are summarized in Table II.
To simplify the extrapolation of cross sections and particle properties, full isospin symmetry is assumed. Therefore, small differences in the masses between isospin partners have been neglected. However, it should be noted that the cross sections for certain processes can indeed depend on isospin (thus breaking isospin symmetry, e.g. in channels like N N → N N * , see Section II D).
We treat all particles as stable which have a width below 10 keV (such as the π, η, K, N, Λ, Σ, Ξ, Ω). All unstable particles ("resonances") are assumed to have a Breit-Wigner shape. We note that this approximation is known to be questionable for the σ meson, our parameters are adjusted to reproduce the ππ elastic cross section.

Spectral functions
In general, the spectral function encodes the dispersion relation for a particle and can depend on the temperature and the density of the system. Medium modifications are currently neglected in SMASH and all spectral functions are described by relativistic Breit-Wigner distributions: Here m is the actual off-shell mass of the resonance and M 0 is the pole mass (i.e. a constant given in Table II). However, the total width Γ is not constant, but given by the mass-dependent width function Γ(m). Each resonance has a minimum mass m min (corresponding to the sum of masses of the lightest decay channels), below which the width, and thus also the spectral function, vanishes. The total width is computed as the sum of all partial widths: Note that the width given in Table II is the total onshell width, i.e. Γ 0 = Γ(M 0 ). The spectral function in relativistic Breit-Wigner form is normalized to one, when integrated from zero to infinity: In practice the integration can start from m min , since the spectral function vanishes below that value. Under the assumption of a constant width Γ, the normalization factor is exactly N = 1. As soon as the width becomes mass-dependent (as it is the case in SMASH), the normalization factor N can deviate from one and needs to be determined numerically. Practically all the normalization constants in SMASH are still rather close to one (within 25%).

Decay widths
All the decay widths in SMASH are currently calculated following the treatment of Manley et al. [42], where in general the width of a two-body decay R → ab is written as Here m is the actual off-shell mass of the resonance R, M 0 is its pole mass, Γ 0 R→ab = Γ R→ab (M 0 ) is the partial width at the pole mass and the function ρ ab is defined as In this formula, m a and m b denote the (off-shell) masses of the particles a and b (which are being integrated over), A a and A b are their spectral functions and | p f | is the absolute value of the final-state momentum of a and b in the center-of-mass frame, which is given by: Finally, L is the orbital angular momentum of a and b in the final state and B L are the so-called 'Blatt-Weisskopf functions' [43]. The parameter R is usually called the 'interaction radius' and is assumed to have a universal value of R = 1 fm for all processes. The form factor F ab is only relevant for unstable decay products and will be discussed later.
The simplest case is that of a resonance R decaying into two stable daughter particles. Popular examples are ∆ → πN or ρ → ππ. In this case, the daughters have fixed masses (i.e. their spectral functions are just δ functions), so that the integrals collapse: As an example, the width for the p-wave (L=1) decays of the ρ and ∆ (mentioned above) becomes . Here m and M 0 are the offshell and pole mass, respectively, while p f and p f,0 denote the final-state momenta in the center-of-mass frame for mass m and M 0 , respectively. Λ = 1/R can be viewed as a cut-off parameter. For an s-wave (L=0) decay like σ → ππ, the width simply becomes since B 2 0 = 1. In the case that one of the daughter particles is itself a resonance, the width calculation becomes more difficult, since the mass of this daughter resonance is not fixed and needs to be integrated over. Examples for this case are N * (1440) → π∆ or ω → πρ. As one of the daughters is stable, at least one of the two integrals collapses: The remaining integral runs from the minimum allowed mass of particle a (i.e. the threshold of its lightest decay channel) up to the maximum possible mass of a in the decay process (given by m − m b ). The form factor F ab (by M. Post [44]) is used only if unstable decay products are involved and is defined as where the cut-off factors given in Table III are used.
It is easy to see that F ab (M 0 ) = F ab ( √ s 0 ) = 1. Note that this form factor was not used by Manley originally, but was added only later in the GiBUU implementation. The effect of the form factor is that it suppresses the highmass tail (m > M 0 ) and slightly enhances the low-mass tail (m < M 0 ). Both of these effects get stronger with decreasing λ (F ab → 1 for λ → ∞). We have decided to follow the GiBUU framework for the width parametrization of resonances, since it has been proven to give a good description of experimental data [26]. All the formulas described above are for the case of resonance decays. For the inverse process, i.e. resonance formation via ab → R, the Breit-Wigner cross section involves the so-called 'in-width' Γ ab→R . For stable particles a and b it is identical to the 'out-width' Γ R→ab . However, the two differ if a or b are unstable. In the Manley formalism, the in-width for unstable particles becomes where m is the off-shell mass of the produced resonance R (i.e. the √ s in the process) and p ab = p cm (m, m a , m b ) is the momentum of a and b in the center-of-mass frame. The difference between the in-and the out-width is essentially due to the fact that for the out-width one integrates over the mass of the unstable particle, while in the in-width this mass is fixed.
In Fig. 10 the theoretical decay width of the N * (1440) resonance is shown as a function of mass. The total width is given as the sum of all partial widths. Each partial width has a threshold that is given by the sum of the minimal masses of the decay products. The branching ratios are fixed at the pole mass. One can see that all partial widths increase as a function of mass, since more phase space is available for heavier resonances. The lifetime correspondingly has an opposite trend and heavy particles decay faster than low-mass resonances. Since the width also enters in the production cross section (Eq. (39)), the production of such low-mass resonances becomes more unlikely. total π + n π 0 p π + ∆ 0 π 0 ∆ + π − ∆ + + σp

D. Collision Term
The collision term includes all different processes (decays and collisions) that can happen to particles within this hadronic transport approach. At this point, unstable particles can decay, 2 particles scatter in-/elastically or excite a resonance. Weak decays are neglected since they have significantly longer lifetimes than processes associated with the strong interaction. Electromagnetic processes are treated perturbatively and will be discussed in detail in a forthcoming publication [45].
We note that the current implementation is limited to the energy regime of a few GeV, where all hadronic cross sections are expected to be dominated by the excitation and decay of resonances. Since the model at present is lacking a string fragmentation mechanism, the cross sections are not sufficient at higher energies. In the following, a detailed description of all implemented processes is given.

Decays
The lifetime of a resonance is defined as τ = 1/Γ(m), where Γ(m) is the mass-dependent total decay width. The probability to decay in a sufficiently small time interval ∆t is This leads to exponential decay, as the survival probability after n time steps is P (alive after n steps) = (1 − Γ(m)∆t) n when ∆t → 0. As noted above, the total width Γ(m) is computed as the sum of all partial widths. When a resonance decays in SMASH, the decay channel is randomly chosen from the list of allowed channels for this particle, based on the off-shell branching ratios Γ i (m)/Γ(m). The decay channels and their on-shell ratios Γ i (M 0 )/Γ(M 0 ) are listed in an input file and can be turned on and off separately.

2 → 1 processes
The cross section formula for 2 → 1 resonance production is based on Eq. (176) in [26]: where: • J is the spin of the particle • Γ ab→R is the partial in-width for the process • A R is the spectral function of the resonance • S ab is a symmetry factor, which is 2 if a and b are identical, and 1 otherwise is the center-of-mass momentum of the initial state Note that in the above, Γ ab→R (s) refers to the isospinspecific channel instead of the isospin-generic channel. Hence there is no need for isospin factors in the cross section formula.
The so-called "in-width" Γ ab→R simply equals the usual decay width Γ R→ab for the case of stable particles a and b, see Section II C 3. For unstable particles however, it is given by Eq. (37), which differs from the decay width.

Elastic collisions
There are different cases of elastic collisions in SMASH. For the meson-baryon and meson-meson collisions, one assumes that the elastic cross sections are fully determined by resonance excitation and decay, e.g. πN → ∆ → πN or ππ → ρ → ππ. For baryon-baryon collisions on the other hand, one typically uses parametrized cross sections. The parametrizations of the elastic pp and pn cross sections in particular are taken from [46], eq. (44) and (45).

2 → 2 processes with one resonance in final state
When there is another particle in the final state, the resonance mass must be integrated over the allowed range: where p i and p f are the center-of-mass momenta of the initial and the final state and A R is the spectral function of the resonance R. The symbol C refers to isospin Clebsch-Gordan factors, which couple the initial and final state to a total isospin I. Here it is assumed that the matrix element |M| 2 is a constant (or only depends on s) without angular dependence, resulting in the factor 4π from the trivial angle integration. If the matrix element does depend on the angle, the factor 4π must be replaced with the proper integration of |M| 2 over the phase space. The lower mass limit for the resonance, m min R , is defined as the sum of the particle masses in the lightest decay channel. This is the lowest mass the resonance can have and still be able to decay into one of the implemented channels.
For the process N N → N ∆, the parametrized energy dependence (with the parameters A = 68, b = 1.104 GeV and c = 1.951) is based on a fit to the Dmitriev one-bosonexchange (OBE) model [47]. For other resonance production processes (i.e. N N → N R and N N → ∆R, with R = N * , ∆ * ), the matrix element is assumed to be a constant (independent of s), but can depend on the total isospin and the pole masses m a and m b of the outgoing particles. It is parametrized as with parameters A I as given in Table IV.

2 → 2 processes with two resonances in final state
Analogously to Eq. (40), one can write down the cross section for a process with two resonances in the final state. In this case both their masses must be integrated over the allowed range: The double-resonance production processes that are currently implemented in SMASH are N N → ∆∆, ∆N * and ∆∆ * . The matrix elements are parametrized in the same way as for single-resonance production, see Eq. (42) and Table IV.

Detailed balance
The cross sections for the inverse resonance-absorption processes are derived from the production cross section by imposing the principle of detailed balance (see Eqs. (B.6), (B.9) and (181) in [26]): This equation holds for both single-and doubleresonance absorption, i.e. c and d can be either two resonances or a resonance and a stable particle. The symmetry factors S xy here are defined such that they are 2 if x and y are in the same isospin multiplet and 1 otherwise. In SMASH all processes are following explicit detailed balance in the whole phase space, as will be demonstrated in Section III C below.

Mass sampling
In any process where a resonance is produced in the final state, its mass needs to be sampled according to the spectral function and the available phase space. The simplest case is that a single resonance is produced in a 2 → 2 collision together with a stable particle. Then the mass of the resonance is sampled from the integrand of Eq. (40): The allowed mass range is from m min R to √ s − m stable , where s is the Mandelstam s of the process and m stable is the mass of the stable final-state particle.
The mass sampling is slightly more complicated for the case of a resonance decay (1 → 2) with one resonance and one stable particle in the final state. In this case an additional Blatt-Weisskopf factor appears, which takes into account the angular momentum in the decay, cf. Eq. (35): For a scattering process with two resonances in the final state, the masses of both resonances have to be chosen according to the function which is the integrand of Eq. (43). It is important to note that both masses cannot be determined independently, but have to be chosen simultaneously according to a common sampling function.
Analogously to the single-resonance case, a decay into two resonances also includes an additional Blatt-Weisskopf factor: Drawing random numbers from these distribution functions is numerically non-trivial. We first draw from a Cauchy distribution which approximates the spectral function and handle the remaining factors by rejection sampling (where the unknown maximum value is determined adaptively).

Angular distributions
We currently have anisotropic angular distributions implemented for N N → N N , N N → N ∆ and N N → N R (with R = N * , ∆ * ). For elastic nucleon-nucleon collisions we follow the prescription by Cugnon et al. [48], using an exponential ansatz dσ/dt ∝ e −bt , with an energy-dependent parameter b which is fit to data. In the second case we also follow Cugnon et al. [48], using the same ansatz as for elastic NN collisions. For the last case of N N → N R we use the ansatz dσ/dt ∝ t −a , with parameters a which have been fitted to HADES data [27]. In Section III B a comparison to elementary data is shown. We note that in the present implementation all resonances decay isotropically in SMASH.

Pauli blocking
Pauli blocking is an effective way to obtain the solution of the quantum BUU (Boltzmann-Uehling-Uhlenbeck) equation from classical molecular dynamics. To understand the way this is achieved one has to compare the classical Boltzmann equation and the BUU equation -its quantum analog: Here the plus sign is for bosons and the minus sign for fermions. One can see that quantum BUU equation differs from classical Boltzmann only in the Uehling-Uhlenbeck factors in the collision term. One can interpret this factors as a multiplication of the cross sections by i (1 ± f i ), where the product is taken over all final states in the reaction and f i ≡ f (r i , p i , t) is the phasespace density of final-state particle i. This means that for bosons cross sections are effectively increased and for fermions cross sections are effectively decreased. This is called Bose enhancement and Pauli-blocking respectively. While Bose enhancement has been attempted to implement recently in a parton cascade [49], Pauli blocking is taken into account in many transport approaches. Since Pauli blocking is important in the energy range under consideration in this work, we describe in the following how it is taken into account in a Monte-Carlo model. The implementation of Pauli blocking consists of two parts: the calculation of the phase-space density and the rejection of reactions with probability 1 − i (1 − f i ). For the latter SMASH loops over all baryons in the final state after a collision has taken place and returns 'true' for blocking, if a uniformly distributed random number r > f i . This means that the reaction is not blocked with probability i (1 − f i ). In this way, no fermion can be produced or scatter into a phase space bin that is already occupied by another fermion.
The implementation of the phase space density calculation basically follows the method used in the GiBUU model, see section D.4.3 in [26]. By definition N (∆V r , ∆V p ) = gf (r, p)∆V r ∆V p , where N is the number of (test)particles in a given phase-space volume ∆V r ∆V p and g is the degeneracy. Theoretically, the size of the phase-space goes to zero ∆V r , ∆V p → 0. In practice ∆V r , ∆V p and the way of averaging are chosen to balance between the smoothness of the obtained distribution function and the resolution of coordinate and momentum space. This implementation relies on a large number of test particles (N test 20).
The phase-space density is calculated according to the following equations: with κ given as Here r j is a vector connecting the point, where f is calculated, and the position of the j-th particle. All these expressions can be analytically further evaluated for r c > r r . This is a reasonable assumption, because the Gaussian cut-off r c has to be large enough, so that the results do not depend on it. If r c < r r the whole method is hardly applicable. In GiBUU these integrals are computed numerically, but we have found analytical expressions for them (see Appendix A). For V p a sphere of radius 80 MeV is taken.
In Fig. 11 the number of collisions that is blocked due to prior phase space occupation has been calculated in central Cu+Cu and Au+Au collisions as a function of beam energy. One can see that at very low energies there are as many blocked collisions as collisions taking place. The ratio drops rather fast and around E kin = 2A GeV only a quarter of the collisions are blocked. It then saturates around 10% for higher beam energies. Fig. 12 demonstrates the need for a decent number of test particles to obtain stable results. If the number of test particles is low the phase-space volume cannot be calculated with enough precision and therefore, there are too many collisions allowed. Saturation sets in around N test = 20 and is very similar for Au+Au and Cu+Cu collisions. In particular the π − p cross section in the upper panel of Fig. 13 shows some very clear resonance structures. The lowest excitation here is the ∆(1232), followed by several N * resonances in the second and third resonance region at around 1.5 and 1.7 GeV, respectively. ∆ * states only play a significant role at higher energies of around 1.9 GeV. In fact SMASH exclusively produces s-channel resonances in this case, which then decay into different final states. In this way, we can saturate the total cross section up to about 2 GeV with only minor deviations, which may be caused by the negligence of non-resonant backgrounds and/or uncertainties regarding resonance parameters.   [41]. Also the elastic cross section only involves contributions from s-channel resonances, which then decay back into π − p, and is reasonably well described over most of the displayed energy range. Only above energies of 2 GeV, SMASH starts to underestimate the total and elastic cross section. Here further production mechanisms, such as string fragmentation, will be necessary to achieve agreement with the data.
The π + p cross section in the bottom panel of Fig. 13, shows a similar dominance of s-channel resonances. However it is limited to ∆-type excitations due to isospin arguments. The resonance contributions in π − p and π + p are related by simple Clebsch-Gordan factors.
The purely mesonic case of the π + π − cross section in Fig. 14 exhibits a similar resonance pattern. Here the dominant resonances are the ρ and f 2 states. There is also a contribution from the scalar σ (or f 0 ) meson. However, it should be noted that the parameters (mass and width) of the σ in SMASH differ significantly from the PDG values [41], in order to achieve a reasonable agreement with the ππ data. Presumably this discrepancy is due to our usage of the Breit-Wigner approximation, which is known to be questionable for a state like the σ meson, for which the width is comparable to the mass.
For the nucleon-nucleon cross sections in Fig. 15, the resonance contributions are less apparent, simply because the resonances do not occur in the s-channel. Instead the prevalent physical picture in this case is a t-channel meson exchange, which may excite one or both of the scattered nucleons into a resonance state that subsequently decays. Both the pp and pn cross sections include a significant elastic contribution that rises towards the threshold. We simply parametrize the √ s dependence in this case, cf. Section II D. The first inelastic channel that opens up is the excitation of a single ∆ resonance. At higher energies it is followed by the excitation of heavier resonance states (N * and ∆ * ) as well as double-resonance excitations. For the nucleon-nucleon case, the resonancebased mechanisms are able to saturate the total cross section up to energies of 4 to 4.5 GeV, above which they need to be supplemented by additional production mechanisms (e.g. string fragmentation). Further it should be noted that the √ s dependence of the total cross section is not described perfectly well here, which may be caused by assuming matrix elements which are independent of s, e.g. in Eq. (40).
The exclusive cross section for single pion production in proton-proton collisions in Fig. 16 shows an overall good agreement with the data. The dominant contribution for single pion production in nucleon-nucleon collisions is the ∆ resonance (compare Fig. 15). Above 2.5 GeV also additional contributions from excited resonance states (N * and ∆ * ) occur. A slight undershoot for the π 0 production at low energies in proton-proton collisions might come from non-resonant background terms that are not included in the model. Fig. 16 also reveals a systematic undershooting for the single pion production in proton-neutron collisions, which could be due to an underestimation of the contributions with total isospin I = 0, see Section II D 4.

B. Angular distributions
In Fig. 17 we show two examples of angular distributions dσ/dt in pp collisions, t being the Mandelstam variable. The upper plot shows a collision at a relatively low energy, where essentially only the elastic and single-∆production channels are open. The angular distribution of the elastic channel is of course symmetric in the allowed t-range and matches the data points rather well, even though the slope at this particular energy appears to be slightly too flat. The distribution for single-∆ production is not symmetric and restricted to a smaller range in t, due to the larger mass of the ∆ in the final state. Unfortunately there is no inelastic data to compare to at this energy.
The lower plot in Fig. 17 shows a pp collision at a somewhat higher energy, where additional resonance production channels are open. In principle the distributions for all these channels are forward/backward-peaked (either exponential or power-law shaped), as mentioned in Section II D 8. This forward/backward peaking is clearly visible for the N N and N ∆ final states at least, while those final states with heavier resonances exhibit a more plateau-like structure, due to the limited phase space and the mass distributions of the resonances. Here the sum of all inelastic channels is compared to data and indeed shows a reasonable agreement, again with a slight tendency of being too flat.

C. Detailed balance
The strong interaction is invariant under time reversal, which implies that for any scattering or decay process the probability of transition w(Γ i , Γ f ) from the point in phase space dΓ i to dΓ f is equal to the probability of the reverse process.
Eq. (53) is embodied in SMASH via the equality of matrix elements of the forward and backward reactions, With this formula one can connect cross sections of the forward and backward 2 → 2 reaction, or the width of the 1 → 2 decay to the backward 2 → 1 reaction. For count reactions example, for 12 → 1 2 scatterings where I = (P 1 · P 2 ) 2 − m 2 1 m 2 2 and the term 1/(1+δ 1 2 ) accounts for identical particles in the final state. Integrating this over momenta one arrives at Eq. (40). For resonances in the final state the transformation from [55] is applied. The corresponding Eq. (39) for decays is derived analogously.
Substituting Eq. (53) back into the Boltzmann equa-tion 1 leads to the principle of detailed balance: In equilibrium the rate of forward reactions dΓ i → dΓ f is equal to the rate of backward reactions [56].
To test, if detailed balance actually holds in our calculations, a periodic box is initialized with multiple particle species. After the matter reaches equilibrium, we check that the numbers of forward and backward reactions are identical. The fact that the box should reach equilibrium is granted by the H-theorem, which is derived assuming Eq. (53) and the hypothesis of molecular chaos (twoparticle distribution function f 2 (Γ 1 , Γ 2 ) = f (Γ 1 )f (Γ 2 ) or, in other words, participants of the reaction are uncorrelated). Strictly speaking, in a transport code both assumptions are valid only in the limit N test → ∞. At finite N test the interactions are non-local due to the geometrical cross sections. In addition, while two particles with space coordinates r 1 and r 2 form a resonance at ( r 1 + r 2 )/2, the products of resonance decay gain the same position as the decaying resonance. This breaks Eq. (53), where for nonlocal interactions the phase space Γ includes coordinate space. This leads to a small violation of detailed balance, which vanishes at large N test as we show in the following.
For the test we are using two configurations: a ρ−π−σ box and a N − π − ∆ box. The first one is initialized with a 100 π + , 100 π − and 100 π 0 in a volume of V = (10 fm) 3 . The reactions ππ ↔ ρ and ππ ↔ σ are allowed, while all the other possible reactions are switched off. From Fig. 18 one observes that the system reaches chemical equilibrium, since the particle multiplicities in the box saturate after around t = 20 fm/c. Starting from this time, forward and backward reactions are counted. The matrix elements of reactions in the same isospin group differ only by Clebsch-Gordan coefficients. Thus one expects, for example, that the number of reactions N (σ ↔ π + π − ) = 2N (σ ↔ π 0 π 0 ). Therefore, the reaction numbers in Fig. 18 are scaled by the isospin and symmetry factors appropriately to make sure that this expectation is fulfilled. Detailed balance is valid not only for the total number of reactions, but it also has to be fulfilled differentially in momentum space. We show in Fig. 18 that detailed balance is indeed fulfilled differentially in each invariant mass bin of the reaction. Let us note that for the ρ − π − σ box detailed balance for the total (but not differential) number of reactions follows trivially from the multiplicity saturation. Indeed, denoting forward and backward reaction rates by r → and r ← , one arrives at For the N − π − ∆ box similar relations become less trivial. We initialize the N −π −∆ box with 100 neutrons and 100 protons and allow reactions ∆ ↔ N π (1), N N ↔ N ∆ (2) and N N ↔ ∆∆ (3), with all the other reactions being forbidden. In chemical equilibrium the following equations are fulfilled: It can be observed that forward and backward rates for N N ↔ N ∆ and N N ↔ ∆∆ being equal does not necessarily follow from multiplicities being saturated. As one can see from Fig. 20, with N test = 100 detailed balance is violated at maximum by 2%. For N test = 1 this violation can reach 10% because of the non-locality effect described above. To see if the numbers of reactions within one isospin group relate as expected from Clebsch-Gordan factors, we multiply every number of reactions N i by a factor α i that compensates for the isospin factors of this reaction. Let us denote N isospin group = 1 k k i=1 α i N i , where k is amount of reactions in the isospin group (forward + backward). If the SMASH result corresponds to the theoretical expectation, then N i / N isospin group should be strictly 1 for every reaction. One can make sure from Fig. 18 and from Fig. 20 that SMASH matches this expectation. Table V shows the origin of compensating coefficients α i . While most of the Clebsch-Gordan factors are simple, for pn ↔ ∆∆ reactions they are less intuitive. The matrix element for N N ↔ ∆∆ reaction is isospin dependent, namely |M (I = 0)| 2 = κ|M (I = 1)| 2 , where κ = 8 3 . Here is one explicit example illustrating the calculation (where states beyond I = 1 have been omitted, since they drop out): pn|∆ − ∆ ++ 2 = 5κ + 9 40 Thus, we have shown that the detailed balance in SMASH for a mesonic system and a more complex situation involving baryons and mesons is fulfilled.

D. Thermodynamics
To investigate the thermodynamic properties of the hadron gas, the energy-momentum tensor T µν ( r) and four-currents j µ ( r) can be calculated from the particle distribution functions. These two quantities provide access to the energy density and particle number density in the corresponding rest frames. Assuming that the potential energies of particles are small compared to their Reaction Clebsch Symmetry Total   kinetic energies and taking into account that collisions happen instantaneously, the corresponding equations for non-interacting particles are applied: where f ( r, p) is the single-particle distribution function. For a discrete set of particles it reads For numerical calculations we substitute the deltafunction by the smearing kernel where ∆ r = r − r part , β = p part /E part is the 3-velocity of the particle and γ = (1 − β 2 ) −1/2 . It is shown in [57] that this kernel has proper Lorentz-transformation properties, is normalized to 1 and represents a simple 3D-Gaussian in the rest frame of the particle. The equations for the numerical evaluation of thermodynamic quantities are then where N ev is the number of events and N test is the test particle number. In the limit of the smearing width σ → 0 and N ev N test → ∞ the full smooth quantities are obtained. This limit is numerically challenging, because when reducing the smearing width σ, one has to increase statistics, keeping σ 3 N ev N test = const. Therefore, we take reasonably small σ = 1 fm and keep in mind the smearing effect, which is demonstrated in Fig. 21 for the density calculation of a Pb nucleus comparing σ = 0.5 fm and 1 fm. The Eckart rest frame density is obtained as ρ Eck = j µ j µ . For net baryon (charge, isospin projection) density a naive weighting of particles in Eq. (70) with their baryon numbers can give rise to j µ j µ < 0. Therefore, we compute ρ = ρ + − ρ − , where + corresponds to positive baryon number (charge, isospin projection) and − corresponds to negative ones. In Fig. 22 the dependence of the net baryon density versus time in the middle of the target in the central Au+Au collision at E kin = 0.8A GeV in the fixed-target frame is shown. The energy density in the Landau frame is depicted in Fig. 23. Both figures show that the ground state baryon/energy density values are reproduced, when the collision term is disabled. Including interactions the baryon/energy density rises to about 4 times the nuclear ground state densities.
In many applications (e.g., connecting non-equilibrium initial states to relativistic hydrodynamics) the Landau rest frame (LRF) quantities are needed. By definition, T 0i LRF = 0, the energy flow in the LRF is zero. To find the LRF we solve the generalized eigenvalue problem (T µν − λg µν )h ν = 0, where g µν is the metric tensor. The eigenvector corresponding to the largest eigenvalue is proportional to the 4-velocity of the LRF and the proportionality constant is fixed by the constraint that u µ are plotted in the x-z-plane in Fig. 24 for a Au+Au collision. One can observe the onset of radial flow after the initial collision of the two nuclei. We note that the LRF energy density before collision reproduces again the nuclear ground state energy density.

IV. RESULTS FOR HEAVY-ION COLLISIONS
In this section we compare particle yields and spectra in heavy-ion collisions calculated with SMASH to experimental data from the HADES and FOPI collaborations. The focus for the current analysis lies on pions, because they contribute the majority of the newly produced particles; and on protons, because they are part of the initial system before the collision.
Some time after the collision, the particles don't interact anymore and thus their momenta are frozen. Therefore, the basic bulk observables to quantify the dynamics of the collision are rapidity and transverse momentum spectra. To obtain Lorentz-invariant spectra, the longitudinal rapidity y and the transverse mass m T are used as momentum coordinates: Usually the rapidity y is rescaled to y 0 such that the nuclei are located at y 0 = ±1 before the collision: where y cm is the rapidity in the center-of-mass frame.
To obtain sensible comparisons between our calculation and experimental data the procedure to select centrality classes needs to be the same. See Appendix C for how this is done for the FOPI data.
First, let us have a look at the total pion multiplicities and their averaged transverse momentum over a broad range of energies. In Fig. 25 the total multiplicities of charged pions in central Au+Au collisions at kinetic energies from 0.4A GeV to 1.6A GeV are compared to FOPI measurements [58]. The upper plot shows the total pion multiplicity, the one in the middle shows the ratio of negative pions to positive pions to indicate the isospin asymmetry. The lower plot shows the average transverse momentum of the pions. The impact parameter b for the SMASH events was sampled from a minimum bias distribution with b < 2 fm onto which the corresponding ERAT cuts have been applied. The simulations were run successively with and without potentials (see Eq. (10)), Fermi motion (see Section II B 1) and Pauli blocking (see Section II D 9). Potentials and Pauli blocking require a sufficient number of test particles to function properly. When any of these features was enabled, 20 test particles were used instead of one.
Without potentials (and Fermi motion and Pauli blocking) the SMASH results agree well with the data, except for the lowest energy at 0.4A GeV. A deviation at low energies is expected, because potentials should have a strong effect there. Running the cascade with 20 instead of one test particle per real particle, there is a slight increase in multiplicity. This effect should be considered a systematic error of the model, since changing the test particle number is not supposed to affect the physics. Additionally enabling the potentials (which are soft, see Section II B 1) decreases the pion multiplicities by a large amount. Adding Fermi motion to the simulation yields the strongest effect and increases the multiplicities. Pauli blocking causes a small decrease in multiplicity. For the more physical scenario with all features enabled there is an overestimation of the number of pions at all energies. Such an overestimation and a decrease of the multiplicities due to soft potentials and Pauli blocking has been observed with one of the first transport models as well [59].
The pion ratios look similar with and without potentials. Only for the lowest energy the results with potentials are a bit closer to the experimental values. Please note that no Coulomb potentials are included in this calculation. In an earlier comparison with the FOPI data for gold-gold collisions at 1.5A GeV, it has been suggested that Coulomb potentials "almost exclusively" account for the difference in the momentum spectra of the charged pion species [58]. The results here do not support this claim, because the total relative multiplicities of the pions are reproduced without any Coulomb potentials.
The transverse momenta do not vary significantly among the different pion species or with and without potentials. It is difficult to pin down the reason for the overestimation of the pion multiplicity. At this energy a lot of implementation details can influence the multiplicity significantly: Fermi momenta, potentials and the N ∆ cross sections (which haven't been measured) introduce some uncertainties. The cross sections can be reduced by in-medium effects [60,61], which is unaccounted for in SMASH. These in-medium effects would reduce the number of produced pions. More work is needed to understand the exact reasons for the discrepancy. On the other hand, SMASH is primarily designed for FAIR energies, where potentials will be less important, and the results are similar to those of other approaches.
Since the multiplicities agree rather well, let us move on to more differential observables. Fig. 26 shows charged pion multiplicities as a function of the scaled rapidity y 0 , comparing the spectra obtained from SMASH to the experimental results of the FOPI collaboration, for Au+Au collisions at a kinetic energy of 1.5A GeV. SMASH reproduces the shape of the rapidity spectra fairly well, overestimating the total multiplicities by a few percent as seen before.
In Fig. 27 the proton rapidity spectrum yielded by SMASH is compared to FOPI measuremnts. The parameters are the same as just discussed. To get rid of spectators, all nucleons that interact only elastically have been ignored. SMASH does not model the production of nuclei formed by clustered nucleons. To be able to compare to the experimental data, a simple coalescence afterburner described in Appendix D has been employed. Any pairs of nucleons with momentum distance ∆p < 0.3 GeV and spatial distance ∆x < 0.9 fm have been ignored. These parameters were chosen to fit the data. The shape is very well reproduced at the tails, but the number of protons is overestimated at mid-rapidity.
In Fig. 28, the multiplicity of charged pions is shown as a function of the transverse mass m T for different windows of normalized rapidity y 0 , for C+C collisions at energies E kin ∈ {1, 2}A GeV as measured by the HADES collaboration [62]. For a purely thermal spectrum one would expect a straight line in the logarithmic plot, with the slope corresponding to the effective temperature. The events were generated with SMASH by sampling the impact parameter distribution provided by HADES (which was reconstructed using another transport model [62]). The calculations were performed with Skyrme and symmetry potentials, Fermi motion and Pauli blocking. It can be seen that SMASH describes the experimental data reasonably well. There are some deviations for large rapidities at 1A GeV and for small transverse mass at 2A GeV. In comparison to the UrQMD transport model [62], SMASH gives a similarly good agreement with the HADES data.
In Fig. 29, the multiplicity of nucleons and charged pions is shown as a function of transverse mass m T and rapidity y, for π − +C collisions at 1.7 GeV. The impact parameter was sampled from a minimum-bias distribution over the full range b ∈ [0, 3] fm. Like before, the SMASH simulation included potentials, Fermi motion and Pauli blocking. 20 test particles were used per real particle. Spectators (particles that only interact elastically) have been ignored. This scenario has been experimentally studied by the HADES collaboration, the results are however not yet public. The spike in the nucleonic rapidity spectrum corresponds to slow participants in the nucleus.
Unlike in experiment or in a simple thermal model, we can look at whole time evolution in a transport code. This enables us to study the different reaction rates. In Fig. 30 various forward (right) and backward (left) reaction rates are shown, as a function of time and as a total per event. Elastic nucleon-nucleon interactions are dominating and are divided by 10 in the plot. These are mostly due to the combination of Fermi motion and potentials, causing the nucleons in the target to interact with each other.
The first inelastic reactions are excitations of N * and ∆ * resonances. Production of ∆ resonances or elastic nucleon-nucleon collisions happen at later times. Producing N * and ∆ * consumes more pions than it directly yields, but these excitations decay mostly into ∆ and ρ, which finally produce pions again. It is remarkable that this system is far from chemical equilibrium, unlike symmetric collisions of heavy nuclei like gold or lead.
All in all, the transport approach presented here matches the experimental data on pion and proton production reasonably well, passes equilibrium/detailed balance tests and compares well to elementary particle production cross sections.     [58]. The experimental data (markers) is compared to the corresponding SMASH results (lines). N is the total number of pions obtained by integrating the spectrum. The normalized rapidity y0 = (y − ycm)/ycm was used. The impact parameter for the simulated events was sampled from the distribution given by the ERAT cuts corresponding to the experimental data (see Appendix C). The SMASH simulations were performed with potentials, Pauli blocking and Fermi motion. modifications of the cross sections, but there is still reasonable agreement in the current approach. Predictions for particle production in π-A collisions are made. In this case the meson-baryon interactions play a more dominant role than in heavy-ion reactions.
In the future, the approach will be enhanced to include the full strangeness production and the cross sections are going to be extended to higher energies by including string excitation and fragmentation. In addition, photon and dilepton production [63] is going to be studied in detail. Here it is of special interest to compare the non-equilibrium hadronic production with the one from thermal rates as currently employed in hydrodynamic approaches. In general, this approach will be very useful to study the effects of hadronic rescattering on flow and correlation observables at RHIC and LHC energies. Infinitematter calculations are going to be employed to study transport coefficients of hadronic matter as a function of temperature and baryo-chemical potential. Also the effects of kinematic cuts and baryon diffusion on higher moments will be investigated [64]. Overall, this approach constitutes a very flexible hadronic transport approach that is going to shed light on the properties of hot and dense strongly interacting matter as created in heavy-ion reactions in a large range of beam energies.    FIG. 28: Transverse mass spectra for pions measured by HADES in carbon-carbon collisions at 1 and 2A GeV [62]. The experimental data (markers) is shown for different longitudinal rapidity bins and compared to the corresponding SMASH results (lines). For readability, the data corresponding to each bin was multiplied with a different power of 10. The impact parameter distribution provided by HADES was used for sampling the events with SMASH.
[70] can be generated and VTK output [71] can be used to visualize the simulation. They were obtained from a SMASH simulation with 20 test particles per real particle, including potentials, Fermi motion and Pauli blocking. Spectators (particles only interacting elastically) were ignored. Data for this scenario has been measured by HADES, but is not yet published. The legend shows the total multiplicity N obtained from integrating the rapidity spectrum.
After this procedure, the remaining events should belong to the same centrality class as the experimental events. Note that ERAT is frame-dependent. For the purpose of this paper, it has been calculated in the fixed-target frame.

Appendix D: Nucleon clustering
A hadronic transport code does not have a concept of nuclei, because it considers only hadronic degrees of freedom. However when comparing to experiment, it is important to know which nucleons are bound in a cluster, because only unbound protons are considered as protons by the detector.
To model clustering we use a simple coalescence after-burner inspired by the work of Li et al. [73] that considers the pairwise distance in position and momentum space. Any pair of nucleons with a relative distance ∆r < r 0 and a relative momentum ∆p < p 0 is considered to be part of a cluster and will be ignored when calculating the nucleon spectra. To make this procedure Lorentz-invariant, before calculating the distances the particles are boosted to the center-of-momentum frame and their position is extrapolated so the boosted four-vectors correspond to the same time.
It is usually experimentally known how many protons are bound in a cluster, so the parameters (r 0 , p 0 ) can be chosen such that the correct multiplicities are obtained. Care has to be taken that the simulation runs long enough, otherwise r 0 strongly depends on the time at which the simulation is stopped. shows the total number of forward and backward reactions per event for various reactions. The same SMASH simulation results as in Fig. 29 were used.