Charged pion production in $\mathbf{Au+Au}$ collisions at $\mathbf{\sqrt{s_{NN}}}$ = 2.4$\mathbf{GeV}$

We present high-statistic data on charged pion emission from Au+Au collisions at $\sqrt{s_{\rm{NN}}}$ = 2.4 GeV (corresponding to $E_{beam}$ = 1.23 A GeV) in four centrality classes in the range 0 - 40$\%$ of the most central collisions. The data are analyzed as a function of transverse momentum, transverse mass, rapidity, and polar angle. Pion multiplicity per participating nucleon decreases moderately with increasing centrality. The polar angular distributions are found to be non-isotropic even for the most central event class. Our results on pion multiplicity fit well into the general trend of the world data, but undershoot by $2.5 \sigma$ data from the FOPI experiment measured at slightly lower beam energy. We compare our data to state-of-the-art transport model calculations (PHSD, IQMD, PHQMD, GiBUU and SMASH) and find substantial differences between the measurement and the results of these calculations.


Introduction
Heavy-ion collisions at relativistic energies allow the study of bulk properties of strongly interacting matter at high temperatures and densities. In these studies, the phase-space distributions of various final-state particles are analyzed and compared to the corresponding distributions in nucleonnucleon interactions in order to disentangle bulk phenomena from the trivial superposition of elementary interactions. The particles present in the final state of relativistic nuclear collisions carry information about the initial state, e.g. the impact parameter, about the properties of the high-density phase of the system, e.g. the pressure and its gradients, and about the expansion and freeze-out conditions of the produced strongly-interacting matter, often called fireball. The final-state particles are either surviving nucleons, nuclear clusters, or newly-produced particles. Being the lightest mesons, pions are the pseudo-Goldstone bosons of SU2 reflecting the approximate spontaneous breaking of chiral symmetry of the QCD. Hence, they are a measure for the degree of excitation in a gas of hadrons [1]. They have an isospin of one and come in the three charge states π + , π − and π 0 , and they are the only abundantly produced particles in the few-GeV energy range. Their yields, phasespace distributions, and multi-particle correlations carry information about all stages of the collision.
In this paper, we present experimental data on chargedpion production in centrality-selected Au+Au collisions at √ s NN =2.4 GeV (corresponding to a beam kinetic energy of E beam = 1.23 A GeV on fixed target). Our results profit from high statistics and thus complement and extend earlier studies of pion production at similar energies and with heavy nuclei [2][3][4][5][6][7][8][9][10][11][12]. They cover rapidity and transverse-momentum (or mass) distributions, as well as derived quantities. Some of the latter are analyzed as a function of the collision centrality. Special emphasis is put on the comparison of the experimental findings with results from microscopic model calculations. Detailed investigations of spectra generated in thermal models in comparison to experimental data are also ongoing and will be discussed in future publications. First results on two-pion correlations have recently been published [13,14]. Multi-pion correlation and collective-flow studies are the subject of separate forthcoming articles.
After describing the experimental setup and the analysis methods in section 2, we present transverse-momentum (p t ) and reduced transverse-mass spectra (m t − m 0 ), as well as rapidity distributions which are used to determine the multiplicities of charged pions in section 3. In section 3.1, our result on the pion yield is compared to the pion excitation function for beam kinetic energies ranging from threshold up to 10 A GeV. It is well established that pions at SIS18 energies are emitted preferably in the forward/backward direction. In section section 3.2, we present the centrality and momentum dependence of the parameter A 2 which quantifies the pion anisotropy. The presentation of our results ends with a detailed comparison of the observables discussed in the previous sections to five state-of-the-art microscopic models [15][16][17][18][19][20] in section 4.
We note in passing that the measured pion yields are important for the normalization of the dielectron data obtained in the same experiment [21]. The low-energy region of the invariant-mass spectra of dielectrons is dominated by the decay products of neutral pions. The charged pions can be used to construct triple differential momentum distributions of neutral pions and their decays which allow to constrain the part of the dielectron spectrum originating from pion decays.

