Baryon-antibaryon annihilation and reproduction in relativistic heavy-ion collisions

The quark rearrangement model for baryon-antibaryon annihilation and reproduction ($B\bar B\leftrightarrow 3M$) - incorporated in the Parton-Hadron-String Dynamics (PHSD) transport approach - is extended to the strangeness sector. A derivation of the transition probabilities for the three-body processes is presented and a strangeness suppression factor for the invariant matrix element squared is introduced to account for the higher mass of the strange quark compared to the light up and down quarks. In simulations of the baryon-antibaryon annihilation and reformation in a box with periodic boundary conditions we demonstrate that our numerical implementation fulfills detailed balance on a channel-by-channel basis for more than 2000 individual $2 \leftrightarrow 3$ channels. Furthermore, we study central Pb+Pb collisions within PHSD from 11.7$A$ GeV to 158$A$ GeV and investigate the impact of the additionally implemented reaction channels in the strangeness sector. We find that the new reaction channels have a visible impact essentially only on the rapidity spectra of antibaryons. The spectra with the additional channels in the strangeness sector are closer to the experimental data than without for all antihyperons. Due to the chemical redistribution between baryons/antibaryons and mesons we find a slightly larger production of antiprotons thus moderately overestimating the available experimental data. We additionally address the question if the antibaryon spectra (with strangeness) from central heavy-ion reactions at these energies provide further information on the issue of chiral symmetry restoration and deconfinement. However, by comparing transport results with/without partonic phase as well as including/excluding effects from chiral symmetry restoration we find no convincing signals in the strange antibaryon sector for either transition due to the strong final-state interactions.


