Particle Production via Strings and Baryon Stopping within a Hadronic Transport Approach

The stopping of baryons in heavy ion collisions at beam momenta of plab = 20 − 160A GeV is lacking a quantitative description within theoretical calculations. Heavy ion reactions at these energies are experimentally explored at the Super Proton Synchrotron (SPS) and the Relativistic Heavy Ion Collider (RHIC) and will be studied at future facilities such as FAIR and NICA. Since the net baryon density is determined by the amount of stopping, this is the pre-requisiste for any investigation of other observables related to structures in the QCD phase diagram such as a first-order phase transition or a critical endpoint. In this work we employ a string model for treating hadron-hadron interactions within a hadronic transport approach (SMASH, Simulating Many Accelerated Strongly-interacting Hadrons). Free parameters of the string excitation and decay are tuned to match experimental measurements in elementary proton-proton collisions. Afterwards, the model is applied to heavy ion collisions, where the experimentally observed change of the shape of the proton rapidity spectrum from a single peak structure to a double peak structure with increasing beam energy is reproduced. Heavy ion collisions provide the opportunity to study the formation process of string fragments in terms of formation times and reduced interaction cross-sections for pre-formed hadrons. A good agreement with the measured rapidity spectra of protons and pions is achieved while insights on the fragmentation process are obtained. In the future, the presented approach can be used to create event-by-event initial conditions for hybrid calculations.


Introduction
Understanding the properties of strongly interacting matter has been a long standing problem that can be addressed by studying the QCD phase diagram. In the case of high temperature and vanishing baryonic chemical potential, it was demonstrated that there can be a crossover phase transition instead of the first-order one, depending on the number of quark flavor and their masses [1]. More recent lattice QCD computations [2,3] show that there is a crossover transition between a hadronic gas and a quark-gluon plasma phase, if one goes to higher temperature (T ) keeping the net baryon chemical potential (µ B ) near zero. On the other hand, the QCD phase transition at T = 0 and finite-µ B has been studied based on effective models, such as the NJL model [4] and the composite-operator formalism [5]. It was demonstrated that the phase transition in this regime is first-order. The existence of a critical end point (CEP) is justified by the fact that the QCD matter exhibits different types of phase transition in two limiting cases [6]. Many theoretical and experimental studies in heavy ion physics have been aiming to find where the CEP is located in the T -µ B plane. On the experimental side, heavy ion collisions at various collision energies and several system sizes are carried out in order to probe a wide range in both temperature and baryon chemical potential. Those include the beam energy scan performed at RHIC [7,8,9] and the CERN-SPS [10,11,12]. In the future, this region will be studied further by CBM at FAIR and at NICA.
To connect the final state observables on particle yields and spectra with the properties of hot and dense QCD matter, sophisticated dynamical approaches are indespensable. The bulk observables in ultra-relativistic heavy ion collisions at RHIC and the LHC are successfully described by solving the hydrodynamic equations [13,14]. Hybrid approaches, which separate the non-equilibrium dynamics in the early stages from a hydrodynamic evolution of the thermalized medium, have proven to give a realistic description of heavy ion collisions also at lower beam energies [15,16]. The dynamical initialization of hybrid approaches is explored as well in [17,18,19].
In particular, the dynamics of baryon stopping has received some attention more recently, since it has been realized that the mean of the proton distribution should be understood before investigating higher order cumulants that are associated with a critical endpoint [7]. The experimentally observed stopping of baryons [12] is still lacking a quantitative description within theoretical models. In principle, there are three different options: (i) Push the gluon saturation picture down in energy and extend it to three dimensions as explored in [20] (ii) Study the source terms of projectile and target into the fireball fluid in a 3-fluid hydrodynamics approach [21,22] (iii) Investigate the details of a string hadron transport approach for the initial nonequilibrium evolution [23,24,25] Here, we follow the last point and apply the hadronic transport model SMASH to understand the stopping of baryons in the SPS energy range. This approach can be employed for the description of the early stages of a heavy ion collision, since microscopic transport is applicable to non-equilibrium circumstances. In this work, SMASH is employed to simulate the full evolution of a heavy ion collision. In the relevant energy range for this work, it is extremely important to have three-dimensional initial conditions for starting a hydrodynamic evolution, since the system cannot be assumed to be boost invariant and the colliding nuclei are too slow to reasonably neglect their longitudinal extend because of length contraction. The paper is structured as follows: Details of the transport model with a focus on the implementation of cross sections and the particle production at intermediate energies are given in section 2. Section 3 continues with calculations for proton-proton collisions, where the influence of varying the free parameters is investigated and the best possible set of parameters is determined. In section 4 we advance to heavy ion collisions, where one has the opportunity to study the interactions of string fragments and their formation process. Finally we provide calculations for the time of the collision where the colliding nuclei just passed through each other in section 5 which can serve as event-by-event initial state profiles for hydrodynamic calculations.