Experimental setup and data analysis
The experiment was performed with the High Acceptance Di-Electron Spectrometer (HADES) at the Schwerionensynchrotron SIS18 at GSI, Darmstadt. HADES, although primarily optimized to measure dielectrons [21][22][23], has also excellent hadron identification capabilities [24][25][26][27][28][29][30]. HADES is a charged particle detector consisting of a sixcoiled magnet producing a toroidal field centered around the beam axis. Six identical detection sections are located between the coil planes and cover polar angles between 18 • and 85 • . Its large azimuthal acceptance varies from 65% at low to 90% at high polar angles. The corresponding losses are corrected for in the analysis. Each sector is equipped with a Ring-Imaging Cherenkov (RICH) detector for electron identification (not relevant for the present analysis) followed by four layers of Mini-Drift Chambers  Fig. 1 Population of charged particles in the β vs. laboratory momentum over charge (p/q) plane for the RPC (a) and TOF (b) detector. Dashed curves correspond to the kinematic correlation for the different particle species as given by equation (1). Solid curves show the selection of charged pions by a 2D kinematic cuts.   Fig. 2 Distribution of the uncorrected (raw) π − (a) and π + (b) yields in the plane spanned by rapidity (y cm ) and reduced transverse mass (m t − m 0 ). The dotted curves depict laboratory polar angles and momenta. The cutoff at p π = 1.3 GeV/c is caused by the PID selection (shown as solid black curve in fig. 1).
(MDCs), two in front of and two behind the magnetic field. The MDCs record space points of the trajectories of charged particles which, together with the known magnetic field, are used to determine the particle momentum. The momentum resolution of the charged pions was found to be ≈ 2.5 % and depends only weakly on laboratory momentum. The arrival time of the charged particles is measured by a scintillator based Time-Of-Flight detector (TOF) covering polar angles from 44 • to 85 • and Resistive Plate Chambers (RPC) covering polar angles from 18 • to 45 • . Their respective time resolutions are 150 and 66 ps. A Pre-Shower detector (behind the RPC, for electron identification) provides additional position information. The experimental setup is described in detail in [31].
The beam consisted of Au 69+ ions and had an intensity of approximately 1.5 × 10 6 particles/s. It impinged on a 15fold segmented gold target, with an integral thickness corresponding to an interaction probability of 1.4 %; the total length of the target assembly was 60 mm. Several triggers were implemented: The minimum-bias trigger was defined by a valid signal in a diamond START detector in front of the target. An online trigger detected interactions and rejected peripheral collisions. It was based on the summed TOF detector multiplicity signal, which selected events with more than about 20 charged particles in this detector. About 2.1 billion Au+Au events, corresponding to the 40% most central collisions, were acquired this way. The centrality of the recorded events has been mapped onto the distribution of the charged particle multiplicity N ch detected in the RPC and TOF detectors (for details see [32]). Interactions of beam particles outside of the target have been rejected by requiring that the main interaction vertex is reconstructed within ±32.5 mm of the center of the target in the direction of the beam. The details of the procedure, which provides a measure of centrality and verifies that the online trigger does not bias the phase-space distributions, are described in [32].
Charged-hadron identification is based on the time-offlight measured with TOF and RPC. Particle velocity is obtained from the measured flight time and flight path. Combining this information with the particle momenta allows to identify charged particles (e.g. pions, kaons or protons) with high significance. Figure 1 shows the population of all charged particles in the plane spanned by their β and laboratory momenta divided by charge for the RPC (left) and TOF (right) detectors. The different particle species are well separated and distributed along the black dashed curves which represent the function: The projections of momentum slices on the β axis are fitted with a Gaussian distribution, the mean of which is fixed to the value given by eq. (1). Their widths are free parameters which are used to select 2 σ bands along the kinematic curves given by eq. (1). In order to avoid contamination of the π + sample by protons and of the π ± sample by high momentum particles with wrong charge assignment, only pions with laboratory momenta below 1.3 GeV/c are accepted for further analysis. The losses due to the 2 σ cut are taken care of by the efficiency and acceptance correction described below. The coverage in the plane spanned by rapidity (y cm ) and reduced transverse mass (m t -m 0 ) of the measured but uncorrected yields is shown in fig. 2. The measured (raw) pion spectra obtained after particle identification must be corrected for the spectrometer acceptance and losses due to the various cuts introduced during track reconstruction. These efficiency and acceptance corrections have been calculated from simulated Au+Au events generated by the UrQMD model [33]. The detector response was simulated using the Geant3 [34] based simulation package including geometry and characteristic of all HADES detectors. Simulations were subjected to the same reconstruction and analysis steps for all centrality classes as the experimental data. The fraction of lost tracks (particles) is quantified by the ratio of the number reconstructed to the number of simulated tracks. The used efficiency and acceptance correction method is described in detail in [25]. The resulting correction factors were calculated in 14 rapidity (∆ y = 0.1), 60 transverse momentum (∆ p t = 25 MeV/c), and 60 transverse mass (∆ m t = 25 MeV/c 2 ) intervals (bins). They are typically on the order of 1.5 -2.0 over our phase space coverage, including the correction for acceptance limitations in azimuthal angle. The bins near the acceptance limits for which the factor was higher than 6.6 (15% efficiency) were excluded from the analysis. The validity of the correction procedure was checked by alternatively using a track-embedding method. Charged pions were generated with a thermal phase space distribution and inverse slope extracted from data. After embedding them into measured events, these events were processed by the standard reconstruction chain and the fraction of lost embedded tracks was calculated. It was found that the resulting correction factors differ by less than 1% from the ones obtained when using plain simulated UrQMD events. Differential pion yields are calculated as the product of raw yields and the correction factors in all accepted bins of y − p t and y − m t and are modeled using the following functions: 1 Examples of the transverse-momentum distributions of  Reduced transverse mass spectra for π − (a) and π + (b) mesons at mid-rapidity (| y cm | < 0.05) for the 10% most central events. The blue curves represent the fit using eq. (3). The red and green curves represent the single Bolzmann function with the parameters T 1 and T 2 (see eq. (3)). The spectra are corrected for efficiency losses and missing acceptance in azimuth. The lower part of each panel shows the ratio of data to the fitted double Bolzmann function. π − and π + mesons (dN/d p t ) at mid-rapidity (a) and forward rapidity (b) are shown in fig. 3 together with fits of the function (2). We have chosen a superposition of two Boltzmann distributions, because the pion reduced transverse-mass spectra deviate from a single exponential as is demonstrated in figs. 4 (a) and (b) which show the m t spectra at mid-rapidity. The parameters T 1 and T 2 account for different slopes at low and high reduced transverse masses, respectively.
The fit procedure starts with independent single Boltzmann fits in separate m t − m 0 ranges (0-300 MeV/c 2 and 500-800 MeV/c 2 ). We use the resulting inverse slope parameters as starting values for T 1 and T 2 for the double Boltzmann fit in the range 0-800 MeV/c 2 . This is done in two steps. First, we require T 2 to be in an interval of few MeV around the value obtained from the single Boltzmann fit in order to improve the fit for T 1 . Then we release the limits and perform the two-slope fit again, extracting the final values of T 1 and T 2 , shown in table 1. The resulting errors are of the order of 1 MeV or smaller but depend on the chosen fit range and are correlated. Therefore, we refrain from quoting them explicitly. The fit function (3) is used to extrapolate the m t spectrum into the low and high regions outside of the acceptance. The particle yield in each rapidity interval is obtained as the sum of the measured and extrapolated yields. The fraction of the latter is a few percent at mid-rapidity and up to 30% towards forward and backward rapidity. The fit ranges were restricted in p t and m t − m 0 to be below 800 MeV/c 2 . The extrapolation into unmeasured regions is subject to systematic uncertainties. A close look to the p t spectra of the positively (negatively) charged pions in fig. 3 (a) reveals that at low momenta the data points are systematically below (above) the fitted curves. We attribute this deviation from a Boltzmann shape to the Coulomb interaction of the pions with the net positive charge of the expanding fireball. This well-established effect [35] causes shifts of the momentum distributions which are different for the positively and negatively charged pions. The former are accelerated leading to a reduced yield at low momenta and the latter are decelerated leading to an increased yield at low momenta. The comparison of the transverse-momentum distributions of π − and π + in fig. 3 (a) illustrate these momentum shifts: the maximum of the π − (π + ) spectrum (represented by red dashed lines) is shifted from their average value of about 125 MeV/c to 100 MeV/c (150 MeV/c). Based on these results, a separate paper on the determination of the Coulomb potential is in preparation. The deviations between the transversemomentum (or mass) distributions and the Boltzmann fits cause a systematic underestimate (overestimate) of the extrapolation into the low-momentum region of the π − (π + ) yield of 4%. The π + and π − rapidity distributions are obtained by integrating the measured transverse-mass distributions and adding the extrapolated yields (see fig. 6 in section 3). The missing yields in the tails of the rapidity distributions are estimated by using the dN/dy shape obtained from five transport models (IQMD, PHSD, GiBUU, SMASH) [15][16][17][18]20] (see section 4). The respective averaged extrapolated contributions to the yields vary from 31% (30%) in central collisions to 36% (33%) in the peripheral ones for negatively (positively) charged pions. The systematic uncertainty of this extrapolation procedure is estimated by considering two extreme assumptions about the polar angle distribution of the pions (see section 3.2): (1) The polar angle distribution of the pions is assumed to be of the form (1 + A 2 cos 2 θ ) with A 2 extracted from our measurement as discussed in section 3.2. This gives a lower limit which is 5% smaller than the estimated yield. (2) The polar angle distribution of the pions is assumed to be of the form (1 + A 2 cos 2 θ + A 4 cos 4 θ with A 2 and A 4 given by the shape of the polar angle distribution of π + in p+p interactions at 2 GeV (see fig. 19 in [36]). This is an upper limit, because the very forward and backward preferences will be more pronounced in p+p than in A + A collisions. This upper limit was found to be 8% and was reduced to 5% to take the higher energy of the p+p data into account.
The statistical errors are negligible due to the large number of analyzed events. The total systematic errors sum up to 7% based on the comparison of the corrected yield in the different sectors (3%), as a measure of the systematic uncertainty of efficiency correction and on the errors coming from the extrapolations in rapidity (5%) and in p t (m t ) (4%). The different systematic uncertainties are added quadratically.

Results
The complete set of the corrected charged-pion measurements is given in fig. 5 together with the fits corresponding to eq. (3). The resulting slope parameters T 1 and T 2 of the transverse-mass spectra at mid-rapidity are listed in table 1. T 1 describes the slope of the low m t part of the spectrum which contains the bulk of the particles and is usually attributed to pions originating from ∆ decays. T 2 stands for the slope at higher m t which is often interpreted as a thermal component [37], but can be also attributed to pions from decays of various broad higher-lying resonances.
The rapidity distributions of m t extrapolated and integrated yields for both charges are presented in fig. 6 for the centrality classes 0-10%, 10-20%, 20-30%, and 30-40%. The 4π yields and their errors listed in table 1 refer to the means and the scatter of the values obtained from the extrapolations in rapidity described in the previous section. Figure 7 shows the centrality dependence of pion yields in Au+Au collisions with the centrality parameterized by the mean number of participants A part . The pion multiplicity per participating nucleon as a function of   Table 1 Measured ("yield") and extrapolated ("yield 4π") particle multiplicities for four (10% wide) centrality classes and for 0-40% range. The statistical errors are negligible. Shown are systematic uncertainties due to the correction procedure and the extrapolation in m t (added quadratically, first error) and extrapolation in rapidity (second error). In addition, the two inverse slope parameters T 1 and T 2 for midrapidity are listed. Here errors are omitted, because both parameters depend on the chosen fit range and are correlated. decreasing centrality and system size increases in our data from 0.13 in the most central Au+Au collisions (0-10%) to 0.14 in the 30-40% interval, see also fig. 9 and fig 10. Note that the statistical errors are negligible and most of the systematic uncertainties partially cancel when comparing the multiplicities in different centrality classes relative to each other. The scaling of the yields with the mean number of participants A part is quantified by the scaling parameter α, where the multiplicity M(π ± ) is ∝ A part α . We find for π − a value of α = 0.93 ± 0.01 and for π + a value of α = 0.94 ± 0.01. Hence, we observe for both pion species a significantly weaker than linear scaling with the mean number of participants A part . Hence, we observe for both pion species a significantly weaker than linear scaling with the mean number of participants A part and the decrease with increasing centrality observed in fig. 9 is indeed significant.
The multiplicity of π − mesons is larger than the one of the π + due to the neutron over proton excess in the Au nucleus. Using the parameterizations from [38] for the energy dependence of the pion production cross sections in the different isospin channels of nucleon-nucleon interactions, the π − /π + ratio can be calculated for our beam energy (1.23 A GeV). The respective value of 1.84 agrees with the experimental finding of 1.83 ± 0.17.

Beam-energy and system-size dependence
The excitation function of pion multiplicities in Au+Au collisions from threshold up to beam energies of 10 A GeV is displayed in fig. 8 and fig. 9. In the past, in this energy range pion data have been collected by the TAPS [6,7], the KaoS [3], the FOPI [4,5] and the E895 [11] experiments. The world data, together with the presented result, are plotted in terms of the normalized total pion  Table 2 Total pion multiplicities for four (10% wide) centrality classes and for the full 0-40% range. Listed are also the mean number of participants A part (taken from [32]) and impact parameter ranges used for centrality selection in the models. multiplicity M(π)/ A part as a function of the beam kinetic energy E beam where M(π) = M(π + ) + M(π − ) + M(π 0 ) and M(π 0 ), when not available, is approximated as M(π 0 ) = 0.5 × [M(π + ) + M(π − )]. The number of participants A part is not a direct observable and the methods used for its estimation vary between different experiments, which puts limits on the accuracy of such a comparison (Note that all published TAPS data were extrapolated to 4π assuming isotropic emission from a source at mid-rapidity).
Our results for M(π) are listed in table 2 (together with values of A part taken from [32]) and fit well into the overall systematics of the world data. We model the energy dependence of the pion multiplicity with a simple second order polynomial (a 0 + a 1 E beam + a 2 E 2 beam ) (the dashed line in fig. 8). The resulting parameters are a 0 = −5.36 × 10 −2 , a 1 = 1.72 × 10 −1 A GeV −1 , a 2 = −3.08 × 10 −3 A GeV −2 . It turns out that the early result from the FOPI experiment at 1 A GeV [5] is significantly below the value suggested by the world data. This has already been discussed by the FOPI collaboration in [4] and attributed to detector effects. Therefore we exclude this measurement in the fit. Furthermore, the FOPI data points at 1.2 and 1.5 A GeV are 25% (2.5 σ ) above our data and the trend of the world data. In order to compare system size dependence of pion production measured by HADES with other experiments at slightly different energies, it is necessary to scale their results to the same energy of 1.23 A GeV. This kind of scaling is done to data in fig. 9 using pion energy excitation function from fig. 8. The centrality dependence of M(π)/ A part from HADES is compared to the energy-scaled FOPI data in fig. 9. Both data sets exhibit a similar weak variation of M(π)/ A part with A part but differ significantly in absolute values. Also shown is the reduced pion multiplicity scaled to 1.23 A GeV obtained by E895 [11] and the BEVALAC Streamer Chamber group for La+La collisions [10]. The FOPI results on pion multiplicity lie even above the La+La data in spite of the known trend that the reduced pion multiplicity increases with decreasing system size at a given energy [41]. This discrepancy might be (partly) explained by the previously mentioned different methods in the estimation of the number of participants. Indeed, if one uses for the FOPI data the charge balance of the reconstructed charged particles as a centrality estimator instead [42], the values for M(π)/ A part come closer to our data points.
We have measured pion production also in the much lighter systems Ar+KCl [23,38] and C+C [25] at similar  beam energies. In order to examine the system-size dependence of the pion production, the yields in Ar+KCl and C+C were scaled to the beam energy of 1.23 A GeV using the curves displayed on fig. 10 Fig. 10 Pion multiplicity per participating nucleon as a function of beam energy for three different systems: C+C (black) [7,22,39], Ar+KCl (blue) [4,[7][8][9]40] and Au+Au (red) [4,6,7,11]. The curves are polynomial fits to these data used to interpolate the multiplicities as a function of bombarding energy for corresponding systems.  tiplicity varies only slightly (less than 10%) in collision systems with 100 participant nucleons and more it increases by up to 30% at 40 participants and is almost a factor of two higher at 6 participants in the light C+C system.

Polar angle distribution
So far we have parameterized the phase space of the pions by rapidity and transverse momentum or reduced transverse mass. In the following, we span the pion phase space using their c.m. polar angle θ cm and momentum p cm . The polar angular distributions (dN/d cos θ cm ) of charged pions in heavy-ion and N+N collisions are known to be non-isotropic with a preference for forward and backward angles. This feature is also illustrated in fig. 12 for the four centrality classes. Fully isotropic emission would correspond to a flat distribution.
In order to quantify the deviation from isotropy, the distributions are fitted with a quadratic function of cos θ cm : where C is a normalization factor and A 2 a parameter which quantifies the forward/backward preference of pion emission. The solid curves in fig. 12 are the results of the fits. The extraction of the anisotropy parameters and their comparison with models is done within the HADES acceptance. The momentum dependence of the A 2 parameter for the most central (0-10%) and the most peripheral (30-40%) class is shown in fig. 13. The corresponding plots for 10-20 % and 20-30 % can be seen in fig. 19 of the appendix (section 6). The overall trend is that in the experimental data for the most central events (0-10%), A 2 is compatible with zero at low momenta (p cm ≤ 50 MeV/c) and increases with momentum for both pion charges, saturating above 400 MeV/c at A 2 1.0. In more peripheral collisions, the saturation value increases up to A 2 1.4. We do not observe the pion energy dependence of A 2 peaking around E π = 200-300 MeV as seen in Ar+KCl data at 1.8 A GeV [43]. The system-size dependence of the mean A 2 (momentum-averaged over the interval 120-800 MeV/c) in nuclear collisions is illustrated in fig. 14. A weak but significant increase of A 2 with decreasing system size is observed for central Au+Au collisions towards Ar+KCl and C+C collisions. One can define the ratio of the anisotropic to isotropic fraction R= A 2 /(3+ A 2 ). In this way, for the most peripheral class measured by HADES (30-40%), one obtains R=0.25, and for the most central R=0.14. Using a linear extraploation to the most central collisions with A part around 400 the anisotropic fraction is reduced to 8%.
Concerning measurements of A 2 in p+p interactions at our beam energy, we found two inconsistent results. Reference [36] reported measurement of dN/d cos(θ cm ) of π + in p + p → pnπ + in a bubble chamber experiment at 2 GeV which allows to determine an A 2 value in the range 0.8 -1.6 which follows the trend given by our A + A data in fig. 14. In reference [43], A 2 higher than 3.0 is reported from an analysis of a counter experiment, which corresponds to R≤0.5. In view of these contradicting results we refrain from presenting a data point of A 2 for p+p interactions.

Comparison with transport models
In the following the data will be compared to five state-ofthe-art microspopic transport models: "Isospin Quantum Molecular Dynamics" (IQMD vi.c8) [15], "Parton Hadron String Dynamics" (PHSD v.4) [16], "Parton-Hadron-Quantum-Molecular Dynamics" (PHQMD v1) [20], "Simulating Many Accelerated Strongly-Interacting Hadrons" (SMASH) [18,19], and "Giessen Boltzmann-Uehling-Uhlenbeck" (GiBUU, release 2019) [17] 1 . We study the differences between model predictions as observed in the rapidity, transverse momentum, and polar angle distributions of charged pions, in the trends of the A 2 parameter as a function of p cm and system size, and in the predicted abundance of resonances as well. In the models, the selection of four centrality classes is done by selecting the corresponding impact-parameter intervals (see table 2) according to the values estimated in [30]. We find that all models over-predict the pion yields for all centralities by factors ranging from 1.2 to 2.1 (see table 3, fig. 7 and fig.  15). The yields of charged pions have already been shown as a function of A part in fig. 7 together with the results of the model calculations. Most of the model calculations give a linear or slightly stronger than linear dependence (α ≥ 1), only PHQMD agrees with the significantly weaker scaling observed in our data (see table 3).
The yields and shapes of the rapidity distributions are compared in fig. 15. In the lower panels the ratio of the experimental and model data quantify both the yield excess and differences in shape. All models tend to have slightly wider (narrower) π − (π + ) distributions than observed in the experiment. The transverse-momentum distributions are  Table 3 Charged-pion multiplicities π − (top part) and π + (bottom part) in full phase space extracted from the indicated models. In the last column the experimental results extrapolated to 4π are given. In the last row the values of the parameter α are listed as obtained from the fits to the experimental data shown in fig. 7 and to the corresponding results of the model calculations.  compared in fig. 16. In the case of negatively charged pions, all models besides IQMD show similar p t dependences. However, the slopes are clearly steeper than observed in the data, which is in particular evident from the ratio plot (lower panels). IQMD, on the other hand, predicts a p t dependence similar to the data and thus the data/theory ratio has a rather flat dependence on p t up to 600 MeV/c. For π + , a rather flat data/theory ratio is observed for all models.
The high p t region resulting from IQMD and PHQMD calculations are steeper than those of the experimental data and from the other models. In all of the presented models, the bulk of the pion yield stems from the decays of ∆ (1232) resonances. Higher-mass states are incorporated at different levels, as specified in table 4 of the appendix (section 6). In the IQMD model, the ∆ (1232) is the only source of pions and, as discussed in [33], all inelastic NN cross sections are projected onto the excitation of this particular resonance.
Since all models overpredict the pion multiplicity significantly, one may ask whether the treatment of the resonances is at the origin of this "pion excess". To this end, we have studied the resonance contributions to the transversemomentum distributions of π − and π + mesons. The result of this investigation is shown in fig. 17. For both π + and π − , PHQMD and SMASH have very similar shapes of the ∆ (1232) component, but PHSD has a steeper slope than the two. Negative pions from ∆ in GiBUU are consistent with SMASH and PHQMD while positive pions have a less steep slope. In IQMD negative pions are shifted towards lower p t . Overall, heavy resonances are a minor source of pions at low p t while they contribute to about 50 % to the yield above 600 MeV/c. Thus, both the pions from the ∆ (1232) and those from higher-lying resonances boost up the pion yield, although in different transverse-momentum regions. Comparing models to an analysis of exclusive channels in elementary reactions, as done e.g. in [44], might be a promising avenue to follow in order to scrutinize the different incorporations of higher-lying baryonic resonances in the models. In section 2, the mid-rapidity transverse-momentum distributions of both, π + and π − were found to deviate significantly from the Boltzmann shape (2) at low p t . These deviation were attributed to the Coulomb interaction with the central positive charge distribution. This subject is addressed again in fig. 18 which shows the reduced transverse-mass dependence of π + and π − ratio for data and model calculations. The ratio remains rather constant for PHSD and SMASH, but not for PHQMD. The small variations might be due to the different energy dependence of π + − p and π − − p inelastic cross sections, as pointed out in [45]. The experimental ratio exhibits a monotonic increase with decreasing m t and a steep rise below 100 MeV/c. A rise at low transverse momentum is also visible in IQMD and GiBUU, but much less pronounced. These two transport models are the only ones which have the Coulomb interaction implemented in their codes. The significant deviations from data at low pt could mean that the Coulomb potential assumed in the models is smaller than in the actual collision system. In section 3.2 the polar angle distributions of charged pions were characterized by the parameter A 2 . Its dependence on p cm for four 10 % centrality classes was shown in fig. 13, together with the results of the model calculations. For the most central (0-10 %) class (upper row) the differences between data and models are moderate except for PHSD which exhibits a significant variation with p cm . The situation becomes more involved if one takes the other three centrality classes into account (see fig. 13). Here, the structures in the data have slightly higher amplitudes. The models, however, develop a strong single oscillation, whose amplitude increases with decreasing centrality.  Fig. 18 The measured π − /π + -ratio at mid-rapidity as a function of the reduced transverse mass in comparison to several transport model predictions. The dashed line corresponds to the pion ratio of 1.94 from a simple isobar model where one assumes that all pions are produced via ∆ resonances.

Summary
In summary, we have presented a comprehensive study of π ± emission in Au+Au collisions at √ s NN = 2.4 GeV. The data is discussed as a function of the transverse momentum (reduced transverse mass) as well as the rapidity in four centrality classes covering the 40% most central events. We find that our results on the pion multiplicity fit well into the overall systematic of the world data, but are lower by 2.5 σ in yield when comparing to data from a FOPI experiment performed at a slightly lower collision energy. The pion multiplicity per participating nucleon increases as a function of decreasing centrality and system size from 0.13 in the 0−10 % most central Au+Au collisions to 0.14 at 30−40 % centrality and to 0.23 in minimum-bias C+C collisions with six participant nucleons only. The polar angular distributions are found to be non-isotropic even for the most central event class. The experimental data are compared to several state-of-art transport model calculations. All of the models substantially overestimate the absolute yields and only PHQMD is able to describe the moderate decrease of pion multiplicity per participating nucleon as a function of increasing centrality. While the shape of the rapidity distribution is fairly well reproduced by all models, the shape of the reduced transverse-mass spectra as well as the anisotropy parameter A 2 (p cm ) are, however, not described satisfactorily by any. Decomposing the contributions to the pion spectra of the various resonances implemented in the different models, we find significant variations in both relative yield and shape of these different pion sources. Data on pion-induced reactions on H 2 and nuclear targets as well as new high-precision measurements of the collision system Ag+Ag at √ s NN =2.55 GeV, taken as part of FAIR Phase0, will become available soon. These upcoming results will be used to further extend the world database on pion production and to constrain model calculations ever more stringently.
Here supplementary information to the main text is provided. Figure 19 extends fig. 13, showing for two more centrality selections (10-20% and 20-30%) the dependence of the anisotropy parameter A 2 on c.m. momentum p cm for negatively and positevely charged pions, respectively. Table 4 lists the baryon and meson resonances which were included in the transport calculations discussed in section 4.