I. INTRODUCTION
Lattice Quantum-Chromo-Dynamics (lQCD) calculations suggest that at vanishing baryon chemical potential (µ B =0) there is a crossover phase transition from hadronic to partonic degrees of freedom [1][2][3][4][5][6] for the deconfinement phase transition as well as for the restoration of chiral symmetry. This leaves the open question whether or not a first-order phase transition might occur at finite baryon chemical potential implying a critical endpoint in the QCD phase diagram [7]. Since lattice calculations so far suffer from the fermion-sign problem, model-independent information on the QCD phase diagram can presently only be obtained from experimental data. It is thus expected that a thorough study of this issue with relativistic heavy-ion collisions (HIC) of different system sizes and at various bombarding energies will provide further information. However, the problem here is the model dependence in the interpretation of the measured particle yields and their relation to the properties of the fireball created in the collision [8]. Among the many observables suggested the strangeness enhancement was already proposed in the 80's of the past century [9] as a probe of the Quark-Gluon-Plasma (QGP). Particularly hyperons and antihyperons should provide an ideal sample for QGP fireballs since in the initial colliding nuclei no net strange quarks are present and a major part of the produced ss pairs should be produced by gluon fusion processes in the QGP [10]. However, due to a partial restoration of chiral symmetry close to the hadron-parton transition the ss production threshold is lowered in a dense hadronic medium and strangeness enhancement might also signal chiral symmetry restoration rather than deconfinement as suggested in Refs. [11][12][13][14]. We note that quark confinement and chiral symmetry breaking are not intimately connected at finite µ B [15]. The multi-strange baryons and antibaryons are expected to be more sensitive to the QGP than single-strange baryons or mesons since multi-strange baryons, and particularly multi-strange antibaryons, are suppressed by high hadronic energy thresholds as well as by long timescales for multi-step processes in a purely hadronic phase [11,16,17]. This, however, holds only for two-body production channels whereas three-body channels (e.g. by three vector mesons) do not suffer from severe energy thresholds. Accordingly, the reaction dynamics for baryon-antibaryon (BB) annihilation and recreation in the hadronic phase have to be under control before solid conclusions can be drawn on the boundary in the QCD phase diagram or on freeze-out conditions in relativistic heavy-ion reactions. Furthermore, all strangeness exchange channels in the hadronic phase have to be taken into account as pointed out in Refs. [18,19].
A first step in this direction has been taken in Ref. [20] where the three-body fusion of nonstrange pseudoscalar and vector mesons to BB pairs has been incorporated in the Hadron-String Dynamics (HSD) transport approach [21] that preferentially describes the hadronic phase.
Here the matrix element squared has been extracted from the experimental data on pp annihilation and the threebody meson channels have been determined on the basis of detailed balance. It was found that in central collisions of heavy nuclei the annihilation of antinucleons is almost arXiv:1710.00665v2 [hep-ph] 23 Jan 2018 compensated by the inverse recreation channels. If this holds true also in the strangeness sector is presently unknown. Furthermore, the former HSD calculations did not incorporate a deconfinement phase transition to the QGP nor effects from chiral symmetry restoration and thus did not allow to draw any conclusions on the phase boundary of QCD. On the other hand the HSD transport approach has been further extended in the last 15 years a) to the formation of an initial partonic phase with quark and gluon quasiparticle properties that are fitted to lattice QCD results in thermodynamic equilibrium, b) to a dynamical hadronization scheme on the basis of covariant transition rates, c) to incorporate further hadronic reactions in the strangeness sector with full baryon-antibaryon symmetry and d) to employ essential aspects of chiral symmetry restoration in the hadronic phase [14]. Whereas the latter developments are important for the lower Super Proton Synchrotron (SPS) energy regime to account for the strangeness enhancement seen experimentally in heavyion collisions, the formation of a partonic phase is mandatory to understand the physics at higher SPS, Relativistic Heavy-Ion Collider (RHIC) and Large Hadron Collider (LHC) energies. Since multistrange baryons and antibaryons at top SPS energies no longer stem from string fragmentation (as in HSD [20]) but preferentially from hadronization at energy densities around 0.5 GeV/fm 3 the issue of three-meson fusion reactions for the formation of baryon-antibaryon (BB) pairs and the annihilation of BB pairs to multiple mesons has to be reexamined. In this work we will present, furthermore, the extension of the quark rearrangement model (QRM) for baryon-antibaryon annihilation and recreation to the strangeness/antistrangeness sector (briefly denoted by SU (3)). We will show the impact of these additional reaction channels for heavy-ion collisions using the Parton-Hadron-String Dynamics (PHSD) transport approach to simulate central Pb+Pb collisions in the bombarding energy regime from 11.7A GeV to 158A GeV. The PHSD [22][23][24], which incorporates in addition to HSD a transition to the partonic phase as well as dynamical hadronization, reproduces many observables for p+p, p+A and A+A collisions ranging from SPS up to Large Hadron Collider (LHC) energies [14,[25][26][27]. Since PHSD is found to also well describe the spectra of strange mesons and baryons from heavy-ion collisions from 2A GeV up to RHIC/LHC energies when incorporating aspects of chiral symmetry restoration in the hadronic phase [14], its performance in the strange antibaryon sector will be tested using the extended QRM and also lead to predictions for rare multi-strange baryons and antibaryons in the lower energy regime where experimental data are scarce or lacking at all. This work is organized as follows: In Sec. II we recapitulate shortly the ingredients of PHSD while in Sec. III we briefly recall and motivate the quark rearrangement model for baryon-antibaryon annihilation and recreation (BB ↔ 3M ). We extend the QRM to the strangeness sector and introduce a strange quark suppression factor for the transition matrix element squared in the strangeness sector. After deriving the transition probabilities on the basis of detailed balance, we present in Sec. IV the validity of our numerical implementation for detailed balance in case of BB ↔ 3M reactions including the strangeness sector within simulations in a finite box with periodic boundary conditions. In Sec. V we present results for antibaryons and multi-strange baryons from PHSD simulations for central Pb+Pb collisions in the SPS energy regime and study the impact of chiral symmetry restoration and deconfinement. We will compare simulations using the baryon-antibaryon annihilation and formation with and without the strangeness sector with each other and to available experimental data for rapidity and transverse mass spectra. Furthermore, we compare the PHSD results for central Pb+Pb reactions at 40A GeV with those from the Ultra-relativistic Quantum Molecular Dynamics model (UrQMD) [28] and the threefluid dynamics model (3FD) using a 2-phase equation of state [29]. We conclude our study with a summary in Sec. VI while more technical details are described in the appendices.

II. THE PHSD TRANSPORT APPROACH
The PHSD is a microscopic covariant transport approach for strongly interacting systems which is based on Kadanoff-Baym equations [30][31][32][33] for the Green's functions in phase-space representation in first order gradient expansion [34,35]. Due to its basis on the Kadanoff-Baym equations it can describe systems in and out-of equilibrium and goes beyond the quasiparticle approximation by incorporating dynamical spectral functions for the partons. It is capable of describing the equilibration process of systems which are far out-of equilibrium to the correct equilibrium state [36]. The PHSD incorporates a partonic as well as a hadronic phase to describe all stages of a relativistic heavy-ion collision with transitions from strings to dynamical partons as well as dynamical hadronization. In the hadronic phase high-energy resonance decays are described by multi-particle string decays. PHSD is capable of simulating the full time evolution of a relativistic heavy-ion collision -from impinging nuclei in their 'groundstates' to the final hadronic particles -ranging from SchwerIonen-Synchroton (SIS), Alternating Gradient Synchrotron (AGS) over Facility for Antiproton and Ion Research (FAIR)/ Nuclotron-based Ion Collider fAcility (NICA) up to Relativistic Heavy-Ion Collider (RHIC) and Large Hadron Collider (LHC) energies and is able to reproduce a large number of observables in these energy regimes for p+p, p+A and A+A reactions [14,25]. The properties of the off-shell partonic degrees-offreedom are determined by the Dynamical-Quasi-Particle-Model (DQPM) [37,38] which provides the masses, widths and spectral functions of the dynamical gluons and quarks/antiquarks [39,40]. The (essentially three) parameters of the DQPM are chosen to reproduce the lQCD equation-of-state at vanishing baryon chemical potential. It has been shown that using PHSD in a box with periodic boundary conditions it reproduces the lQCD results for transport coefficients such as the shear and bulk viscosity as well as the electric conductivity for the partonic phase [41][42][43]. In the PHSD simulation of a nucleus-nucleus collision the primary hard nucleon-nucleon scatterings produce strings which are color-singlet states described by the FRITIOF Lund model [44]. As the strings decay they produce "prehadrons" that have a formation time of τ f ≈ 0.8 fm while "leading hadrons", which originate from the string ends, may interact instantly without formation time but with reduced cross-sections in line with the constituent quark model [21]. If the local energy density is above the critical value of c ≈ 0.5 GeV/fm 3 , as provided by lQCD calculations, the unformed hadrons dissolve into dynamical quarks with properties defined by the DQPM at given energy density. In the partonic phase these partons propagate in the scalar self-generated mean-field potential and scatter with each other with cross sections extracted from the dynamical widths of partons. The expanding system then leads to a decreasing local energy density until it is close to or below the critical value c . At this point the partonic degrees-of-freedom hadronize to colorless offshell mesons and baryons by the fusion of massive quarkantiquark pairs or the fusion of three quarks (antiquarks) conserving energy, three-momentum and quantum numbers in each event [22]. In the hadronic phase -as found in the corona and at late reaction times -the hadrons interact with each other in elastic and inelastic collisions with cross sections taken from experimental data or evaluated within effective hadronic Lagrangian models. The detailed balance relation for each reaction channel is incorporated and ensures the correct backward reaction rates. In particular the strangeness exchange reactions are included in meson-baryon/antibaryon, baryonbaryon and antibaryon-antibaryon collisions following Refs. [45][46][47]. Furthermore, retarded electromagnetic fields as generated by the electric charge currents (from charged hadrons and quarks) are incorporated [48].

III. QUARK REARRANGEMENT MODEL FOR B +B PRODUCTION AND ANNIHILATION
In this section we present the quark rearrangement model and the most relevant equations for the two-and threebody scattering rates. An extensive description for the invariant reaction rates for general particle number changing processes as well as the motivation for the quark rearrangement model is given in Ref. [20].  [20].

A. Concept
As discussed in Ref. [20] one experimentally finds a dominant annihilation of pp into 5 pions at invariant energies 2.3 GeV ≤ √ s ≤ 4 GeV, see Fig. 1. The final number of 5 pions may be interpreted as an initial annihilation into πρρ with the ρ mesons decaying subsequently into two pions each. The channel ππρ then leads to 4 final pions, the channel πωρ to 6 final pions, the channel ρωρ to 7 final pions etc. Accordingly, the baryon-antibaryon annihilation in the first step is a two-to-three reaction with a conserved number of quarks and antiquarks. This is the basic assumption of the quark rearrangement model which is also illustrated in Fig. 2. The annihilation reaction pp → πρρ is the dominant process in pp annihilation for invariant masses below 4 GeV, typical for the hadronic phase of a heavy-ion collision. By allowing the mesons M i to be any member of the 0 − or 1 − nonets one can describe an arbitrary BB annihilation and recreation by rearranging the quark and antiquark content. An implementation of baryon-antibaryon annihilation in such a manner misses the annihilation into one or two mesons, however, higher numbers of final mesons are implemented through the subsequent decay channels. This approach gives a realistic description for pp annihilation and we assume that for other baryon-antibaryon pairs than pp a similar annihilation pattern holds. Since there are no measurements of annihilation cross sections other than pn and pp this is our best guess which might be falsified by experiment.

B. Covariant transition rates
The quark rearrangement model only contains reactions of the kind 2 ↔ 3. The detailed balance based Lorentz invariant on-shell collision rate for the reaction BB → 3M in a volume element of size dV and time-step size dt is written as [20]: In (1) c denotes all BB pairs with the properties c = (m c 1 , m c 2 ; ν c ); c are all the possible meson channels with c = (m c 3 , m c 4 , m c 5 ; λ c ), with m being the masses of the respective particles, and ν and λ the quantum numbers signifying the channel (charge, parity, spin and strangeness). We assume that the transition matrix element squared W 2,3 does not significantly depend on the outgoing momenta and just on the invariant mass of the reaction, which holds approximately true for pp as we will see later. A formulation based on the matrix element will ensure detailed balance. The on-shell n-body phase-space integral is defined by and in case of a constant transition matrix element dominates the interaction rate of the system. The factor N c fin is the multiplicity of the meson triple c and results from the summation over the spin s and possible isospin projections F iso compatible with charge conservation of the meson channel c: N c fin = (2s 3 + 1)(2s 4 + 1)(2s 5 + 1) The division by N id !, with N id denoting the number of identical mesons, ensures that each charge configuration is only considered once for a given meson triple. The functions f are the distribution functions of the BB pair in momentum and coordinate space. When looking at a specific BB pair one has to make sure that only meson channels are considered which conserve charge, energy and parity. The probability of this specific BB pair c to annihilate into any of these possible meson channels c is related to the total annihilation cross section of the BB-pair σ c ann [49]: where dV and dt are taken finite. The probability for a specific final stateP c →c in case of an annihilation is then given by the available phase space and the multiplicity of all possible meson channels c: In a similar manner one finds for the probability of a specific meson channel c fusing together and forming a specific BB pair c , with N c B = (2s 1 + 1)(2s 2 + 1) denoting the multiplicity of the BB pair. A more detailed derivation of these formulae is given in Appendix A.

C. Annihilation cross sections
For the calculation of actual collision probabilities, Eq. (4) and (8), we are still missing the cross sections. As already mentioned above we assume the cross sections to depend only on the invariant energy, not the outgoing momenta. This assumption is approximately fulfilled for pp and pn annihilation, see Fig. 3. Other channels have not been measured so far. Since there are no experimental data available we assume a similar behavior for different spin combinations like p∆. In this work we investigate in particular the strangeness sector. We model the cross sections of particles with strangeness by where ς andς are the number of strange and antistrange quarks in the BB pair c and λ ∈ [0, 1] is a factor suppressing the transition matrix element for particles taking part in the quark rearrangement model and effectively suppressing the cross section. This parametrization is motivated by PYTHIA [50] simulations where one sees a similar suppression for particles with strangeness compared to non-strange particles at the same energy above threshold. In the final implementation in PHSD the suppression factor has the value λ=0.5 which is in rough agreement with the PYTHIA simulations embedded in PHSD. We choose a dependence on not just the net strangeness S but the sum of strange and antistrange quarksς + ς due to their higher mass and a subsequent suppression of the rearrangement. The implementation with the strangeness |S| = |ς − ς| instead ofς + ς has no practical influence on the final results in case of relativistic heavy-ion collisions (cf. Appendix E).

IV. SIMULATIONS IN A FINITE BOX
This section addresses the implementation of the 2 ↔ 3 reactions formulated above in PHSD and checks the consistency of the numerical implementation. We use transport simulations in a box with periodic boundary conditions to investigate the behavior of the quark rearrange-  ment model in equilibrium. We recall that in equilibrium -according to detailed balance -the reaction rate for BB → 3M should be the same as for 3M → BB. Furthermore, for a consistent implementation detailed balance should not only be fulfilled for the sum of all reaction channels but on a channel by channel basis. In the box simulations only the BB ↔ 3M reactions are considered now and all particles are taken as stable such that no decays occur. The particles incorporated are the 0 − , 1 − meson nonets, the 1/2 + baryon octet and the 3/2 + baryon decuplet. Additionally, we consider N(1440) and N(1535) baryonic resonances. Furthermore, we take into account the strangeness content of η and φ with 50% and 83.1% ss content, respectively. With this the number of possible mass channels amounts to more than 2000. Hence, an initialization with every possible channel is not feasible. Therefore, we look at systems which are initialized by a single type of baryon and antibaryon adding up to 100 systems for the consistency check. The box simulations have the following initial conditions: • Box volume lies around 18000 fm 3 with periodic boundary conditions • The initial momentum distribution is of Boltzmann-shape • For the box simulations a suppression of channels including strangeness is neglected.
Before the actual calculations the three-body phase-space integrals R 3 have been calculated and fitted by proper functionals to save enormous CPU time. In detail, the three-body phase-space integrals R 3 , depending on invariant energy and three masses m i , with R 2 defined in Eq. (B10) are fitted by with t = √ s − m 1 − m 2 − m 3 and a i > 0. The fit parameters a i have been evaluated for each combination of meson masses m 1 , m 2 , m 3 and stored on file. For further details on the phase-space integrals we refer the reader to appendix B.
We recall that the fusion of three mesons can not be described in a Lorentz invariant way by geometrical collision criteria between the particles due to the three inertial systems. To find a solution we employ the in-cell method introduced by Lang et al. [52] and adopted in Ref. [20]. This method is also employed for 2 ↔ 3 reactions in partonic cascade calculations [53]. The in-cell method can be used for any number of colliding particles since there is no problem with time ordering due to the locality of the formulation. In the in-cell method spacetime is divided into four dimensional cells with widths ∆x, ∆y, ∆z, ∆t and only particles inside the same cell may collide with each other. One calculates the reaction probabilities of each particle with every other one inside the same cell. The actual collision and the final state is chosen via Monte-Carlo. The possible final states and multiplicities in Eqs. (4) and (8) are precalculated to save computational time during the transport simulation. The cell size and the time step ∆t are optimized for the problem under investigation such that the total probability of a transition in a local cell does not exceed unity but is also not too small. For the actual calculations shown below we use dt = 4 fm/c and dV =40 fm 3 which ensures that the transition probabilities are always below unity. We now discuss results for a few selected systems. We  present randomly picked ensembles that cover the qualitative range of possible systems, i.e. systems consisting of only initial light quarks, only initial strange/antistrange quarks as well as a variety of combinations of light and strange quarks/antiquarks. In Fig. 4 the time evolution of the particle densities for a system initialized with protons and antiprotons is shown to demonstrate the production and annihilation of different particle species in a system consisting initially only of protons and antiprotons. After the first timestep of the simulation a lot of new mesons like pions, ρ and ω mesons are formed. At later times also strange mesons and baryons are formed because of the partial ss content of φ and η. In equilibrium the system has a significant amount of mesons and baryons with strange and antistrange quarks. However, the generation of strange quarks even for the meson sector takes a long time (≈ 60 fm/c) to produce significantly high strange particle densities; thus the generation via φ and η should have negligible influence on actual heavyion collisions since large densities are needed for a significant contribution from the meson fusion. In a 5% central Pb+Pb collision at 158A GeV the meson fusion dies out at ∼ 13 fm/c which is insufficient for having a major influence on the strangeness sector, see Fig. 7 (discussed in section V below).
We show in Fig. 5 the total reaction rate as a function of time for two exemplary initializations which were initialized with p +p and ∆ 0 +Λ, respectively. Both systems share a similar evolution of the total reaction rate. All systems reach detailed balance much faster (≈ 40 fm) than they reach equilibrium (≈ 1000 fm).
Detailed balance should also be valid for the total reaction rate as function of the invariant mass. For this we show in Fig. 6 the total reaction rate as a function of the invariant mass √ s in the plateau region of Fig. 5 which is associated with the equilibrium state. From Fig. 6 we see that detailed balance is also fulfilled for this quantity. Note that the maximum achievable invariant mass of particles participating in annihaltion or recreation (in equilibrium) is lower in systems initialized with lighter baryons than for systems initialized with heavier ones. The last most crucial check for detailed balance is the fulfilment on a channel by channel basis. To this end we define the deviation from detailed balance for each channel by We calculate δ for each of the more than 2000 channels and look at the channels with the largest reaction rates in all 100 investigated systems. In Tab. I the 10 most important channels with the largest reaction rates are shown from highest to lowest for 3 of the exemplary systems as well as the average for all 100 investigated systems and the average over all channels. The average over all 100 investigated systems shows that detailed balance is fulfilled better than 97% on a channel-by-channel basis for the 100 most dominant channels. This verifies the correct implementation of the baryon-antibaryon annihilation and recreation within the quark rearrangement model in the PHSD transport approach. Some channels of a system may deviate by more than 5% from detailed balance, however, this is a relict of too low statistics. We found only few channels (≈ 20 for the 10 most dominant channels) that had a deviation of up to 9%. In general these deviations may be neglected as can be seen in the averaged values and the dominant number of channels being very close to detailed balance which gives a proof for the working principle of the implementation presented.

V. PHSD SIMULATIONS FOR HEAVY-ION COLLISIONS
In this section we show the influence of the additional channels in the strangeness sector for BB ↔ 3M reactions on heavy-ion collisions in the energy regime of 11.7-158A GeV. Before coming to the actual results we compare in Fig.  7 the reaction rate for the total baryon-antibaryon annihilation (solid line) and formation (dashed line) from PHSD in 5% central Pb+Pb collisions at 158A GeV. Whereas the meson-fusion rate dominates at early times (< 13 fm/c) the annihilation takes over for larger times during the final expansion of the system. Although the time integrals of both rates are about the same there is no appreciable time interval in which both rates are identical. This indicates a strong nonequilibrium dynamics of baryon-antibaryon annihilation and reproduction in actual heavy-ion reactions.
We note that a similar analysis has been performed in the earlier study in Ref. [20] (Fig. 7) on the basis of the HSD transport model (version 2.3) for the same system, however, without averaging over the ensembles. The earlier rates differ substantially from the present results from PHSD (version 4.0) due to the different degrees of freedom in the initial phase of the collision. In order to quantify the differences we have recalculated the rates within HSD2.3 (from the year 2002) and compared the numbers with those from PHSD4.0, which is the most recent version including also the effects from chiral symmetry restoration [14] (PHSD3.3) and nonperturbative charm dynamics as well as extended 2 ↔ 3 reactions. We found that both rates (from HSD2.3 and PHSD4.0) differ only slightly for times ≥ 6 fm/c (after contact of Pb+Pb at b=2fm) but the huge rates (from HSD2.3) at the first few fm/c are essentially missing in PHSD4.0. This is due to the fact that at the top SPS energy the initial energy conversion goes to interacting partons in PHSD4.0 and not to strings decaying to hadrons (and partly to BB pairs) in HSD2. 3 sities are sizeably lower. In this dilute regime the 3-body channels first dominate and decrease fast in time whereas the 2-body annihilation reactions still continue for some time. As addressed in the Introduction we thus expect also differences in the antibaryon rapidity spectra as compared to the early results from HSD2.3 [20]. However, in both transport calculations -incorporating the 2 ↔ 3 reactions -the time integrated rates for annihilation and reproduction turn out to be about equal. The actual PHSD calculations for relativistic nucleusnucleus collisions are carried out in the parallel ensemble method, i.e. in case of the cascade mode a typical number of 100 -300 ensembles are propagated in time fully independent from each other. However, the calculation of net-baryon densities, scalar densities and energy densities -needed for the full PHSD dynamics -is carried out by averaging over all ensembles. This results in a crosstalk between ensembles due to the propagation of particles in the self-generated mean fields (for partons and baryons/antibaryons) as well as in the baryon/antibaryon formation in the hadronization. A systematic study of all particle spectra in rapidity and transverse mass shows that the results for mesons and baryons well scale with the number of ensembles whereas the antibaryon sector shows small variations with the number of ensembles. This scaling violation is essentially due to the numerical approximations that have to be presently introduced in order to keep the huge number of reaction channels manageable. This introduces a systematic error in our calculations for the antibaryon sector which is accounted for by hatched bands in the following figures. The solid or dashed lines correspond to the standard ensemble number of 150 used as default in PHSD calculations in the energy range of interest.

A. Rapidity and transverse mass spectra
We now discuss the influence of the BB ↔ 3M reactions on observables measured in actual experiments from 11.7 -158A GeV. We first focus on rapidity spectra and mention that the BB ↔ 3M reactions have practically no influence on baryon and meson spectra [14] and, hence, we only show the results for the relevant antibaryons and Ξ − to demonstrate that the influence on baryons is barely visible. For results on meson and baryon spectra we refer the reader to the review [25] and Ref. [14]. As mentioned above the full, dashed and dotted lines show the results for 150 ensembles; the blue and red hatched areas result when employing different ensemble numbers in a wide range. We first focus on the influence of the newly incorporated strangeness sector. In the following, we compare the implementation with only light quark channels (SU(2)) with the new one including also the strangeness sector (SU (3)). The rapidity spectra ofp,Λ +Σ 0 , Ξ − ,Ξ + , Ω − + Ω + for central Pb+Pb collisions from 11.7 to 158A GeV are shown in Figs. 8 and 9. The rapidity spectra of  (3)) while the dashed lines results from discarding strange or antistrange quarks in the reaction channels (denoted by SU (2)). The error bands indicate the systematic uncertainty of the calculations due to a different ensemble size. The dotted lines show the results with BB ↔ 3M reactions switched off. The data points are taken from Refs. [54][55][56].
the anti-hyperons are overall closer to the experimental data when taking into account the strangeness sector for the BB ↔ 3M reactions. However, the antiproton spectra are faintly influenced by the incorporated sector and describe the data only moderately well. In general the investigations suggest that the BB ↔ 3M reactions have the largest impact at energies below 80A GeV. This result shows that the consideration of the strange quarks helps improving the description of a heavy-ion collision in the framework of PHSD. For particles likeΞ + , Ω − andΩ − at lower energies, where currently no experimental data are available, our results should be taken as predictions.
In Fig. 8 we, furthermore, show results from calculations neglecting the BB ↔ 3M reactions. We find that the rapidity distribution forp has a higher peak and is narrower compared to calculations with BB ↔ 3M , while the total number of antiprotons is about the same. The results for the antihyperons -starting from 20A GeVlie on top of the SU(2) simulations. At 11.7A GeV the hyperon spectra are closer to the SU(3) calculations and for Ω − +Ω + lie even below those. Another interesting observable measured in experiment is the transverse mass (m t ) spectrum at midrapidity, i.e. dN/(m t dydm t ) as displayed in Fig. 10. Here the additional strangeness sector has qualitatively the same impact as for the rapidity spectra. Accordingly, we only show results for central Pb+Pb collisions in the energy regime from 20 to 158A GeV including the strangeness sector for the BB ↔ 3M reactions. For the Ξ − we find that PHSD describes the low m t regime for energies below 158A GeV rather well. However, for higher m t the data points are missed due to a harder experimental slope of the spectrum. At 158A GeV some Ξ − 's are missed in the low m t regime. TheΛ +Σ 0 spectrum is close to the experimentally measured data for all energies, however, at 158A GeV it falls off too fast. The transverse mass spectra of the antiprotons are overall in very good agreement with experiment, the only drawback is the overproduction at midrapidity which is most visible for 20 and 30A GeV. Also, theΞ + are in close vicinity to the experimental data for energies smaller than 158A GeV, but fall off too quickly at 158A GeV. The production of Ω − and Ω + was underestimated already in the rapidity spectra, see Fig. 9, but looking at the transverse mass spectra at 158A GeV the results are in reasonable agreement with experiment for m t < 0.8 GeV.