Model Description
In this work, we investigate baryon stopping within the transport model SMASH [26]. The code is publicly available on Github, see https://smash-transport.github.io/. The degrees of freedom within the calculation are hadrons. The properties of the hadrons are adopted from the Particle Data Group 2018 [27], where the more established resonances up to a mass of m ≈ 2 GeV are included. SMASH has been tested against an analytic solution of the Boltzmann equation [28] and strangeness as well as dilepton production has been confronted with experimental data at GSI-SIS energies [29,30].
The inelastic scatterings between hadrons at low energies are described via resonance formations and decays. Since there are no resonances with masses larger than m ≈ 2 GeV in the calculation, the cross section for resonance formations fades out when the center of mass energy of the interacting hadrons grows larger. This can be seen as an example for the proton-pion cross section shown in figure 1.
Let us start with an overview of the general setup of our approach while more details are presented in the following sub-sections. In order to investigate baryon stopping at higher incident energies, in this work, a string model is employed, where colliding hadrons are excited to strings which then fragment producing new particles. In the transition region between resonances and string processes, the respective cross-sections are weighted with a linear function to achieve a smooth interpolation between both regimes. This is important to avoid artificial high mass resonances that are suppressed in this way.
For most combinations of particle species, the transition region starts at the sum of the masses of the colliding hadrons plus 0.9 GeV and has a width of 1 GeV. For two very important special cases, the transition region is individually specified. The first one is nucleon-pion collisions, where the transition from resonances to strings takes place between 1.9 GeV and 2.2 GeV. In this case, the transition region is shorter, because the contribution from resonances is too small above √ s = 2.2 GeV. The second special case is in collisions between two nucleons. Here, the transition region spans from 4 GeV to 5 GeV. Compared to the default, the transition region is shifted to higher √ s because up to 4 GeV, the total cross section from resonances reproduces the measurement and the exclusive cross sections from resonances are more realistic at low energies.
The calculation for the string excitation is split into hard and soft processes. The hard string processes are relevant for very high energetic binary interactions as can be seen in figure 1, where perturbative QCD is applicable. For the description of the pQCD scatterings, the string excitation and the string fragmentation, Pythia 8.235 [31,32] is used. The hard string routine is described in more detail in section 2.2. The hard string process, where pQCD interactions are involved, is based on the p T -ordered multiparton interaction (MPI) framework with initial and final state radiations [33]. Given that pQCD is not applicable at low momentum scale, the lower p T threshold of those partonic interactions is chosen to be 1.5 GeV and the pQCD cross section is computed accordingly.
In the transition region where the energy is too large to have cross sections via resonances but too low to apply pQCD, a phenomenological model for the excitation of strings is implemented. In single diffractive, double diffractive and non-diffractive processes, strings are excited in hadronic interactions. Using the calculated mass and momentum of the string as well as the flavor of the leading quarks, Pythia is employed only for the fragmentation of the string. Details of the string excitation at intermediate energies can be found in section 2.3. Figure 1 also shows the contribution of elastic scatterings to the total cross section. Elastic collisions play an important role at all beam energies, since in heavy ion collisions, a large amount of nucleons will scatter elastically. Especially for baryon stopping, the angular distribution of the final state particles in elastic scatterings plays an important role as shown in section 4. Elastic Resonance Soft string Hard string Total Figure 1. Cross section of a proton interacting with a negatively charged pion as a function of the center of mass energy of the colliding hadrons within SMASH. The total cross section is split into contributions from elastic collisions, resonance formations, soft string excitations and hard string excitations via Pythia.

Cross Sections for String Processes
The total cross section for each pair of hadrons is required within SMASH to determine, if two particles will scatter. Subsequently, the actual process is decided randomly according to the underlying partial cross-sections for the different channels. The total σ tot and the elastic σ el cross section, which is shown for the example of p-π − collisions in figure 1, is parameterized to fit the experimental data. The inelastic cross section σ inel is the difference between the two The parametrizations for the total and elastic cross sections are taken from [34] and [35,36] respectively. For the process selection, cross sections for both the single σ SD and double diffractive σ DD processes are necessary. They are estimated in [37] and implemented in Pythia. From σ SD and σ DD the non-diffractive cross section σ ND is derived The non-diffractive cross section includes the hard and soft non-diffractive processes. The cross section for hard non-diffractive processes is based on the pQCD cross section σ hard from partonic interactions. It is given by where σ k i,j |p T,min is the cross section for a subprocess k between two partonic flavors i and j with minimum transverse momentum transferp T,min , which is chosen to be 1.5 GeV. The parton distribution function f i (x) provides the average number of flavor i carrying the momentum fraction x of the incoming hadron. The NNPDF 2.3 parton distribution function with QED correction [38] is used in this work. The sum takes each possible combination of partons from each ingoing hadron into account.
This pQCD cross section can therefore be larger than σ ND , incorporating the information of multi-parton scattering. We take the multi-partion interaction (MPI) picture [39] and interpret the ratio σ hard /σ ND to be the number of partonic interactions involved in a hadronic interaction, rather than the probability to have a hard nondiffractive interaction. In addition, the number of parton interactions is assumed to follow a Poissonian distribution, where the average is given by σ hard /σ ND . The probability of having no hard interaction is calculated according to the Poissonian distribution as In this case, the process is assumed to be soft non-diffractive, leading to a soft nondiffractive cross section σ ND,soft of Finally, the cross section σ ND,hard for the hard string process follows as Since the pQCD cross sections can only be applied at sufficiently large energies, there is no contribution from hard non-diffractive processes below collision energies of 10 GeV. If the energy is smaller, all non-diffractive processes will be soft. Due to the fact that Pythia 8 accepts only (anti-)nucleons and pions as the incoming hadrons, it is necessary to extrapolate these processes to handle arbitrary pairs of incoming hadrons. This is done by mapping hadronic species onto pions and nucleons and then rescaling cross sections based on the additive quark model. If a baryon has positive electric charge, it is mapped onto a proton. Otherwise it is mapped onto a neutron. Similarly, if a meson has positive or vanishing electric charge, it is mapped onto π + . Otherwise it is mapped onto π − . The additive quark model [40] is implemented in a similar manner as in the UrQMD model [41,42]. The total, elastic, diffractive and pQCD cross sections are multiplied by a constant factor, which depends on the valence quark/antiquark contents of the incoming hadrons.
where n q/q is the number of quark/antiquark constituents, while n s/s is the number of strange quark/antiquark constituents. h and h stand for respectively the incoming and mapped hadronic species.

