1 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 nucleon–nucleon 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 in quantum chromodynamics (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 \(\pi ^{+}\), \(\pi ^{-}\) and \(\pi ^{0}\), and they are the only abundantly produced particles in the few-GeV energy range. Their yields, phase-space distributions, and multi-particle correlations carry information about all stages of the collision.

In this paper, we present experimental data on charged-pion production in centrality-selected Au + Au collisions at \(\sqrt{s_\mathrm{NN}} = 2.4~\hbox {GeV}\) (corresponding to a beam kinetic energy of \(E_{beam} = 1.23~\hbox {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 Sect. 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 Sect. 3. In Sect. 3.2, our result on the pion yield is compared to the pion excitation function for kinetic beam energies ranging from threshold up to 10 A GeV. In Sect. 3.3, we present the centrality and momentum dependence of the parameter \(A_2\) which quantifies the pion production 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 Sect. 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.

2 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,24, 40], has also excellent hadron identification capabilities [25,26,27,28,29,30,31]. HADES is a charged particle detector consisting of a six-coiled 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^\circ \) and \(85^\circ \). 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 (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 \(\approx 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^\circ \) to \(85^\circ \) and Resistive Plate Chambers (RPC) covering polar angles from \(18^\circ \) to \(45^\circ \). 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 [32].

Fig. 1
figure 1

Population of charged particles in the \(\beta \) 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 Eq. (1). Solid curves show the selection of charged pions by a 2D kinematic cuts

Fig. 2
figure 2

Distribution of the uncorrected (raw) \(\pi ^-\) (a) and \(\pi ^+\) (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_{\pi } = 1.3~\hbox {GeV}/c\) is caused by the PID selection (shown as solid black curve in Fig. 1)

The beam consisted of \(\hbox {Au}^{69+}\) ions and had an intensity of approximately \(1.5\times 10^6\) particles/s. It impinged on a 15-fold 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, corresponding to the 0–43% most central events. About 2.1 billion Au + Au events, corresponding to the 40% most central collisions not influenced by the trigger selection, 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 [33]). Interactions of beam particles outside of the target have been rejected by requiring that the main interaction vertex is reconstructed within \(\pm \, 32.5~\hbox {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 [33].

Charged-hadron identification is based on the time-of-flight 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 \(\beta \) 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:

$$\begin{aligned} \beta = \left( \left( {\textstyle \frac{m \, c}{p_{lab}}} \right) ^{2}+1 \right) ^{-1/2}. \end{aligned}$$
(1)

In order to determine the \(2\sigma \) band around the theoretical curve a Gaussian function is fitted to the \(\beta \) distributions in momentum slices. The mean of the Gauss is fixed to the value given by Eq. (1), and the width is a free parameter. In order to avoid contamination of the \(\pi ^+\) sample by protons and of the \(\pi ^{\pm }\) sample by high momentum particles with wrong charge assignment, only pions with laboratory momenta below \(1.3~\hbox {GeV}/c\) are accepted for further analysis. The remaining contamination (impurities) by other particles within our detector acceptance are on average (i.e. averaged over \(\hbox {p}_t\)) below 2% for both \(\pi ^{+}\) and \(\pi ^{-}\). The losses due to the 2 \(\sigma \) cut and the impurities are taken care of by the efficiency and acceptance correction described below. The losses due to the 2 \(\sigma \) 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 [34]. The detector response was simulated using the Geant3 [35] 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 [26]. The resulting correction factors were calculated in 14 rapidity (\(\Delta y = 0.1\)), 60 transverse momentum (\(\Delta p_t = 25~\hbox {MeV}/c\)), and 60 transverse mass (\(\Delta m_t = 25~\hbox {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:

$$\begin{aligned}&\frac{d^{2}N}{dp_{t}\, dy} = C_{1}(y)\, m_t\, p_t\, \exp {\textstyle \left( -\frac{m_{t}c^{2}}{ T_{1}(y) }\right) }\nonumber \\&\quad +\,C_{2}(y)\, m_t\, p_t\, \exp {\textstyle \left( -\frac{m_{t}c^{2}}{ T_{2}(y) }\right) },\ \end{aligned}$$
(2)
$$\begin{aligned}&\frac{1}{m_{t}^{2}}\, \frac{d^{2}N}{dm_{t}\, dy}=C_{3}(y)\, \exp {\textstyle \left( -\frac{m_{t}c^{2} }{ T_{1}(y) } \right) }\nonumber \\&\quad +\, C_{4}(y)\, \exp {\textstyle \left( -\frac{m_{t}c^{2} }{ T_{2}(y) } \right) }. \end{aligned}$$
(3)
Fig. 3
figure 3

Mid-rapidity (a) and forward rapidity (b) transverse momentum distributions (\(p_{t}\)) for \(\pi ^+\) and \(\pi ^-\) mesons for the 10% most central events. The black curves represent the fit of function (2). The spectra are corrected for the efficiency losses and the missing acceptance in azimuth. Red dashed lines mark the peak position of the respective transverse momentum spectrum. The lower part of each panel shows the ratio of data to the fitted function

Fig. 4
figure 4

Reduced transverse mass spectra for \(\pi ^-\) (a) and \(\pi ^+\) (b) mesons at mid-rapidity (\(\mid \hbox {y}_{cm} \mid < 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 Boltzmann 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 Boltzmann function

Examples of the transverse-momentum distributions of \(\pi ^-\) and \(\pi ^+\) mesons (\(dN/dp_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 Fig. 4a, 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 \(\hbox {MeV}/c^{2}\) and 500–800 \(\hbox {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 \(\hbox {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 \(\hbox {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. 3a 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 [36] 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 \(\pi ^{-}\) and \(\pi ^{+}\) in Fig. 3a illustrate these momentum shifts: the maximum of the \(\pi ^{-}\) (\(\pi ^{+}\)) spectrum (represented by red dashed lines) is shifted from their average value of about \(125~\hbox {MeV}/c\) to \(100~\hbox {MeV}/c\) (\(150~\hbox {MeV}/c\)). Based on these results, a separate paper on the determination of the Coulomb potential is in preparation. The deviations between the transverse-momentum (or mass) distributions and the Boltzmann fits cause a systematic underestimate (overestimate) of the extrapolation into the low-momentum region of the \(\pi ^{-}\) (\(\pi ^{+}\)) yield of 4%. The \(\pi ^+\) and \(\pi ^-\) rapidity distributions are obtained by integrating the measured transverse-mass distributions and adding the extrapolated yields (see Fig. 6 in Sect. 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 Sect. 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 Sect. 3.3): (1) The polar angle distribution of the pions is assumed to be of the form (\(1+A_2\ \cos ^{2}\theta \)) with \(A_2\) extracted from our measurement as discussed in Sect. 3.3. 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}\theta + A_4 \cos ^{4}\theta \)) with \(A_2\) and \(A_4\) given by the shape of the polar angle distribution of \(\pi ^{+}\) in p + p interactions at 2 GeV (see fig. 19 in [37]). This is an upper limit, because the very forward and backward preferential emission 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. We assume that the contributions to the systematic uncertainty from the phase space extrapolations are correlated for different event centrality classes and are not relevant when considering centrality-dependent trends in the data.

Fig. 5
figure 5

Reduced transverse mass distributions for negatively (a) and positively (b) charged pions in rapidity bins of \(\Delta y_{cm} = 0.1\) width between − 0.65 and 0.75 for the 0–10% most central events. The most backward rapidity data is shown unscaled while the following rapidity slices are scaled up by successive factors of 10. The solid curves correspond to the two-slope Boltzmann function given by Eq. (3)

Fig. 6
figure 6

Rapidity distribution of negatively (a) and positively (b) charged pions for four (10% wide) centrality classes. Full points are the measured data and open points are the data reflected at \(\hbox {y}_{cm} = 0.\) The error bars indicate the systematic uncertainties. The statistical errors are negligible

3 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 \(\Delta \) decays. \(T_2\) stands for the slope at higher \(m_t\) which is often interpreted as a thermal component [38], but can be also attributed to pions from decays of various broad higher-lying resonances.

Table 1 Measured (“yield”) and extrapolated (“yield \(4\pi \)”) 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 mid-rapidity are listed. Here errors are omitted, because both parameters depend on the chosen fit range and are correlated

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\pi \) 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.

Fig. 7
figure 7

Multiplicities of \(\pi ^{-}\) (a) and \(\pi ^{+}\) (b) as a function of the mean number of participants \(\langle A_{part}\rangle \). The vertical size of the open boxes stands for the systematic uncertainties due to the correction factors and extrapolations. The horizontal size of the open boxes indicates the error on the mean number of participants (for details see [33]). Colored solid lines represent the results of various model calculations for \(\pi ^-\) (\(\pi ^+\)) (see Sect. 4)

3.1 Centrality dependence

Figure 7 shows the dependence of pion yields in Au + Au collisions on the centrality parameterized by the mean number of participants \(\langle A_{part}\rangle \) (see [33]). The yield increases towards more central events. The pion multiplicity per participating nucleon as a function of centrality 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 Fig. 8. Note that the statistical errors are negligible and the systematic ones are not random point-to-point errors, but are strongly correlated common uncertainties in the efficiency correction and some model dependence of the \(\langle A_{part}\rangle \) determination. From the scatter of the data points, we estimate the random point-to-point part of the systematic error to be about 1%, which is smaller than the full 7% systematic error discussed in Sect. 2, but much larger than our statistical errors (typically < 0.1%). This random point-to-point error of 1% is used when fitting the data in order to characterize its evolution with \(\langle A_{part}\rangle \). The scaling of the yields with the mean number of participants \(\langle A_{part}\rangle \) is quantified by the scaling parameter \(\alpha \), where the multiplicity \(M(\pi ^{\pm })\) is \(\propto \langle {A_{part}} \rangle ^\alpha \). We find for \(\pi ^-\) a value of \(\alpha = 0.93\pm 0.01\) and for \(\pi ^+\) a value of \(\alpha = 0.94\pm 0.01\). Hence, we observe for both pion species a significantly weaker than linear scaling with the mean number of participants \(\langle A_{part}\rangle \) and the decrease with increasing centrality observed in Fig. 8 is indeed significant.

The multiplicity of \(\pi ^{-}\) mesons is larger than the one of the \(\pi ^{+}\) due to the neutron over proton excess in the Au nucleus. Using the parameterizations from [39] for the energy dependence of the pion production cross sections in the different isospin channels of nucleon-nucleon interactions, the \(\pi ^{-} / \pi ^{+}\) 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 \pm 0.17.\)

3.2 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 Figs. 9 and 10. 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 multiplicity \(M(\pi )/\langle A_{part}\rangle \) as a function of the beam kinetic energy \(E_{beam}\) where \(M(\pi ) = M(\pi ^{+}) + M(\pi ^{-}) + M(\pi ^{0})\). Unfortunately, the data published on \(\pi ^{0}\) production are much scarcer than for charged pions. When not available, we take the neutral pion multiplicity as \(M(\pi ^{0}) = 0.5 \times \left( M(\pi ^{+}) + M(\pi ^{-}) \right) \) This is the procedure proposed in Ref. [4] making use of approximate isospin symmetry of pion production seen in NN collisions [39]. 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\pi \) assuming isotropic emission from a source at mid-rapidity). Our results for \(M(\pi )\) are listed in Table 2 (together with values of \(\langle A_{part}\rangle \) taken from [33]) and fit well into the overall systematics of the world data.

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 \(\langle A_{part}\rangle \) (taken from [33]) and impact parameter ranges used for centrality selection in the models

We model the energy dependence of the pion multiplicity with a simple second order polynomial (\(a_{0} + a_{1}E_{beam} +a_{2}E_{beam}^{2}\)) (the dashed line in Fig. 9). The resulting parameters are \(a_{0}=-5.36 \times 10^{-2}\), \(a_{1}=1.72 \times 10^{-1}~\hbox {A GeV}^{-1}\), \(a_{2}=-3.08 \times 10^{-3}~\hbox {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~\sigma \)) above our data and the trend of the world data.

Fig. 8
figure 8

Comparison of the centrality dependence of \(M(\pi )/\langle A_{part}\rangle \) in Au + Au collisions to earlier measurements at similar energies. The results from FOPI, E895, and from the BEVALAC Streamer Chamber group (the latter for La + La collisions) have been scaled to 1.23 A GeV; note the suppressed zero on the ordinate. Error bars on HADES data points represent the systematic uncertainties of the pion multiplicity and \(\langle A_{part}\rangle \) added in quadrature

Fig. 9
figure 9

Pion multiplicity \(M(\pi )\) per mean number of participating nucleon \(\langle A_{part}\rangle \) in Au + Au collisions as a function of the kinetic beam energy \(E_{beam}\). The dashed curve is a fit to the data points except for the one labeled “FOPI (Pelte et al.)”, as suggested in [4]. The inset magnifies the energy region around the HADES point

Fig. 10
figure 10

Pion multiplicity per participating nucleon as a function of beam energy for three different systems: C + C (black) [7, 22, 40], Ar + KCl (blue) [4, 7,8,9, 41] 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

Fig. 11
figure 11

Pion multiplicity per participating nucleon as a function of centrality given by \(\langle A_{part}\rangle \). Here only HADES data are shown. The data points for Ar + KCl [24] (open cross) measured at 1.76 A GeV as well as the C + C [26] at 1 A GeV (closed diamond) and at 2 A GeV (open diamond) data points were scaled to be at a beam energy of 1.23 A GeV (for details see text)

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. The scaling to the common energy was done using the fit to the pion excitation function shown as dashed line in Fig. 10. Considering that \(M(\pi )/\langle A_{part}\rangle \) does not change strongly with \(\langle A_{part}\rangle \), the same energy scaling is used independent of the centrality selection.

This kind of scaling is done to data in Fig. 8 using pion energy excitation function from Fig. 9. The centrality dependence of \(M(\pi )/\langle A_{part}\rangle \) from HADES is compared to the energy-scaled FOPI data in Fig. 8. Both data sets exhibit a similar weak variation of \(M(\pi )/\langle A_{part}\rangle \) with \(\langle A_{part}\rangle \) 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 [42]. 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 [43], the values for \(M(\pi )/\langle A_{part}\rangle \) come closer to our data points.

We have measured pion production also in the much lighter systems Ar + KCl [24] and C + C [26] 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. The scaling factors are 0.63 for Ar + KCl at 1.76 A GeV, 1.33 for C + C at 1 A GeV and 0.54 for C + C at 2 A GeV. Figure 11 compares our centrality dependent normalized pion yields from Au + Au collisions with those of the lighter systems. While the normalized pion multiplicity varies only slightly (less than 10%) in collision systems with 100 participating 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.

3.3 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 \(\theta _{cm}\) and momentum \(p_{cm}\). The polar angular distributions (\(dN/d\cos \theta _{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.

Fig. 12
figure 12

Center-of-mass polar angle distributions of negatively (a) and positively (b) charged \(\pi \) mesons. Shown are pions with center-of-mass momenta of 120–800 \(\hbox {MeV}/c\). The curves represent fits with the function given by Eq. (4)

Fig. 13
figure 13

Dependence of the anisotropy parameter \(A_{2}\) on the pion momentum in the center-of-mass system for negatively (upper row) and positively (lower row) charged pions for centrality classes of 0–10% (left) and 30–40% (right), respectively. Points with error bars are the results of fits to the experimental polar angle distributions. The curves represent the results of model calculations (see Sect. 4). Extraction of the anisotropy parameter from models is done within HADES polar angle coverage

In order to quantify the deviation from isotropy, the distributions are fitted with a quadratic function of \(\cos \theta _{cm}\):

$$\begin{aligned} \frac{dN}{d(\cos \theta _{cm})} = C \, (1 + A_{2} \, \cos ^{2}\theta _{cm}), \end{aligned}$$
(4)

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 (Sect. 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} \le 50~\hbox {MeV}/c\)) and increases with momentum for both pion charges, saturating above \(400~\hbox {MeV}/c\) at \(A_{2} \simeq 1.0\). In more peripheral collisions, the saturation value increases up to \(A_{2} \simeq 1.4\). We do not observe the pion energy dependence of \(A_2\) peaking around \(E_{\pi } = 200{-}300\) MeV as seen in Ar + KCl data at 1.8 A GeV [44]. The system-size dependence of the mean \(\langle A_{2}\rangle \) (momentum-averaged over the interval 120–\(800~\hbox {MeV}/c\)) in nuclear collisions is illustrated in Fig. 14. A weak but significant increase of \(\langle A_{2}\rangle \) with decreasing system size is observed for central Au + Au collisions towards Ar + KCl and C + C collisions [26, 45]. One can define the ratio of the anisotropic to isotropic fraction \(\hbox {R}=\langle A_2 \rangle /(3+\langle A_2 \rangle )\). In this way, for the most peripheral class measured by HADES (30–40%), one obtains R \(=\) 0.25, and for the most central \(\hbox {R}=0.14\). Using a linear extrapolation to the most central collisions with \(\langle A_{part} \rangle \) around 400 the anisotropic fraction is reduced to 8%.

Fig. 14
figure 14

Mean anisotropy parameter \(\langle A_{2}\rangle \) as a function of the mean number of participants \(\langle A_{part}\rangle \). Circles represent experimental Au + Au data (closed \(\pi ^{-}\), open \(\pi ^{+}\)). Also shown are earlier HADES results on \(\langle A_{2}\rangle \) as average of \(\pi ^{-}\) and \(\pi ^{+}\) from C + C (triangles) and Ar + KCl (cross) collisions at similar energies. The up-pointing triangle (down-pointing) stands for C + C at the beam energy of 2.0 (1.0) A GeV

Table 3 Charged-pion multiplicities \(\pi ^{-}\) (top part) and \(\pi ^{+}\) (bottom part) in full phase space extracted from the indicated models. In the last column the experimental results extrapolated to \(4\pi \) are given. In the last row the values of the parameter \(\alpha \) are listed as obtained from the fits to the experimental data shown in Fig. 7 and to the corresponding results of the model calculations

Concerning measurements of \(\langle A_{2}\rangle \) in p + p interactions at our beam energy, we found two inconsistent results. Reference [37] reported a measurement of \(dN/d\cos (\theta _{cm})\) of \(\pi ^+\) in p + p \(\rightarrow \) p n \(\pi ^+\) in a bubble chamber experiment at 2 GeV which allows to determine an \(\langle A_{2}\rangle \) value in the range 0.8–1.6 which follows the trend given by our \(A+A\) data in Fig. 14. In reference [44], \(\langle A_{2}\rangle \) higher than 3.0 is reported from an analysis of a counter experiment, which corresponds to \(\hbox {R}\le 0.5\). In view of these contradicting results we refrain from presenting in Fig. 14 a data point of \(\langle A_{2}\rangle \) for p + p interactions.

4 Comparison with transport models

In the following the data will be compared to five state-of-the-art microscopic transport models: “Isospin Quantum Molecular Dynamics” (IQMD vi.c8) [15], “Parton Hadron String Dynamics” (PHSD v4) [16], “Parton-Hadron-Quantum-Molecular Dynamics” (PHQMD v1) [20], “Simulating Many Accelerated Strongly-Interacting Hadrons” (SMASH v1.5.1) [18, 19], and “Giessen Boltzmann–Uehling–Uhlenbeck” (GiBUU, release 2019) [17].Footnote 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 [31]. 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, Figs. 7, 15 and 16).

The yields of charged pions have already been shown as a function of \(\langle A_{part}\rangle \) 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 (\(\alpha \ge \ 1\)), only PHQMD agrees with the significantly weaker scaling observed in our data (see Table 3).

Fig. 15
figure 15

Rapidity distributions of negatively (a) and positively (b) charged pions. The experimental data and the results of five transport models are shown. Full points are the measured data and open points are the data reflected at \(y_{cm} = 0\). The comparison is done for the 10% most central events. The ratio of experimental over model data is displayed in the lower panels

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) \(\pi ^-\) (\(\pi ^+\)) distributions than observed in the experiment.

Fig. 16
figure 16

The mid-rapidity \((y_{cm} \pm 0.05)\) experimental transverse-momentum distributions of negatively (a) and positively (b) charged pions in comparison to models for the \(10\%\) most central events

Fig. 17
figure 17

Decomposition of the transport model calculations into the decay contributions from \(\Delta (1232)\) (solid lines) and from all higher-lying resonances (HR, dashed lines) as function of transverse momentum at mid-rapidity for \(\pi ^{-}\) (a) and \(\pi ^{+}\) (b)

The transverse-momentum distributions are compared in Fig. 16. In the case of negatively charged pions, all models besides IQMD show similar \(p_{t}\) dependence. 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~\hbox {MeV}/c\). For \(\pi ^{+}\), 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 \(\Delta (1232)\) resonances. Higher-mass states are incorporated at different levels, as specified in Table 4 of the appendix (Sect. 6). In the IQMD model, the \(\Delta (1232)\) is the only source of pions and, as discussed in [34], all inelastic NN cross sections are projected onto the excitation of this particular resonance.

Since all models over-predict 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 transverse-momentum distributions of \(\pi ^{-}\) and \(\pi ^{+}\) mesons. The result of this investigation is shown in Fig. 17. For both \(\pi ^{+}\) and \(\pi ^{-}\), PHQMD and SMASH have very similar shapes of the \(\Delta (1232)\) component, but PHSD has a steeper slope than the two. Negative pions from \(\Delta \) 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~\hbox {MeV}/c\). Thus, both the pions from the \(\Delta (1232)\) and those from higher-lying resonances boost up the pion yield, although in different transverse-momentum regions. From this differential comparison with various transport calculations, the origin of the displayed overshoot of pion production remains unclear. As transport models typically make use of measured vacuum NN cross sections (in particular \(\Delta \) production) our observation might indicate that the latter are modified inside the nuclear medium and/or that the re-absorption processes are underestimated in the calculations. In order to scrutinize the different incorporations of \(\Delta \) as well as higher-lying baryonic resonances in the models an analysis of exclusive channels in elementary reactions, as done e.g. in [46], could be a promising avenue to follow.

In Sect. 2, the mid-rapidity transverse-momentum distributions of both, \(\pi ^{+}\) and \(\pi ^{-}\) 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 \(\pi ^{+}\) and \(\pi ^{-}\) 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 \(\pi ^{+} -\) p and \(\pi ^{-} -\) p inelastic cross sections, as pointed out in [47]. 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 \(p_{t}\) could mean that the Coulomb potential assumed in the models is smaller than in the actual collision system.

Fig. 18
figure 18

The measured \(\pi ^{-}/\pi ^{+}\) 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 \(\Delta \) resonances

In Sect. 3.3 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 Figs. 13 and 19, 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.

5 Summary

In summary, we have presented a comprehensive study of \(\pi ^{\pm }\) emission in Au + Au collisions at \(\sqrt{s_{\mathrm{NN}}}=2.4~\hbox {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~\sigma \) 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\) versus \(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 \(\hbox {H}_2\) and nuclear targets as well as new high-precision measurements of the collision system Ag + Ag at \(\sqrt{s_\mathrm{NN}} = 2.55~\hbox {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.