B. Impact of chiral symmetry restoration and deconfinement
We now address the question with respect to traces of chiral symmetry restoration and deconfinement in the antibaryon and multi-strange baryon spectra from central Pb+Pb collisions at SPS energies. We recall that clear signals have been found before in the strange meson and baryon rapidity distributions [12,14] and one might speculate if a similar signal can be seen in the antibaryon sector. To this aim we perform transport calculationsincluding the BB ↔ 3M channels specified above -with different settings: • HSD calculations without chiral symmetry restoration (CSR) and deconfinement since HSD does not include a partonic phase • HSD calculations with chiral symmetry restoration (CSR) in the hadronic phase but without deconfinement • PHSD calculations without chiral symmetry restoration (CSR) in the hadronic phase but with a deconfinement transition • PHSD calculations with chiral symmetry restoration (CSR) in the hadronic phase and with a deconfinement transition.
The systems addressed are central Pb+Pb collisions at 30 and 158A GeV. The rapidity spectra for antibaryons and Ξ − are displayed in Fig. 11 and show that at 158A GeV the impact of chiral symmetry restoration is very small in the HSD calculations (without deconfinement) as well as for PHSD (including deconfinement) except for thē Λ +Σ 0 spectrum. When comparing HSD and PHSD results including CSR we find a slight reduction of thep spectra, a moderate enhancement for theΛ+Σ 0 spectrum and only a small enhancement for Ξ ± and Ω − +Ω + when including a partonic phase. Since the reproduction of the multistrange sector by PHSD is very poor one cannot conclude on the presence of a deconfinement transition on the basis of the rapidity spectra shown in Fig. 11. Note, however, that a clear signal has been found in the elliptic and triangular flow before in Ref. [27] at this energy. At 30A GeV the situation is not much better. The PHSD calculations with CSR perform best for Ξ − andΞ + , however, overestimate thep andΛ +Σ 0 yield. The HSD calculations are too low in the strange antibaryon sector including/excluding CSR providing some hint that a partonic phase should be present in a moderate space-time volume at this energy. Accordingly, the antibaryons and in particular the multi-strange sector do not give additional information on chiral symmetry restoration or deconfinement within the framework of PHSD calculations.