Hard String Routine
Hard non-diffractive string excitations dominate the hadronic cross section at large center of mass energies. As mentioned in section 2.1, Pythia 8 accepts only a limited number of species as incoming hadrons, and it is necessary to extrapolate hard nondiffractive scattering handled by Pythia 8 to all other hadronic species. This is particularly crucial in high-energy heavy ion collisions, where plenty of hadrons other than (anti-)nucleons and pions are produced by primary nucleon-nucleon collisions. To do this extrapolation, we rely on the assumption that the structure functions (or parton distribution functions) of all mesons and (anti-)baryons look similar to respectively those of pions and (anti-)nucleons once we swap the valence quark flavors. Technically, this is achieved by mapping different hadron species to (anti-)nucleons and pions, where the quantum numbers of the original and mapped particle are as similar as possible. This is done in the same way as the mapping for the cross sections, which is described in section 2.1. Before the produced strings are fragmented within the Pythia calculation, light (anti-)quarks are exchanged with quarks of the original flavor. The momenta of all particles are rescaled in order to conserve the energy of the system, since the constituent masses are affected by the flavor exchange. Due to annihilation processes, it is not always possible to find a quark with the flavor of the mapped quark. In this case, a gluon is split into a quark-antiquark pair with the flavor of the mapped quark or anti-quark.

Soft String Routine
The soft string excitations are the most abundant processes in the intermediate energy range as can be seen in figure 1. As in UrQMD [41,42], the excitation of a soft string can be performed according to one of three subprocesses, the single diffractive (see section 2.3.1), the double diffractive (see section 2.3.2) and the most common subprocess the non-diffractive case (see section 2.3.3). All soft string processes rely on Pythia for the fragmentation of the strings into final-state hadrons.

Single Diffractive
The single diffractive process describes the interaction between two hadrons, where exactly one of the two colliding hadrons A and B is excited to form a string X The excited string X has a larger mass than the incoming hadron. The differential cross section, as a function of the string mass M X from diffractive excitation, is given by [43] dσ SD dM 2 Once the string mass is sampled, the three-momenum p CM of the string in the centerof-mass frame can be evaluated by solving where m H is the mass of the incoming hadron, which remains intact. Following the UrQMD approach [41], the transverse momentum transfer p ⊥ between the incoming hadrons is assumed to follow a Gaussian distribution where σ T is a free parameter that is constrained by observables in proton-proton collisions in section 3.3. To completely determine kinematics of the string-hadron system, we sample the transverse momentum transfer p ⊥ with a maximum of p ⊥,max = p CM . The string has a longitudinal momentum p = (p 2 CM − p 2 ⊥ ) 1/2 , which is parallel to the collision axis. Knowing the mass and the momenta of the reaction products, one can calculate the velocity of the string in order to boost into its rest frame, where the fragmentation machinery from Pythia is employed to obtain the particles in the final state of the interaction.

Double Diffractive
The double diffractive subprocess is a collision in which the two incoming hadrons A and B are both excited to strings The dynamics of the interaction is determined in the center of mass frame of the incoming hadrons, where the collision axis is set to be the longitudinal direction.
Kinematics of the double-diffractive excitation is modeled via pomeron exchange between gluons from two incoming hadrons. Those gluons exchange transverse and lightcone momenta, such that they remain on-shell after the momentum exchange [41,42]. The light cone momentum fraction x of each gluon is sampled from the parton distribution function for gluons, which is assumed to be of the form where β is set to be 0.5. The light cone momenta p ± of the hadrons are given by where E is the energy of the hadron and p is the projection of the momentum on the collision axis of the colliding hadrons. The light cone momenta of the exchanged gluons are calculated as p ± g = x ± p ± . The distribution for the transverse momentum transfer p ⊥ between gluons is taken to be Gaussian, whose width is the same as in the single-diffractive case (see equation (11)). The lightcone momenta q − A and q + B , which are gained by the gluons from hadrons A and B, are given by Note that the collision axis is defined as the direction in which the hadron A is moving. The lightcone momenutum Q ± transferred from the hadron B to A is given by and it leads to evaluation of the energy ∆E and longitudinal momentum ∆p transferred from B to A as The mass of both excited strings can be calculated individually using the energymomentum relation. Each string is then fragmented in the rest frame of that string using the implementation of the fragmentation within Pythia.

Non-Diffractive
The non-diffractive string excitation is the most probable soft process, and therefore has the largest impact on the dynamics of the produced particles in the SPS energy region. During the interaction, each hadron emits one valence quark, which is adopted by the other hadron. The exchanged valence quark carries a fraction of the longitudinal momentum of the hadron it is emitted from. The light cone momentum fraction carried by the exchanged quark is sampled according to the parton distribution function for quarks, which is assumed to have the following functional form: where α and β are in general free parameters. In section 3.1, the dependence of the particle production on the PDF is studied in detail and the parameters are adjusted such that the measured dynamics is reproduced as well as possible while supporting the physical picture of a valence quark exchange. The momentum transfer in the transverse direction is sampled according to the same Gaussian as in the single diffractive and double diffractive case, using equation (11). With the light cone momentum fraction each exchanged quark carries and the transferred transverse momentum, the light cone momentum transfer is written as [41,42] where x A and x B are the light cone momentum fractions for the exchanged quarks, p ± A and p ± B are the light cone momenta of the colliding hadrons before the collision and p T is the transferred transverse momentum. The exchanged energy and longitudinal momentum can be calculated using equation (19). The masses of the strings are obtained using the relativistic energy-momentum relation and each string is fragmented individually in the rest frame of the string using Pythia.

String Fragmentation
Once the mass of the excited string and the flavor of the quarks spanning the string is determined, the string is fragmented into hadrons by employing Pythia. Within Pythia, the species of the fragmented hadron follows from the flavour of the quarkantiquark or diquark-antidiquark pair that is produced. While the light quarks have the same probability to be produced, there are empirical suppression factors for producing heavier quarks and diquarks.
The transverse momentum of each string fragment is sampled from a Gaussian distribution with a width of σ T,string which is a free parameter that is tuned to experimental data in section 3.3. The longitudinal momentum of each string fragment is determined using the fragmentation function. Pythia is based on the symmetric Lund fragmentation function [31], which has the following shape m T is the transverse mass of the string fragment while a and b are free model parameters.
For the fragmentation of leading baryons produced in soft non-diffractive processes, the parameters a and b are chosen differently from Pythia. The consequences of this treatment are discussed in more detail in section 3.2.