C. Comparison to other dynamical models
In this subsection we compare our current PHSD results to those from other dynamical models which have been employed for heavy-ion reactions in the SPS energy regime, in particular from the Ultra-relativistic Quantum Molecular Dynamics model (UrQMD) [59,60] and the three-fluid dynamics model (3FD) [61]. The UrQMD is a hadronic transport model including a multitude of hadronic resonances as well as strings that are responsible for multi-particle production. The 3FD is a fluid dynamical model describing -within the framework of hydrodynamics -the transition from the initial baryonic fluids (projectile and target) to the newly produced fluid (around midrapidity). For details we refer the reader to the original literature [59][60][61]. We show in Fig. 12 our actual results in case of the rapidity spectra for a central Pb+Pb collision at 40A GeV with the BB ↔ 3M reactions including the strangeness sector in comparison to results from the UrQMD [28] and the 3FD using a 2-phase equation of state [29]. The 3FD model, like PHSD, overshoots the antiproton yield whereas UrQMD is close to the experimental data. TheΛ +Σ 0 spectrum is described by PHSD and the 3FD model similarly close to the experimental data whereas UrQMD produces too few. For the Ξ − all models show different behav- iors; whereas the 3FD model overpredicts the production, PHSD produces slightly too few Ξ − at midrapidity but describes otherwise the shape well. UrQMD predicts (just like forΛ +Σ 0 andΞ + ) too few antibaryons since BB annihilation is incorporated, however, not the backward channels thus violating detailed balance. PHSD and the 3FD model are close to the experimental data forΞ + , with the 3FD slightly underpredicting the yield. Depending on the particle species of interest one model describes some yield better than the other at higher SPS energies. In general, the 3FD model and PHSD appear to be similarly capable of roughly describing the dynamics of baryons and antibaryons with strangeness content in this energy range.