Particle Formation
A string fragments into hadrons by producing quark-antiquark pairs. In a dynamical picture, the pair production does not happen simultaneously but at different points in time. Figure 2 illustrates how a string fragments in coordinate space within the yoyo model. The straight lines indicate the trajectories of (anti)quarks or (anti)diquarks.  In practice, all particles in SMASH are immediately produced once the colliding hadrons reach their point of closest approach. Until the formation time has passed, the cross section of the string fragments are multiplied by a cross section scaling factor f σ . For most string fragments, this factor is initially 0. However, since the leading string fragments contain quarks that do not originate from a pair production but from the initially colliding hadrons, the initial cross section scaling factor is not zero for leading string fragments. The initial cross section scaling factor for each string fragment is set to be the number of quarks from the initially colliding hadrons contained in the fragment divided by the total number of quarks of that fragment. For example, a leading baryon that contains a diquark from the initially colliding hadrons is assigned a scaling factor of f σ = 2/3 and a meson at the other end of the string that contains another quark from the initially colliding hadrons is assigned a scaling factor of f σ = 1/2.
Instead of having a constant cross section scaling factor until the time of formation, where the particle suddenly is allowed to interact, it is possible to mimic a continuous formation process by increasing the cross section scaling factor with time. To realize a continuous formation, the cross section scaling factor becomes a function of time f σ = f σ (t). This function needs to have the initially assigned scaling factor f 0 as described above at the time t prod when the particle is produced and f σ (t form ) = 1 at the formation time t form . Between the two points, the cross section scaling factor grows with a given power α in order to have a simple but flexible functional shape. Using the three conditions, the function f σ (t) is written as follows This function is only used for t prod < t < t form , since it has no meaning before the particle is produced and the scaling factor is f σ (t) = 1 for t > t form , when the particle is fully formed. The cross section scaling factor as a function of time for different values Cross section scaling factor f σ as a function of time for different powers α with which the cross section grows in time. In this example, the initial cross section scaling factor is set to be f σ (t prod ) = 0.
of α is shown in figure 3. The initial cross section scaling factor is set to f 0 = 0 in this figure. In the limit of α going to infinity, one recovers a step function, while for small positive values of α, the particles form immediately. In section 4, the effect of the details of the particle formation on particle spectra in heavy ion collisions is investigated.

Elastic Collisions
Elastic collisions play an important role in describing heavy ion collisions, since the contribution of elastic collisions to the total hadronic cross sections is at all energies significant as can be seen for example for proton-π − collisions in figure 1. In SMASH, all hadrons have a finite cross section to interact elastically. The most important elastic cross sections are parametrizations of the experimental data, if available. If the elastic cross section is not measured for a pair of particle species, the additive quark model is applied to obtain a cross section for that pair of particles [41] as shown in equation (7).
While the angular distributions are close to being isotropic at low collision energies, they are more forward-backward peaked at larger collision energies [44]. Especially going to higher collision energies, including anisotropic angular distributions for elastic collisions is, therefore, of major importance for describing the dynamics of heavy ion collisions. This is shown in a comparison between calculations with isotropic and anisotropic angular distributions provided in section 4.
For the elastic collisions of nucleons at relatively low collision energies, parametrizations for the angular distributions have been proposed in [45]. Since there is few experimental data for elastic collisions of hadronic pairs other than nucleonnucleon ones, we implement the Cugnon parametrizations [45] for nucleon-nucleon elastic collisions and extrapolate them to all other hadronic pairs. Given that the Cugnon parametrization is a function of p lab of fixed-target experiments, there is an ambiguity for incoming particles with different masses. To circumvent this, we first compute the center-of-mass momentum p cm of the collision. Then a new lab-frame momentum p * lab is evaluated from a new mandelstam variable s * , which yields the original center-of-mass momentum when the nucleon mass is assumed Lastly, the differential cross section of elastic proton-proton collisions is extrapolated to all hadronic pairs such that where θ is the scattering angle. In addition, these angular distributions are extrapolated to arbitrary collision energies in order to obtain forward-backward peaked angular distributions for elastic collisions at large collision energies.

Proton-Proton Collisions
In a heavy ion collision, the stopping of baryons is mostly determined by the first interactions of the participants. Therefore, experimental data for elementary protonproton collisions is used to adjust the parameters of the string approach. Even though most of the parameters influence multiple observables, this section introduces the most important parameters and demonstrates their effect on the particle production of (anti)protons, pions and kaons to give some insight on how the value of each parameter is chosen. Within this section, each parameter is varied separately, while all others are kept constant. If not further specified, the value of each parameter can be found in table 1. We mainly concentrate on the highest SPS energy, since the experimental data set is the largest at that energy and in the end of the section show results for all other beam energies as well. Popcorn rate Probability for popcorn processes as described in section 3.6 0.15 Table 1. Default set of parameters of SMASH-1.6 tuned to reproduce the experimental data for proton-proton collisions at SPS energies with brief description.