VI. SUMMARY
In this work we have recapitulated and extended the quark rearrangement model for baryon-antibaryon annihilation (BB ↔ 3M ) in the course of heavy-ion collisions. The approximate validity of this model was motivated by the distribution in the number of final state pions in pp annihilation for 2.3 GeV≤ √ s ≤ 4 GeV (cf. Fig. 1), where the 3-body channel πρρ e.g. leads to 5 pions (on average) in the final state. Additionally to the HSD calculations in Ref. [20], we have included in the 2 ↔ 3 channels the strangeness sector with a suppression factor for the matrix elements of particles having strange and anti-strange quarks. We have shown, using simulations in a box with periodic boundary conditions, that the numerical implementation of the quark rearrangement model including the strangeness sector satisfies the detailed balance 2 ↔ 3 relation on a channel-by-channel basis as well as differentially as a function of the invariant energy √ s. We found that the earlier rates from HSD2.3 [20] differ substantially from the present results from PHSD (version 4.0) due to the different degrees of freedom in the initial phase of the collision. Both rates (from HSD2.3 and PHSD4.0) differ only slightly for times ≥ 6 fm/c (after contact of Pb+Pb at b=2fm) but the huge rates (from HSD2.3) at the first few fm/c are essentially missing in PHSD4.0. This is due to the fact that at the top SPS energy the initial energy conversion goes to interacting partons in PHSD4.0 and not to strings decaying to hadrons (and partly to BB pairs) in HSD2.3. Thus in PHSD4.0 (at the top SPS and higher energies) there are initially no BB pairs that might annihilate nor mesons that might fuse! Due to the very high hadron densities in HSD2.3 (after string decay) both the annihilation and reproduction rates are very high and about equal whereas in the hadronic expansion phase the densities are sizeably lower. In this dilute regime the 3-body channels first dominate and decrease fast in time whereas the 2body annihilation reactions still continue for some time. However, in both transport calculations -incorporating the 2 ↔ 3 reactions -the time integrated rates for annihilation and reproduction turn out to be about equal. The influence of the newly implemented channels in the strangeness sector on actual heavy-ion collisions has been investigated in PHSD simulations (version 4.0) of central Pb+Pb collisions from 11.7-158A GeV. The rapidity spectra of antibaryons -using the quark rearrangement model with and without the strangeness sectorhave been compared to experimental data where available. Changes could only be seen for the antibaryons in the investigated energy regime whereas the meson and baryon sector are practically unchanged [14]. Due to the chemical rearrangement between the baryons and mesons considered in BB ↔ 3M reactions an overall higher antiproton production was observed for all energies which pushed the PHSD results up thus overestimating the experimental data. The other antibaryons got closer to the experimental data when the strangeness sector was included for the BB ↔ 3M reactions. For the energies investigated the strangeness sector has the largest impact at the lowest energies. The results show that the BB ↔ 3M reactions indeed need the strangeness sector to describe the heavy-ion collisions more properly. We note, however, that the quark rearrangement model might be too crude to allow for robust conclusions. We still need experimental information on baryonantibaryon annihilation cross sections other than pp and pn to achieve a better description and understanding of heavy-ion collisions. In addition to the rapidity spectra, we have shown the transverse mass spectra for various antibaryons and have seen that the low m t region is well described for all antibaryons with the exception of the antiprotons that are overpredicted at energies lower than 80A GeV. For higher transverse masses some spectra fall off too fast thus underestimating the experimental data to some extent. Accordingly, our understanding of antibaryon dynamics is far from being complete and we might still miss essential ingredients. We have additionally addressed the question if the antibaryon spectra (with strangeness) from central heavyion reactions at SPS energies provide further information on the issue of chiral symmetry restoration and deconfinement. By comparing results from HSD (without partonic phase) with those from PHSD (with partonic degrees of freedom) as well as including/excluding effects from chiral symmetry restoration (Fig. 11) we did not find convincing signals for either transition due to the strong final-state interactions.
where N c B denotes the multiplicity of the final state and the two-body phase-space integral R 2 is given by with λ defined in Eq. (5). The transition matrix element squared W 2,3 is not known but using Eq. (9) for our special problem of 2 ↔ 3 processes one gets, where we have taken W 2,3 out of the sum over the baryonantibaryon pairs and end up with the expression for the normalisation constant for the invariant energy √ s from Eq. (7). Inserting Eq. (A3) for the transition matrix element squared into (A1) gives the result for the transition probability for the meson fusion in Eq. (8). Note that all energies and momenta in the calculations of the transition probabilities are in the lab frame.