Parton Distribution Function
In the SPS energy range, the soft non-diffractive string processes are the dominant interactions in proton-proton collisions. As described in section 2.3.3, the amount of exchanged longitudinal momentum is determined by the momentum fraction x the exchanged valence quark carries. This does not only affect the dynamics of the string before the fragmentation, but also the mass of the string and therefore the energy available for producing string fragments. The PDFs used for the calculations are shown in figure 4 (left). Figure 4 (right) shows the dependence of the longitudinal momentum of protons in proton-proton collisions on the value of the parameter β of the parton distribution function as defined in equation (20). The longitudinal momentum distribution of protons is clearly the most relevant quantity to understand baryon stopping in heavy-ion collisions. The larger momentum transfer in the longitudinal direction is reflected in the x F = p z /p z,beam distribution of protons as they are shifted to higher x F for larger values of β quark . With a softer PDF, there is more energy available for the production of new particles. This reflects in the proton yield at low x F , which increases by a factor of 2 when using a very soft PDF. Even though the proton x F distribution is not reproduced quantitatively, a reasonable agreement is found for β quark = 7.
Note that the description of the longitudinal momentum of protons in proton-proton collisions proves to be very challenging within theoretical models [47]. Modifications to improve the agreement between a different model and the data have been suggested [48]. One of the suggestions from [48], is to modify how a proton is split into a quark and a diquark. An option to specify the probability to split a proton into uu + d rather than ud + u is implemented in SMASH but does not improve the overall agreement with the measurement. Setting β quark = 7 and α quark = 2, the mean value of the PDF is 2/9, which is close to the expectation of 1/3 which is assuming that there are three valence quarks sharing the full momentum of the proton.

Fragmentation Function
The longitudinal momentum of each string fragment is determined by the shape of the fragmentation function. Starting at the forward and backward ends of the string, the fraction of the remaining lightcone momentum is sampled from the fragmentation function. Pythia employs the symmetric Lund fragmentation function defined in equation (23). Figure 5 (left) shows the Lund fragmentation function for two different values of the parameter b. On the right panel, the distribution of x F = for protons is plotted for two different settings. The curve labeled Lund fragmentation refers to a calculation where the softer fragmentation from the left part of the figure is used consistently throughout the fragmentation. Comparing to the experimental data shows that protons obtain too little longitudinal momentum that way. Therefore, the protons require a different fragmentation mechanism. To increase the longitudinal momentum of protons, without producing harder light mesons than before, the harder fragmentation function shown in the left panel is used for leading baryons in soft non-diffractive string processes. The other curve on the right panel of figure 5 shows the result of the calculation after that adjustment. A drastic improvement of the agreement with the experimental data is observed.
Even though the fragmentation function for non-leading hadrons is considered as an intrinsic property of a string which does not depend on what happens outside, the leading diquark holds information on the initial state kinematics. That legitimates having a separate fragmentation function to determine lightcone momenta of leading baryons.   The influence of using a separate fragmentation function for leading baryons on the transverse momentum is shown in figure 6, where the mean transverse momentum as a function of x F is shown for protons and pions. The shape of the transverse momentum of protons as a function of the longitudinal momentum fraction is better reproduced with the unmodified Lund fragmentation function. The curve employing a separate fragmentation function reflects the expected change that has been observed in figure  5. To understand the stopping in heavy-ion collisions, the match of the transverse momentum at midrapidity corresponding to low values of x F is most important for our present work. Please refer to section 3.3 for a more detailed discussion of transverse momenta.
Let us now demonstrate in detail how the parameters for the fragmentation function for leading baryons from soft non-diffractive processes (a leading , b leading ) and for all other π + NA49 data SMASH 1.6 SMASH Lund fragmentation 27 GeV compared to experimental data [46,49]. The curve labeled Lund fragmentation is obtained using the standard Pythia fragmentation mechanism employing the same fragmentation function for each string fragment, while for the other curve, a separate harder fragmentation function is used for leading baryons from soft non-diffractive string processes.
particles (a string , and b string ) have been determined. Figure 7 shows the longitudinal momentum distribution for protons and pions in pp collisions at the highest SPS energy for different values of b leading . In general higher values of b leading are prefered by the proton x F distribution, but there needs to be enough energy for particle production as well. Therefore, b leading = 2.0 GeV −2 provides the best compromise to generate hard enough protons, while still producing a reasonable amount of pions. In addition, higher values of b leading lead to a double-peal structure in the x F distribution that is not supported by the experimental data. x F distribution of protons is observed in the height of the bump at x F ≈ 0.3. For larger a leading , the peak is more pronounced, since small values of a leading correspond to harder fragmentation functions. The similarity between the curves with the hardest fragmentation function is caused by the fact, that the string fragmentation fails numerically, if there is not enough energy left to produce new particles. Therefore, very high momentum fractions are rejected more often and therefore, the difference in the fragmenattion function is not visible anymore in the observable. In order to achieve a good agreement with data for protons without requiring many attempts to determine the kinematics, the value of a leading is set to 0.2.
The fragmentation function that is used for all other particles has a strong effect on the production of light mesons. The value of the parameter b, which will be referred to as b string in the following, is varied in figure 9. The rapidity spectra of positively and negatively charged pions are sensitive to small changes in the parameter b string . A softer fragmentation function will lead to more low-energetic pions, while a harder fragmentation function produces pions with larger momenta.
While the mid-rapidity yield of positively charged pions is overestimated for all three values of b string , the production of negatively charged pions is well described. At lower beam energies, the pion multiplicity is slightly lower compared to the data as can be seen in figure 18. Therefore, the overall best agreement is obtained with b string = 0.55 GeV −2 .
The final parameter a of the fragmentation function used for particles that are not leading baryons is called a string in the following. It is varied in figure 10 which shows the x F distribution of protons and positively charged pions. Higher a string corresponds to a harder fragmentation function which is reflected in the soft proton sector and the pion x F distribution. Overall, the best agreement, considering both distributions, is found with a string = 2.0.
For completeness, figure 11 shows the fragmentation function used for all string fragments apart from leading baryons in soft non-diffractive processes. For the other string fragments, softer fragmentation functions are applied, where the difference between the particle species originates exclusively from the m T dependence in equation (23).