Appendix B: Phase-space integrals
The on-shell phase-space integrals occurring throughout this work inhibit most of the dynamics of the system. As they play a major role this section is dedicated to some more details of phase-space integrals. We recall that the n-body phase-space integral is generally defined by with ρ denoting the spectral function of the respective particle. Since the phase-space integrals are Lorentz invariant we will always work in the center-of-mass system. In the on-shell case the spectral function takes the form To show (as an example) the behavior of the different nbody phase-space integrals it is instructive to look e.g. at the consecutive decays pp → πρρ → 3πρ → 5π which are essentially the motivation for the QRM. Also, this example connects the 3-,4-and 5-body phase-space integrals as a function of the invariant energy above threshold (see below).
For the sake of completeness, we start with the 1-body phase-space integral, where E is the on-shell energy E = m 2 + p 2 and the mass m of the particle is equal to the invariant energy √ s. This result shows that the 1-body phase-space decreases with increasing √ s. The 2-body phase space can also be  evaluated analytically, The zeroes of the delta function are given by where only the positive value has to be taken in our calculation. Rewriting the delta function as and plugging Eqs. (B8) and (B9) into Eq. (B7) we obtain the two-body phase-space integral with E 1 + E 2 = √ s from the original delta function. The typical shape of R 2 ( √ s, m 1 , m 2 ) is shown in Fig. 13 for the masses m 1 = 1 GeV and m 2 = 2 GeV as a function of the invariant energy above threshold. The upper limit is independent of the masses and is given by 1/8π. The on-shell three-body phase-space integral is the most important one for our work and a good example for the evaluation of phase-space integrals of higher order since the n-body decay can be considered as consecutive 2-body decays, see Fig. 14 for an illustration. Note that in Fig. 14 k n = p and k 1 = p 1 . A prerequisite in calculating the phase-space integral is that we do not have any incoming momenta in between the first and final 2-body decay. For the calculation of the process we employ the recursion relation for phase-space integrals, and also insert two identities The first identity from Eq. (B12) gives the mass of the first cluster from which the 4-momentum p n splits. The second identity ensures energy-momentum conservation in the splitting process. Plugging both identities into Eq.
(B11) we find With this expression any n-particle phase-space integral can be calculated in a straight forward fashion as long as the masses m i are known. Note that the last R 2 , which one gets after applying Eq. (B15) several times, has no additional factor 1/(2π). In Fig. 15 the phase-space integrals for 3, 4 and 5 particles are shown as a function of the invariant energy above threshold for our example of initial πρρ with a subsequent decay into 3πρ and a final decay to 5 pions. All phase-space integrals share a similar shape, only the magnitudes close to threshold vary substantially with the number of particles.
Appendix C: In-cell method: cell-size dependence We here show the stability of our approach with respect to the equilibrium state when changing the size of the cells. For this investigation we keep the time step dt constant but enhance the cell volume ∆V by 20% and compare the reaction rate as a function of time to the default calculations in Fig. 16. We observe that the change in the cell size does not have any impact on the equilibration at all. For all times both cell sizes produce the same results giving testimony to the stability of the numerical implementation.
Appendix D: In-cell method versus next-neighbor interaction The in-cell method used for the description of the BB ↔ 3M reactions has been implemented cutting effectively the space-time into cells of cell-size ∆V × ∆t and letting only particles of the same cell interact with each other. Another possibility for the implementation of the baryonantibaryon annihilation (and recreation) is by defining the volume ∆V by a sphere around the first particle and letting all particles in the sphere interact with each other; this implementation we denote by next-neighbor (NN) algorithm in the following. In Fig. 17 we compare the results of these two choices. Due to the large finite size effects for the NN method the volume of the box had to be enhanced and filled with the same density as the standard box but letting only the particles inside the standard box volume be the particles from whose sphere the partners are selected. After employing this minimization of finite size effects we find that both methods give the same reaction rates for times larger than ≈30 fm. A small deviation between both methods is seen for smaller times. As expected one might use in general also the NN method. The disadvantage of the numerical implementation of the NN method is the larger computational time in comparison to the discretization of space-time. Thus PHSD uses the in-cell method not for the individual cells from the NN method but for the fixed cells of the spacetime discretization.

Appendix E: Strangeness suppression
A further point to discuss in our model is whether to use the sum or the difference of the number of strange and anti-strange quarks in Eq. (9) for the strangeness suppression. Fig. 18 illustrates the deviation between the two suppression models for the total reaction rate. For the system consisting initially only of light quarks, p +p, we see no sizeable differences between the sum and the difference of strange and anti-strange quarks in Eq. (9).
The system with an initial large difference between the number of strange and anti-strange quarks, Ω − +Ω + , converges to rather different equilibrium states for the two assumptions. The suppression with the sum leads to an overall larger total reaction rate and its equilibrium value is twice as large as the suppression with the difference assumption. However, both models produce rather similar results for times t < 50 fm, which is of relevance for the heavy-ion collisions considered in this work. Accordingly we use in PHSD the suppression with the sum of the number of strange and anti-strange quarks since both models give practically identical results in PHSD simulations of relativistic heavy-ion reactions.
Another issue relates to the actual value of the strangeness suppression factor λ which had been taken as λ = 0.5. In order to demonstrate the impact of the parameter λ on antibaryon spectra we show in Fig. 19 the rapidity distributions for central Pb+Pb collisions at 30 A GeV for λ=0.5 (dashed lines) and λ=1 (solid lines). Without strangeness suppression in 2 ↔ 3 reactions for hadrons with strange/antistrange quarks we find at 30A GeV that the rapidity spectra ofΛ +Σ 0 andΞ + are slightly shifted to lower values and broadened in comparison to the standard value of λ = 0.5. The spectrum for Ω − +Ω + very slightly broadens and thep spectrum is basically not influenced by the change of λ.