Transverse Momentum Production
Transverse momentum is produced in two steps in proton-proton collisions: first in the excitation process adjusted by σ T and afterwards during the fragmentation tuned by changing σ T,string . The initial transverse momentum transfer between the interacting hadrons is sampled according to a Gaussian with a width of σ T as described in section 2.3. Figure 12 shows the dependence of the mean transverse momentum of protons and pions as a function of x F on the value of σ T . 27 GeV compared to experimental data [49,46], note that the zero is suppressed to zoom in to the region of interest.
The mean transverse momentum of pions shows only a weak dependence on the value of σ T . At low x F , the protons are fragmented from a string. Therefore they show the same behavior as pions. Most large x F protons in proton-proton collisions are however not fragmented from a string, but only took part in a singe diffractive process. Their transverse momentum is directly sampled from the Gaussian with a width of σ T , which explains the strong dependence on σ T . σ T = 0.42 GeV is therefore fixed to match the mean transverse momentum at large x F .
The production of transverse momentum during the fragmentation of a string is regulated by σ T,string . In Pythia, a Gaussian with a width of σ T,string is used to sample the transverse momentum of each individual string fragment. The influence of changing σ T,string on the mean transverse momentum of pions and protons is shown in figure 13. The transverse momentum distribution for the fragmentation is much more important for particle production than the one in the string excitation process. 27 GeV compared to experimental data [49,46].
In the case of pions, the transverse momentum is scaled up for all bins of x F . For protons, the opposite behavior to what is seen when varying σ T is observed. At small x F , where the protons originate from a string fragmentation, the proton p T is strongly dependent on σ T,string , while at x F ≈ 1 all curves lie on top of each other. The value of σ T,string therefore needs to be tuned to multiple particle species simultaneously. The best agreement is found for σ T,string = 0.5 GeV, where for protons too little transverse momentum is produced, while the pions at low x F obtain too much p T .
Since the NA61 collaboration recently measured the transverse mass of protons at mid-rapidity as a function of collision energy, calculations from SMASH are compared with the data and other transport approaches in figure 14. Please note that the rapidity ranges do not match exactly, since the HSD and UrQMD calculations were performed before the experiment was carried out, therefore the comparison is not fully quantitative. The UrQMD calculation overshoots the data at low √ s due to the transition from resonances to strings which is located at higher energies for the binary collisions than in HSD and SMASH. The shape of the HSD curve and the SMASH calculation qualitatively follow the trend of experimental data, while both underpredict the mean transverse momentum slightly.

Strangeness Production
The production of strange quarks heavily relies on the probability of producing an ss pair compared to the probability of producing a light qq pair during the string fragmentation.  Figure 14. Difference between mean transverse mass of protons at mid-rapidity and the proton mass in proton-proton collisions as a function of the center of mass energy √ s of the collision compared to experimental data and other hadronic transport calculations [51,42,41,52] To suppress the production of strange quark pairs according to their higher mass, the strangeness suppression factor λ s is introduced: where P (uū), P (dd) and P (ss) denote the probabilities to produce a up-antiup, downantidown and a strange-antistrange quark pair respectively. The impact of varying this parameter on the kaon rapidity spectra is shown in figure 15. Without affecting  the dynamics of the system much, the strangeness suppression factor regulates the multiplicity of strange hadrons. Since the rapidity distribution in our calculation is slightly steeper than the measured one, the strangeness suppression factor is set to λ s = 0.16 in order to obtain a good agreement for the total kaon multiplicity. For tuning λ s , only the positively charged kaons are considered, since the energy dependence of the negatively charged kaons is not as well understood as can be seen in the bottom right panel of figure 18.

Diquark Production
Similar to the description for the production of strange quarks, a diquark suppression factor λ diquark is introduced to quantify the likelihood of producing diquarks: A diquark and an antidiquark always combine to a baryon and an antibaryon, since a meson cannot contain two (anti)quarks but only one quark and one antiquark. Since diquarks are present in a much larger fraction than the newly produced pairs as valence quarks in the excited baryons, the antiproton production constrains the diquark suppression factor much more directly. The comparison of the rapidity spectrum of antiprotons for different values of λ diquark is shown in figure 16. The antiproton multiplicity is very low, resulting in a small diquark suppression factor and a value of λ diquark = 0.036 yields the best agreement with the measured antiproton rapidity spectrum. Even though the data point at mid-rapidity suggests a larger λ diquark , all other points are reproduced very well and at lower energies a higher antiproton production contradicts the measurement as shown in figure 18.

Popcorn Rate
When a diquark-antidiquark pair is produced, they will recombine with surrounding quarks and antiquarks, forming new baryons. Since the diquark and the antidiquark are produced in a pair production, they are connected via their color charge. This will in many cases lead to the two fragmented baryons to be produced next to each other in phase space. It is however also possible to create another quark-antiquark pair in the color field spanned by the diquark and the antidiquark. This leads to the production of a meson between the two baryons [53]. In the case of a baryonic string, it is, within a popcorn process, possible to fragment a meson at the diquark end of the string. Given that a diquark-antidiquark pair is created, the probability of such a process is given by the popcorn rate, which is a Pythia parameter that can be varied in order to reproduce the experimental data. The effect of changing the popcorn rate on the dynamics of protons is investigated in figure 17. Increasing the popcorn rate leads to protons being shifted to low x F . In the transverse direction, the proton p T increases with a growing popcorn rate. While a large popcorn rate results in a better agreement for the transverse momentum, the shape of the x F distribution is not compatible with the data. Because the x F distribution is, like the data, flat in the low x F region in the experimental data and a fair agreement in the transverse momentum can be obtained, the popcorn rate is set to 0.15. Compared to the effect on the proton dynamics, the other particle species are only slightly affected by changing the popcorn rate.

Proton+Proton Results Overview
In this Section, we present the full set of final results for proton-proton collisions including mean transverse momentum as a function of x F and rapidity spectra for protons, antiprotons, positively and negatively charged pions and kaons. This serves to benchmark the whole calculation in elementary collisions and provides the baseline for heavy ion calculations. Figure 18 shows the rapidity distribution for the mentioned particle species for different collision energies. A good agreement over the entire SPS energy range is achieved for pions, positively charged kaons and antiprotons. The energy dependence of the negatively charged kaon multiplicity is too strong inside the string model compared to the measurement. The rapidity spectrum is therefore only well  Figure 18. Rapidity spectra of (anti)protons, positively and negatively charged pions and kaons at different collision energies compared to experimental data [50].
reproduced at p lab = 80 GeV, while the calculation overshoots at larger and undershoots at lower collision energies. The proton rapidity spectrum follows the shape of the data roughly and the quantitative agreement is reasonable, but the longitudinal momentum of protons in proton-proton collisions proves to be quite challenging to describe within the string model as we expected regarding the difficulties to obtain a reasonable x F distribution. Figure 19 shows the mean transverse momentum of the different particle species as a function of x F . The mean transverse momentum of protons at low and intermediate K − s = 17.27 GeV NA49 data SMASH 1.6 Figure 19. Mean transverse momentum of (anti)protons, positively and negatively charged pions and kaons in proton-proton collisions at √ s = 17.27 GeV compared to experimental data [46,49,54].
x F is underestimated in the calculation while the data is reproduced at large x F . The opposite behavior is observed for the other particle species, where the mean transverse momentum at low x F is slightly overestimated while the pion mean transverse momentum undershoots the data at large x F . The mean transverse momentum for all particle species does not deviate much from the data so that overall a sufficient agreement with the measurement is found in proton-proton collisions for advancing to heavy ion collisions.

Heavy Ion Collisions
After adjusting the parameters for particle production in string processes to experimental observations in proton-proton collisions, we present calculations for heavy ion collisions. The evolution of the shape of the proton rapidity spectrum from a single peak structure at low collision energies to a double peak structure at high energies is observed in central collisions indicative of the dynamics of baryon stopping. In addition to understanding the net baryon content in the fireball, it is possible to gain insight on the formation process of string fragments.
Since the first collisions mainly take place via string excitation and fragmentation in the considered energy regime, let us first discuss the impact of the formation times and cross-section scaling factors on the results. During the fragmentation process, the particles are not immediately fully interacting hadrons but rather some pre-formed states that interact with lower cross-sections (see description in section 2.5).
The influence of the formation time on the particle spectra is first studied in the simplest case, where the cross section scaling factor f σ , as introduced in section 2.5, is a step function in time and does not increase continuously. For most string fragments, this implies that they instantly form, when the formation time has expired.   While all three calculations reproduce the shape of the measured rapidity spectrum, all curves fail to describe the number of stopped protons at mid-rapidity. This reflects the fact that the formation times are too large for the string fragments to form while the nuclei still overlap. If the cross sections continuously grow with time, there is a small probability for string fragments to immediately interact. As shown in figure 20 (right), this enhances the amount of stopped protons significantly. Figure 20 (right) shows the calculation for a fixed formation time of τ form = 1 fm for different powers in which the cross section scaling factor grows in time. Using a step function (α = −1) gives similar results as the quadratic increase. When the cross section scaling factor grows with the square root in time, the string fragments interact too much at early times and the protons are stopped far too much. Only for the linearly growing cross section scaling factor, the amount of stopping can be reproduced.
A deeper understanding about how the power α translates to more stopping can be gained by studying the interaction rate as a function of time for the different scenarios. Figure 21 shows the rate of different interactions as a function of time.
The rates are compared between a calculation with linearly growing cross section scaling factor (α = 1) and a cross section scaling factor, that does not continuously grow in time (α = −1). The interaction rate at early times is dominated by string processes. At later times, the energy is not sufficiently large to excite strings and the strongest contribution to the interaction rate stems from resonance formations and decays. For the stopping, the most interesting period is right after the initial collisions. There one can see that if the cross section immediately starts growing, the rate of elastic collisions and resonance formations is significantly increased. These additional interactions are responsible for the higher proton yield at mid-rapidity. At √ s N N = 17.27 GeV, the formation time affects the results only slightly even when the cross sections grow linearly in time. Going to lower collision energies, changing the formation time directly reflects in the rapidity spectra as shown in figure 22. In the case of slower, and less Lorentz contracted, nuclei the passing time is on the order or larger than the formation time. With shorter formation times, the cross section scaling factors grow faster in time which leads to more stopping of protons and larger pion multiplicities at mid-rapidity. A formation time of τ = 1 fm is best suited for reproducing the proton and pion rapidity spectrum. Figure 23 shows the rapidity spectrum of protons in central heavy ion collisions compared to experimental data and UrQMD calculations for different collision energies for the final set of parameters for particle formation (τ = 1 fm and α = 1). Over the entire SPS energy range, a good agreement between the SMASH calculation and the experimental data is found. Even though the proton multiplicity at low SPS energies is overestimated in the SMASH results, the evolution of the shape of the proton rapidity spectrum from a single peak at low energies to a double peak structure at large collision energies is well reproduced. At low beam energies, a small fraction of protons is bound in light nuclei and should not be counted in the proton spectra. The clustering is not taken into account in the shown SMASH results, which might be part of the reason for the overshoot of protons at low beam energies. To put our results into context, we compare to UrQMD calculations, where a very similar treatment of string processes is applied. In general, the protons within the UrQMD calculation are stopped more at mid-rapidity than the protons in the SMASH calculation. In UrQMD, the cross section of an unformed particle is kept constant until the formation time of that particle is reached. This corresponds to the SMASH calculation with α = −1 shown in figure 20 (right), where the least stopping is observed in the case of α = −1, so the details of the   Figure 23. Rapidity spectra of protons in central lead-lead collisions at different beam energies compared to experimental data [55,12] and UrQMD calculations [57] formation process of string fragments is not the main source of the difference between SMASH and UrQMD. An important ingredient for understanding the shape of the proton rapidity spectrum are the anisotropic angular distributions for elastic collisions between all hadrons as described in section 2.6. Figure 24 shows the rapidity spectrum of netprotons in a calculation of lead-lead collisions at √ s N N = 8.765GeV at different times comparing isotropic and anisotropic angular distributions in elastic collisions. At t = 2 fm, the nuclei are still in the process of passing through each other, explaining why a large portion of protons are still located at beam rapidity. Up to this point, the dynamics of the system are dominated by primary interactions between nucleons. A significant difference is already at that time observed between the calculations with isotropic and anisotropic elastic scatterings. Advancing in time, the net-proton number increases, mostly due to resonance decays. The difference between the calculations with isotropic and anisotropic elastic collisions is not washed out during the evolution of the system but can be observed even after all resonances have decayed. A double peak structure only builds up, when the anisotropy of elastic scatterings is properly taken into account.
To conclude the study of particle production in heavy ion collisions, figure 25 shows the rapidity spectrum of negatively charged pions for different collision energies. As shown in figure 22 (right), the pion production is relatively well understood at intermediate SPS energies. Similar to the intermediate energies, a good agreement with the data can be observed at the lowest collision energies in the SPS range. At top SPS energies, the multiplicity of negatively charged pions is underestimated but still a reasonable agreement is found. In comparison to the SMASH results, pions are more   Figure 25. Rapidity spectrum of negatively charged pions in central lead-lead collisions at different beam energies compared to experimental data [56,58] and UrQMD calculations [57].
abundantly produced in the UrQMD calculations. Compared to the data, UrQMD describes the pion production very well at high energies while SMASH gives a better description for low collision energies.

Initial State Calculations
In this last Section, we would like to show how the results from our approach can in the future be employed for initial conditions for hydroydnamic calculations. This has been very successfully done in a hybrid approach based on UrQMD initial conditions [15,59]. More recently, a toy model for initial conditions including a 3D Glauber model with energy loss according to a string picture is developed [18]. Further, dynamical initial states are constructed based on UrQMD [17]. In a similar fashion, the hadronic transport approach JAM is combined with relativistic viscous hydrodynamics [19]. The advantage of a dynamical approach is to include event-by-event fluctuations of all relevant quantities and having full 3D distributions for all quantum numbers available. Figure 26 shows the energy density and the net baryon density in a single event for a slice at z = 0 at the time when the two colliding nuclei have just passed through each other. Due to secondary interactions and the produced transverse momentum, some strings are not aligned with the beam axis, which reflects in small lines of large energy density. Since the baryon density of the string is located at the end, this structure cannot be observed on the right panel of figure 26. The scale on the left hand side of figure 26 ranges up to large energy densities, well in the regime where a quark-gluon plasma should be formed. Therefore, a description of the dynamical evolution in the hot and dense system in terms of hydrodynamics seems more appropriate than a pure hadronic transport approach.
Advancing to the longitudinal dynamics of the system, figure 27 shows on the left panel the spacetime rapidity distribution of net baryons compared to the momentum space distribution on the right hand side in heavy ion collisions at different collision energies. Again, the distributions are plotted at the time, when the nuclei have just passed through each other. The momentum space distribution at √ s N N = 6.27 GeV shows a peak at mid-rapidity, while, with increasing energy, a flat plateau is developed. Even though the spacetime rapidity spectrum shows a similar behavior with increasing energy, momentum space and coordinate space rapidity spectra differ drastically. This supports the finding, that the Bjorken assumption breaks down at lower beam energies and a full 3-dimensional initial state is more realistic.

Summary
Baryon stopping in the SPS energy range is studied within the hadronic transport model SMASH. Going to high collision energies, string excitations and fragmentations are the most important processes, since the contribution from resonance excitations fades out. The string model inroduced in this work is split into soft and hard processes, where the soft processes dominate the cross section at intermediate energies while the hard processes are most important at very high energetic interactions. The soft string processes follow the UrQMD approach while hard processes are handled via Pythia.
To take the dynamics of particle production in a string model into account, a formation time is introduced during which the cross sections of string fragments are reduced. In order to mimic a continuous particle formation process, a mechanism is introduced to smoothly increase the cross section of forming particles over time. The model has been benchmarked against experimental data from NA49 and NA61 in elementary proton-proton collisions and all parameters and their default values are explained. This comparison evidently shows that a distinct fragmentation function for leading baryons from soft non-diffractive string processes needs to be employed to get a reasonable agreement with the measured distribution of longitudinal momentum of protons. Studying how the different parameters influence the set of observables leads to a good overall agreement with the available data.
After fixing the parameters of the string routine to proton-proton collisions and achieving a good agreement over the entire SPS energy range, heavy ion collisions are investigated. Since secondary interactions only play an important role for the dynamics of heavy ion collisions, the formation process of string fragments is studied in leadlead collisions. Comparing to experimental data, the best agreement was found, if the formation time is τ form = 1 fm in the rest frame of the respective string fragment and the cross section grows linearly in time during that period. The proton and pion rapidity spectra closely follow the data but the proton multiplicity is overestimated at lower collision energies. Further, the importance of non-isotropic elastic collisions is shown. More forward-backward peaked angular distributions in elastic collisions are essential for reproducing the experimentally observed double peak structure in heavy ion collisions at top SPS energies.
Finally, SMASH is used to obtain event-by-event initial conditions for starting the evolution of the system in terms of hydrodynamics. An energy density and net baryon density profile at the time right after the colliding nuclei have passed through each other is provided. These profiles indicate that realistic fluctuating initial conditions for all conserved charges can be obtained in the future. Another avenue for further research includes to explore dynamical initialization of hydrodynamics via source terms fed by the hadronic transport approach.