1 Introduction

Two-particle intensity interferometry of hadrons is widely used to study the spatio-temporal size, shape and evolution of their source in heavy-ion collisions or other reactions involving hadrons (for a review see Ref. [1]). The technique, pioneered by Hanbury Brown and Twiss [2] to measure angular radii of stars, later on named HBT interferometry, is based on the quantum-statistical interference of identical particles. Goldhaber et al. [3] first applied intensity interferometry to hadrons. In heavy-ion collisions, the intensity interferometry does not allow to measure directly the reaction volume, as the emission zone, changing shape and size in the course of the collision, is affected by dynamically generated space-momentum correlations (e.g. radial expansion after the compression phase or resonance decays). Thus, intensity interferometry generally does not yield the proper source size, but rather an effective “length of homogeneity” [1]. It measures source regions in which particle pairs are close in momentum, so that they are correlated as a consequence of their quantum statistics or due to their two-body interaction. At energies in the GeV region, the measured particles can originate from many different processes. Therefore, the intensity interferometry may provide additional information to the understanding of reaction mechanisms which finally determine the particle emission sources.

In general, the sign and strength of the correlation is affected by (i) the strong interaction, (ii) the Coulomb interaction if charged particles are involved, and (iii) the quantum statistics in the case of identical particles (Fermi-Dirac suppression for fermions, Bose–Einstein enhancement for bosons). In the case of \(\pi \pi \) correlations, the mutual strong interaction appears to be negligible [4] compared to the effects (ii) and (iii).

It is worth emphasizing that, before the measurement presented here, only preliminary data [5] of identical-pion HBT data exist for a large symmetric collision system (like Au + Au or Pb + Pb) at a beam kinetic energy of about \(1A~\text{ GeV }\) (fixed target, \({\sqrt{s_{\mathrm{NN}}}} = 2.3~\hbox {GeV}\)).Footnote 1 For the somewhat smaller system La + La, studied at \(1.2A~\hbox {GeV}\) with the HISS spectrometer at the Lawrence Berkeley Laboratory (LBL) Bevalac, pion correlation data were reported by Christie et al. [6, 7]. An oblate shape of the pion source and a correlation of the source size with the system size were found. Also, pion intensity interferometry for small systems (Ar + KCl, Ne + NaF) was studied at 1.8\(A~\hbox {GeV}\) at the LBL Bevalac using the Janus spectrometer by Zajc et al. [8]. Both groups made first attempts to correct the influence of the pion-nuclear Coulomb interaction on the pion momenta. The effect on the source radii, however, was found negligible for their experiments.

In this article we report on the investigation of \(\pi ^-\pi ^-\) and \(\pi ^+\pi ^+\) correlations at low relative momenta in Au + Au collisions at 1.23\(A~\text{ GeV }\)  (fixed target, \(p_\mathrm{proj}/A=1.96\) \(\text{ GeV }/c\), \(\sqrt{s_\mathrm{NN}} = 2.41~\hbox {GeV}\)), continuing our previous femtoscopic studies of smaller collisions systems [9,10,11]. In a recent letter [12] presenting only the results of the azimuthally-integrated HBT analysis for central collision, we reported a substantial charge-sign difference of the source radii, particularly pronounced at low transverse momenta, and a smooth extrapolation of the \(\sqrt{s_\mathrm{NN}}\) dependence of the source parameters towards low energies. Here, we present the complete HBT analysis, including azimuthal-angle and centrality dependences. In Sect. 2 we shortly describe the experiment. In Sect. 3 we define the correlation function and discuss possible distortions to it. In Sect. 4 we present the three-dimensional pion emission source resulting from the correlation analysis, its dependences on collisions centrality and kinematic quantities and compare our observations to the findings of other experiments. Finally, we summarize our results in Sect. 5 and give an outlook.

2 The experiment

The measurement was performed with the High Acceptance Di-Electron Spectrometer (HADES) at the Schwerionensynchrotron SIS18 at GSI, Darmstadt. HADES, although primarily optimized to measure di-electrons [13], offers also excellent hadron identification capabilities [14,15,16,17,18,19]. The setup of the HADES experiment is described in detail in [20]. HADES is a charged particle detector consisting of a six-coil toroidal magnet centered around the beam axis and six identical detection sections located between the coils and covering polar angles between \(18^{\circ }\) and \(85^{\circ }\). Each sector is equipped with a Ring-Imaging Cherenkov (RICH) detector followed by four layers of Mini-Drift Chambers (MDCs), two in front of and two behind the magnetic field, as well as a scintillator Time-Of-Flight detector (TOF) (\(44^{\circ }\) – \(85^{\circ }\)) and Resistive Plate Chambers (RPC) (\(18^{\circ }\) – \(45^{\circ }\)). TOF, RPC, and Pre-Shower detectors (behind RPC, for e\(^\pm \) identification) were combined into a Multiplicity and Electron Trigger Array (META). Charged hadron identification is based on the time-of-flight measured with TOF and RPC, and on the energy-loss information from TOF as well as from the MDC tracking chambers. Electron candidates are in addition selected via their signals in the RICH detector. Combining this information with the momentum, p, as determined from the deflection of the tracks in the magnetic field, allows to identify charged particles (e.g. pions, kaons or protons) with a high significance. Particles are assumed to be identified as pions if their velocity, \(\beta \), is found within a \(3\sigma \) window around the theoretical expectation, \(\beta _\pi =p/\sqrt{p^2+m_\pi ^2}\). The corresponding cut parameters \(\sigma \) are derived from Gaussian fits to the velocity distribution in slices of momentum. (We use units with \(\hbar =c^2=1\).)

Several triggers are implemented. The minimum bias trigger is defined by a signal in a diamond START detector in front of the 15-fold segmented gold target. In addition, online Physics Triggers (PT) are used, which are based on hardware thresholds on the TOF signals, proportional to the event multiplicity, corresponding to at least 5 (PT2) or 20 (PT3) hits in the TOF. Events are selected offline by requiring that their global event vertex is inside the target region, i.e. within an interval of 65 mm along the beam axis. About 2.1 billion Au + Au collisions corresponding to the 43% most central events are taken into account for the present correlation analysis.

The centrality determination is based on the summed number of hits detected by the TOF and the RPC detectors. The measured events are divided in centrality classes corresponding to intervals of the integrated cross section. Note that, as a result of the analogue hardware threshold of the online trigger, the most peripheral class can be extended to 45% (or even further), i.e. beyond the above given effective value.

The respective average impact parameters and mean numbers of nucleons participating in the formation of the nuclear fireball, \(\langle A_\mathrm{part} \rangle \), as deduced from extensive Monte-Carlo (MC) calculations [21, 22] with the Glauber model [23] are summarized in Table 1.

The determination of the reaction plane angle, \(\phi _{\mathrm{RP}}\), is based on the measurement of the charged projectile spectator fragments (mostly \(Z=1,\,2\)) by their position, flight time and energy deposit. They are detected by a three square meter scintillator hodoscope 7 m downstream the target, consisting of 288 scintillator cells. It covers polar angles from 0.3 to 7.2\(^{\circ }\). Due to the dispersion of the event plane a resolution correction needs to be applied (cf. Sect. 4.4).

Finally, we want to note that, as result of both the neutron excess of the collision system and the different detector acceptances of negatively and positively charged particles, about five to six times more \(\pi ^-\pi ^-\) pairs than \(\pi ^+\pi ^+\) pairs are measured.

Table 1 The mean number of participants (2nd column) and the average impact parameter (3rd column) calculated [21] with Glauber MC simulations [22] for the centrality classes listed in the 1st column. The 4th and 5th columns give the corresponding 1st and 2nd order reaction plane resolutions (cf. Sect. 4.4), respectively, calculated according to Ref. [24]
Full size table

3 The correlation function

Generally, the two-particle correlation function is defined as the ratio of the probability \(P_2(\varvec{p}_1,\varvec{p}_2)\) to measure simultaneously two particles with momenta \(\varvec{p}_1\) and \(\varvec{p}_2\) and the product of the corresponding single-particle probabilities \(P_1(\varvec{p}_1)\) and \(P_1(\varvec{p}_2)\) [1],

$$\begin{aligned} C(\varvec{p}_1, \varvec{p}_2) = \frac{P_2(\varvec{p}_1,\varvec{p}_2)}{P_1(\varvec{p}_1) P_1(\varvec{p}_2)}. \end{aligned}$$
(1)

Experimentally this correlation is formed as a function of the momentum difference between the two particles of a given pair and quantified by taking the ratio of the yields of ‘true’ pairs (\(Y_\mathrm{true}\)) and uncorrelated pairs (\(Y_\mathrm{mix}\)). \(Y_\mathrm{true}\) is constructed from all particle pairs in the selected phase space interval from the same event. \(Y_\mathrm{mix}\) is generated by event mixing, where particle 1 and particle 2 are taken from different events. Care was taken to mix particles from similar event classes in terms of multiplicity, vertex position and reaction plane angle. The events are allowed to differ by not more than 10 in the number of the RPC + TOF hit multiplicity (i.e. corresponding to the typical multiplicity uncertainty as deduced from simulations for fixed impact parameter [21]), 1.2 mm in the z-vertex coordinate (amounting to less than one third of the spacing between target segments), and 10\(^{\circ }\) in azimuthal angle relative to the reaction plane (significantly below the event plane resolution, cf. Table 1), respectively. The momentum difference is decomposed into three orthogonal components as suggested by Podgoretsky [25], Pratt [26] and Bertsch [27]. The three-dimensional correlation functions are projections of Eq. (1) into the (out, side, long)-coordinate system, where ‘out’ means along the pair transverse momentum, \(\varvec{p}_{\mathrm{t},\,12} = \varvec{p}_{\mathrm{t},\,1}+ \varvec{p}_{\mathrm{t},\,2}\), ‘long’ is parallel to the beam direction z, and ‘side’ is oriented perpendicular to the other two directions. The particles forming a pair are boosted into the longitudinally comoving system, where the z-components of the momenta cancel each other, \(p_{z_{1}}+p_{z_{2}}=0\). This system choice allows for an adequate comparison with correlation data taken at very different, usually much higher, collision energies, where the distribution of the rapidity, \(y=\tanh ^{-1}{(\beta _\mathrm{z})}\), of produced particles is found to be not as narrow as in the present case but largely elongated. (Here, \(\beta _\mathrm{z}=p_\mathrm{z}/E\), \(E=\sqrt{p^2+m_0^2}\) and \(m_0\) are the longitudinal velocity, the total energy and the rest mass of the particle, respectively. ) Hence, the experimental correlation function is given by

$$\begin{aligned} C(q_{{\text {out}}},q_{{\text {side}}},q_{{\text {long}}}) = \mathcal{N} \,\frac{Y_{\mathrm{true}}(q_{{\text {out}}},q_{{\text {side}}},q_{{\text {long}}})}{Y_{\mathrm{mix}}(q_{{\text {out}}},q_{{\text {side}}},q_{{\text {long}}})}, \end{aligned}$$
(2)

where \(q_i=(p_\mathrm{1,\,i}-p_\mathrm{2,\,i})/2\) (i = ‘out’, ‘side’, ‘long’) are the relative momentum components, and \(\mathcal{N}\) is a normalization factor which is fixed by the requirement \(C \rightarrow 1\) at large relative momenta, where the correlation function is expected to flatten out at unity. The statistical errors of C are dominated by those of the true yield, since the mixed yield is generated with much higher statistics. Analogously to Eq. (2), the experimental one-dimensional correlation function,

$$\begin{aligned} C(q_{{\text {inv}}}) = \mathcal{N}' \,\frac{Y_{\mathrm{true}}(q_{{\text {inv}}})}{Y_{\mathrm{mix}}(q_{{\text {inv}}})}, \end{aligned}$$
(3)

is generated by projecting Eq. (1) onto the Lorentz-invariant relative momentum,

$$\begin{aligned} q_{{\text {inv}}}=\frac{1}{2}\sqrt{(\varvec{p}_1 - \varvec{p}_2)^2-(E_1-E_2)^2}, \end{aligned}$$
(4)

where \(E_i=\sqrt{p^2_i+m^2_i}\) and \(m_i\) (\(i=1,\,2\)) are the total energies and the rest masses of the particles forming the pair, respectively. Note that for two particles of equal mass, \(q_{{\text {inv}}}\) is identical in magnitude to the single particle momenta in the rest frame of the pair.

3.1 Treatment of close-track effects

Two different methods to correct for the possible bias due to close-track effects introduced by the limitations of the HADES detector + track finding procedures are investigated and found to agree within statistical fluctuations.

Method A, called double-ratio method, is based on UrQMD [28] transport model simulations of Au + Au collisions at 1.23\(A~\text{ GeV }\) fully transported through the HADES implementation HGeant of the ‘Detector Description and Simulation Tool’ GEANT3.21 [29]. Since UrQMD does not incorporate quantum-statistical or Coulomb effects, any identical-pion correlation function generated from the corresponding HGeant output according to Eq. (3) is expected to be flat at unity. Deviations from that value are considered to result from close-track effects. Hence, dividing the experimental and the simulated correlation functions should allow for a reasonable correction of the former one. However, this method is statistically limited because of the present availability to us of only 100 million UrQMD events, i.e. 20 times less than the experimental data sample. Thus, the statistical errors of the corrected correlation function would be dominated by those of the simulation, a drawback becoming particularly important when investigating the correlation function in a multi-dimensional parameter space. Also, it could not be fully guaranteed that HGeant is realistic enough in reproducing the behaviour of nearby tracks. These simulations also showed that there are no significant long-range correlations, usually attributed either to energy-momentum conservation in correlation analyses of small systems or to minijet-like phenomena at high energies. Therefore, this method serves for cross checks but is not involved in the results presented in Sect. 4.

Method B implements appropriate selection conditions on the META-hit and MDC-layer level, i.e. by discarding pairs which hit the same META cell, and by excluding for particle 2 three successive wires symmetrically around the MDC wire fired by particle 1. This method was tested with simulations based on UrQMD + HGeant and a detailed description of the detector response, to firmly exclude any close-track effect. Also broader exclusion windows have been investigated, but no significant improvement was found. Though there is a certain amount of pairs with small relative momenta getting lost due to this condition (about 50% for \(q_{{\text {inv}}}<40~\text{ MeV }/c\)), its superior statistical significance still clearly favors this method over the double-ratio method. Consequently, Method B, applied to both the true and the mixed-event yields, is used throughout the analysis [30] presented in this article.

3.2 Parameterization of one-dimensional correlation functions

The fits to the one-dimensional correlation function are performed with the function

$$\begin{aligned} C(q_{{\text {inv}}}) = N \left[ 1-\lambda _{{\text {inv}}}+ \lambda _{{\text {inv}}}K_\mathrm{C}(q_{{\text {inv}}},R_{{\text {inv}}}) C_\mathrm{qs}(q_{{\text {inv}}})\right] , \end{aligned}$$
(5)

using a Gaussian function for the quantum-statistical (Bose–Einstein) part,

$$\begin{aligned} C_\mathrm{qs}(q_{{\text {inv}}})=1+\exp {(-(2q_{{\text {inv}}}R_{{\text {inv}}})^2)}. \end{aligned}$$
(6)

The influence of the mutual Coulomb interaction in Eq. (5) is separated from the Bose–Einstein part by including in the fits the commonly used Coulomb correction by Sinyukov et al. [31]. The Coulomb factor \(K_\mathrm{C}\) results from the integration of the two-pion Coulomb wave function squared over a spherical Gaussian source of fixed radius. This radius is iteratively approximated by the result of the corresponding fit to the correlation function. The parameters N and \(\lambda \) in Eq. (5) represent a normalization constant and the fraction of correlated pairs, respectively.

All fits performed to one-dimensional and three-dimensional (cf. Sect. 3.3) correlation functions involve a log-likelihood minimization [32]. No significant differences are observed when using a \(\chi ^2\) minimization.

3.3 Parameterization of three-dimensional correlation functions

The three-dimensional experimental correlation function is fitted with the function

$$\begin{aligned}&C(q_{{\text {out}}},\,q_{{\text {side}}},\,q_{{\text {long}}}) \nonumber \\&\quad = N \big [1-\lambda _{{\text {osl}}}+ \lambda _{{\text {osl}}}K_\mathrm{C}(\hat{q},R_{{\text {inv}}}) C_\mathrm{qs}(q_{{\text {out}}},\,q_{{\text {side}}},\,q_{{\text {long}}})\big ], \end{aligned}$$
(7)

where, for azimuthally-integrated analyses (Sect. 4.3) at midrapidity (\(y_\mathrm{cm}=0.74\)),

$$\begin{aligned}&C_\mathrm{qs}(q_{{\text {out}}},\,q_{{\text {side}}},\,q_{{\text {long}}}) \nonumber \\&\quad =1+\exp {(-(2q_{{\text {out}}}R_{{\text {out}}})^2-(2q_{{\text {side}}}R_{{\text {side}}})^2}\nonumber \\&\qquad -\,(2q_{{\text {long}}}R_{{\text {long}}})^2) \end{aligned}$$
(8)

represents the quantum-statistical part of the correlation function and

$$\begin{aligned} \hat{q}= q_{{\text {inv}}}(q_{{\text {out}}},\,q_{{\text {side}}},\,q_{{\text {long}}},\bar{k}_\mathrm{t}) \end{aligned}$$
(9)

is the average value of the invariant momentum difference for given intervals of the relative momentum components and \(k_\mathrm{t}\). Actually, Eq. (8) can be written in a more general way,

$$\begin{aligned} C_\mathrm{qs} = 1+\exp {\left( -4 \sum _{i,\,j} q_i R^2_{ij} q_j\right) }, \end{aligned}$$
(10)

as used for the azimuthally-sensitive HBT analyses presented in Sect. 4.4. For symmetry reasons [33] the non-diagonal (\(i \ne j\)) elements of \(R_{ij}\), comprising the combinations ‘out’-‘side’ and ‘side’-‘long’, vanish when azimuthally and rapidity integrated analyses are performed [34, 35], as is done in Sect. 4.3. The ‘out’-‘long’ component, however, should have a finite value depending on the degree of symmetry of the detector-accepted rapidity distribution w.r.t. midrapidity. We studied this effect by including in Eq. (8) an additional term \(-\,2 q_{{\text {out}}}(2 R_{{\text {{out long}}}})^2 q_{{\text {long}}}\), where the prefactor accounts for both non-diagonal terms, ‘out-long’ and ‘long-out’. We found only marginal differences in the fits which resulted, for all centrality and transverse-momentum classes, in rather small values of \(R^2_{{\mathrm{out long}}}< 1\,\hbox {fm}^2\), which in many cases are consistent with zero within uncertainties. This finding is not surprising, since the rapidity distribution of charged pions accepted by HADES is almost centered at midrapidity and rather narrow, i.e. extending just over about \(\pm 0.7\) rapidity units.

Finally, for all results presented in Sect. 4, we restricted the pair rapidity to an interval \(\vert y - y_\mathrm{cm} \vert <0.35\), within which dN/dy does not vary by more than 10%.

3.4 Momentum resolution correction

The effect of the finite momentum resolution of the HADES tracking system [20] is studied with dedicated simulations. For that purpose, pairs of identical pions (resulting from the two-body decay of a fictive mother particle with its mass chosen to produce relative momenta of interest, e.g. \(q_{{\text {inv}}}=10,\,20,\,30,\ldots ~\text{ MeV }/c\)) are simulated with the event generator Pluto [36] and subsequently tracked through HGeant, the latter one modelling the HADES detector with its granularity and momentum resolution. The relative-momentum distribution of the pion pairs delivered by the simulation is fitted with a Gaussian function. The resulting widths \(\sigma _q(q_\mathrm{i})\) (i = ‘inv’, ‘out’, ‘side’, ‘long’) are folded into the fit functions (Eqs. (5) and (7)) to account for the slight resolution-induced decrease of both R and \(\lambda \). A typical resolution amounts to \(\sigma _q(q_{{\text {inv}}}=20~\text{ MeV }/c) \simeq 2~\text{ MeV }/c\) corresponding to a radius shift of \(\varDelta R/R \simeq +2\%\) after the correction.

3.5 Systematic error estimate

The main contribution to the systematic uncertainties of the results presented in the subsequent section is due to the fluctuation of the fit results when varying the fit range over the respective relative-momentum quantity. The standard fit range was chosen to an interval extending from 4 to \(80~\text{ MeV }/c\). We varied the lower limit to \(8~\text{ MeV }/c\) and the upper one from 60 to \(100~\text{ MeV }/c\). Typical radius changes of 0.1–0.3 fm are observed.

The contribution of the close-track effects discussed in Sect. 3.1 was estimated by varying the size of the wire-exclusion window (3 vs. 5 wires) in Method B. The resulting radius uncertainties did not exceed 0.2 fm; typically they were smaller than 0.1 fm.

The influence of possible impurities entering the charged pion samples was tested with stronger cuts on the quality parameters of particle identification (cf. Sect. 2), i.e. by a \(\pm 30~\text{ MeV }\) mass window around the most probable pion mass. About two third of the pairs survived this cut. No systematic differences w.r.t. to the full data sample are found within the statistical errors.

The effect of varying \(R_{\mathrm{inv}}\) in the Coulomb correction in Eq. (7) results in systematic uncertainties of \(\sim 0.01\)  fm; the effect of the finite size of the averaging intervals in Eq. (9) yields systematic uncertainties of even smaller size.

The uncertainty of the momentum resolution correction described in Sect. 3.4 appears to be an order of magnitude smaller than the absolute source radius shift, i.e. typical values of 0.01–0.03 fm are considered.

Another systematic uncertainty was estimated from studies of the forward-backward symmetry of the fit results w.r.t. midrapidity. Selecting the rapidity windows \(-\,0.35<y-y_\mathrm{cm}<0\) and \(0<y-y_\mathrm{cm}<0.35\) and taking similar (unbiased by the detector acceptance) transverse momentum intervals, typical systematic variations of the fit radii of 0.03–0.1 (0.2) fm for \(R_{\mathrm{inv}}\), \(R_{\mathrm{side}}\), \(R_{\mathrm{long}}\) (\(R_{\mathrm{out}}\)) are observed.

For the fit of the azimuthally integrated three-dimensional correlation function, the slight differences of the results when switching on/off the ‘out’-‘long’ component in the fit function (cf. Sect. 3.3) are taken as further systematic uncertainty. Typical values of these differences are 0.03–0.15 fm.

As an additional cross check, the stability of the results w.r.t. a reversed setting (for about 10% of the beam time) of the magnetic field has been investigated. Within the larger statistical errors, the results for \(\pi ^-\pi ^-\) (\(\pi ^+\pi ^+\)) in reversed field are found identical to the \(\pi ^-\pi ^-\) (\(\pi ^+\pi ^+\)) results in default field, although oppositely charged pions are subjected to different overall detector acceptances.

Finally, all systematic error contributions are added quadratically.

4 Results

4.1 Experimental correlation functions

The data are divided into centrality classes (cf. Sect. 2) and into classes of pair transverse momentum, \(p_{\mathrm{t},\,12}\), from which the pair transverse mass, \(m_\mathrm{t}=\sqrt{k_\mathrm{t}^2+m_{\pi }^2}\) with \(k_t=p_{\mathrm{t},\,12}/2\), is derived. As width of the underlying \(p_{\mathrm{t},\,12}\) bins 100 \(\text{ MeV }/c\) was selected. In case of azimuthally-sensitive HBT analyses (Sect. 4.4), the bin size of the pair angle w.r.t. the reaction plane, \(\varPhi =\phi _{12}-\phi _\mathrm{RP}\), is chosen to be \(\pi /4\).

A representative one-dimensional \(\pi ^{-}\pi ^{-}\) correlation function fitted with Eq. (5) is shown in Fig. 1. Note that we usually exclude from the fits pairs with \(q_{{\text {inv}}}< 6~\text{ MeV }/c\), since at very small relative momenta slight remnants of close-track effects could not completely be avoided when applying the correction Method B described in Sect. 3.1. This exclusion from the fit also renders a potential correction of the center of gravity of the first \(q_{{\text {inv}}}\) bin unnecessary, within which the Coulomb factor \(K_\mathrm{C}(q_{{\text {inv}}})\) exhibits a rather steep increase. (However, already within the 2nd bin function \(K_\mathrm{C}\) is fairly smooth. Moreover, since both the true and mixed yields naturally show different slopes at low \(q_{{\text {inv}}}\), a unique positioning of the corresponding bin centers required by Eq. (3) would be quite sophisticated.) For the quality of the three-dimensional fits we refer to [12].

Two-dimensional projections of the Coulomb-corrected three-dimensional correlation function (Eq. (7)) are presented in Fig. 2 for \(\vert \varPhi \vert <\pi /8\). The low-momentum enhancement due to the quantum-statistical (Bose–Einstein) effect is visible in all three directions. Comparing these in-plane correlations with the corresponding out-of-plane (\(\vert \varPhi -\pi /2\vert <\pi /8\)) correlations displayed in Fig. 3, a modification of the q range of the Bose–Einstein part of the correlation becomes visible. (Due to the permutability of particles 1 and 2, one of the q projections can be restricted to positive values.) Detailed results of the azimuthally-dependent correlation functions are presented in Sect. 4.4.

Fig. 1
figure 1

Upper panel: The distribution of the invariant relative momentum \(q_{{\text {inv}}}\) for \(\pi ^-\pi ^-\) pairs with transverse momentum of \(p_{\mathrm{t},\,12}=100{-}200~\text{ MeV }/c\) for central (0–10%) Au + Au collisions at 1.23\(A~\text{ GeV }\). The black (red) histogram displays the true (mixed) yield. The grey-shaded area represents the yield used for normalization. Lower panel: The one-dimensional \(\pi ^-\pi ^-\) correlation function as function of \(q_{{\text {inv}}}\). Red circles display the ratio of the true and mixed distributions, respectively. The green dashed curve represents the Coulomb correction function \(K_C\) as described in Sect. 3.2. The black squares correspond to the Coulomb-corrected correlation function. The red full (black dotted) curve shows the fit function (Eq. (5)) before (after) the Coulomb correction. The blue long-dashed curve gives the pure Bose–Einstein part (Eq. (6)) of the correlation function

Full size image
Fig. 2
figure 2

Two-dimensional projections of the Coulomb-corrected three-dimensional correlation function (Eq. (7)) in the ‘out’ vs. ‘side’ (upper left panel), ‘out’ vs. ‘long’ (upper right), and ‘side’ vs. ‘long’ (lower left) planes, respectively, for \(\pi ^-\pi ^-\) pairs with transverse momentum of \(p_{\mathrm{t},\,12}=100{-}400~\text{ MeV }/c\) and angles w.r.t. the reaction plane of \(\vert \phi _{12}-\phi _{\mathrm{RP}}\vert <\pi /8\) for centralities of 10–30%. The respective third direction is integrated over \(\pm \,\, 20~\text{ MeV }/c\)

Full size image
Fig. 3
figure 3

The same as Fig. 2, but for \(\vert \phi _{12}-\phi _{\mathrm{RP}}-\pi /2\vert <\pi /8\)

Full size image

4.2 Separation of central charge bias: construction of neutral pion radii

To quantify a potential source radius bias introduced by the Coulomb force the charged pions experience in the field of the charged fireball, we follow the ansatz used in Refs. [37, 38],

$$\begin{aligned} E(\varvec{p}_\mathrm{f}) = E(\varvec{p}_\mathrm{i}) \pm V_\mathrm{eff}(\varvec{r}_\mathrm{i}), \end{aligned}$$
(11)

where \(\varvec{p}_\mathrm{i}\) (\(\varvec{p}_\mathrm{f}\)) is the initial (final) momentum, E the corresponding total energy, and \(\varvec{r}_\mathrm{i}\) the initial position of the pion in the Coulomb potential \(V_{\mathrm{eff}}\) with positive (negative) sign for \(\pi ^+\) (\(\pi ^-\)). With

$$\begin{aligned} \frac{R_{\pi ^{\pm }\pi ^{\pm }}}{R_{\pi ^0\pi ^0}} \approx \frac{q_\mathrm{i}}{q_\mathrm{f}}&= \frac{|\varvec{p}_\mathrm{i}|}{|\varvec{p}_\mathrm{f}|} \nonumber \\&= \sqrt{ 1 \mp 2 \frac{V_{\mathrm{eff}}}{|\varvec{p}_\mathrm{f}|} \sqrt{1 + \frac{m_{\pi }^2}{\varvec{p}_\mathrm{f}^2}} + \frac{V^2_{\mathrm{eff}}}{\varvec{p}_\mathrm{f}^2}}, \end{aligned}$$
(12)

where \(q_\mathrm{i}\) (\(q_\mathrm{f})\) is the initial (final) relative momentum, and with \(V_{\mathrm{eff}} / k_\mathrm{t} \ll 1\), it turns out that the squared source radius for pairs of constructed neutral pions (denoted by \(\tilde{\pi }^0\tilde{\pi }^0\) in the following) is simply the arithmetic mean of the corresponding quantities of the charged pions,

$$\begin{aligned} R_{\tilde{\pi }^0\tilde{\pi }^0}^2 = \frac{1}{2} \big (R_{\pi ^+\pi ^+}^2 + R_{\pi ^-\pi ^-}^2\big ). \end{aligned}$$
(13)

Finally, the constructed \(\tilde{\pi }^{0}\tilde{\pi }^{0}\) correlation radii, discussed in the subsequent sections, are derived from cubic spline interpolations of the \(k_\mathrm{t}\) and \(A_\mathrm{part}\) dependences of the corresponding experimental \(\pi ^{-}\pi ^{-}\) and \(\pi ^{+}\pi ^{+}\) data.

4.3 Azimuthally-integrated HBT analysis

4.3.1 Central collisions

We start with the study of source radii for central (0–10%) collisions. Figure 4 shows the \(m_\mathrm{t}\) dependence of the one-dimensional (invariant) and three-dimensional source radii for \(\pi ^-\pi ^-\) (black squares) and \(\pi ^+\pi ^+\) (red circles) pairs. While for low transverse mass the Coulomb interaction with the fireball leads to an increase (a decrease) of the source size derived for negative (positive) pion pairs, at large transverse momentum the Coulomb effect apparently fades away. The effect is smallest for \(R_{{\text {out}}}\) and most pronounced for \(R_{{\text {side}}}\). Note that the charge splitting of the source radii as well as its difference in magnitude in out and side direction was early on predicted by Barz [39, 40] who investigated the combined effects of nuclear Coulomb field, radial flow, and opaqueness on two-pion correlations for a large collision system such as Au + Au in the \(1 A~\text{ GeV }\) energy regime. Earlier experimental works at the Bevalac employing a three-body Coulomb correction found the effect negligible for their studies of smaller systems [6,7,8]. We note that the extent of the source in direction of \(\varvec{k}_\mathrm{t}\) is potentially enlarged by a finite emission time duration, which is expected to last a few fm /c [1], cf. Sect. 4.4. It would be very interesting to study the charge-sign effect quantitatively with the help of dynamical models of the collision evolution. At present, such theoretical investigations concentrate on the higher energies where the charge splitting is hardly visible.

The parameter \(\lambda _{{\text {osl}}}\) derived from the fits with Eq. (7) appears to be rather independent of transverse mass and charge sign and decreases only slightly with increasing transverse mass, cf. lower right panel of Fig. 4. It fits well into a preliminary evolution with \(\sqrt{s_\mathrm{NN}}\) established previously [41], except the lowest E895 data point at 2\(A~\text{ GeV }\). In contrast, \(\lambda _{{\text {inv}}}\) resulting from the fits to the one-dimensional correlation function, exhibits a significant decrease with \(m_\mathrm{t}\) (cf. lower left panel), probably pointing to the fact that the one-dimensional fit function is not adequate.

Fig. 4
figure 4

Source radii as function of pair transverse mass, \(m_\mathrm{t}\), for central (0–10%) Au + Au collisions at 1.23\(A~\text{ GeV }\). The corresponding \(p_{\mathrm{t},\,12}\) range amounts to \(100-800~\text{ MeV }/c\). The upper left, upper right, center left, and center right panels display the invariant, ‘out’, ‘side’, and ‘long’ radii, respectively. The lower left and lower right panels show the corresponding \(\lambda \) parameters resulting from the fits to the one- and three-dimensional correlation functions, respectively. Black squares (red circles) are for pairs of negative (positive) pions. Blue dashed curves represent constructed radii (Eq. (13)) of neutral pion pairs. Error bars and hatched bands represent the statistical and systematic errors, respectively

Full size image
Fig. 5
figure 5

Source radii as function of the cube root of the number of participants for \(\pi ^{-}\pi ^{-}\) (upper left) and \(\pi ^{+}\pi ^{+}\) (upper right) pairs with transverse momenta of \(p_{\mathrm{t},\,12}=300{-}400~\text{ MeV }/c\). The lower left panel shows corresponding radii of constructed neutral pions as deduced from interpolation of the charged-pion data according to the procedure described in Sect. 4.2. Black circles, red triangles, green inverted triangles, and blue squares display the ‘invariant’, ‘out’, ‘side’, and ‘long’ radii, respectively. Dotted lines represent \(\chi ^2\) fits to the data with straight lines

Full size image

4.3.2 Centrality dependence

Figure 5 exhibits the \(A_\mathrm{part}\) dependence of source radii of \(\pi ^{-}\pi ^{-}\) pairs (upper left panel) for the transverse-momentum interval \(p_{\mathrm{t},\,12} = 300{-}400~\text{ MeV }/c\). The volume of the overlap region (i.e. the participant zone) of the colliding nuclei is expected to be proportional to \(A_{\mathrm{part}}\). Consequently, the involved length scales should be proportional to \(A_{\mathrm{part}}^{1/3}\), which therefore could provide a good scaling variable for emission source radii. Close to linear dependences of the pion source radii as function of \(A_{\mathrm{part}}^{1/3}\) are observed, as demonstrated by the straight-line fits to the data. The upper right panel of Fig. 5 shows the same for \(\pi ^{+}\pi ^{+}\) source radii, following similar linear dependences on \(A_{\mathrm{part}}^{1/3}\). Note that the \(R_{{\text {inv}}}(A_\mathrm{part})\) dependences (black dotted lines) of the charged pion pairs, when extrapolated down to the \(A_\mathrm{part}\) value (not shown) of the system Ar + KCl previously investigated by HADES [10], match the corresponding radii within uncertainties. This observation supports the idea that the variation of the participant volume via both, the choice of the size of the collision partners and the selection of a certain collision centrality, is equivalent. Following the recipe given in Sect. 4.2 to remove the Coulomb effect from the above dependences of \(\pi ^{-}\pi ^{-}\) and \(\pi ^{+}\pi ^{+}\) radii, we present the \(A_\mathrm{part}^{1/3}\) dependence of constructed \(\tilde{\pi }^{0}\tilde{\pi }^{0}\) source radii in the lower left panel of Fig. 5.

All results of the azimuthally integrated one-dimensional \((R_{{\text {inv}}},\,\lambda _{{\text {inv}}})\) and three-dimensional \((R_{{\text {out}}},\,R_{{\text {side}}},\,R_{{\text {long}}},\,\lambda _{{\text {osl}}})\) fits in dependence on collisions centrality and mean transverse momentum for \(\pi ^-\pi ^-\), \(\pi ^+\pi ^+\) and constructed \(\tilde{\pi }^{0}\tilde{\pi }^{0}\) pairs are summarized in Tables 2, 3 and 4, respectively.

Table 2 Source parameters resulting from fits with Eqs. (5) and (7) for \(\pi ^-\pi ^-\) pairs in dependence of centrality and average transverse momentum, \(\bar{k}_\mathrm{t}\). Values in the 1st (2nd) brackets represent the corresponding statistical (systematic) uncertainties in units of the last digit of the respective quantity
Full size table
Table 3 The same as Table 2, but for \(\pi ^+\pi ^+\) pairs
Full size table
Table 4 The same as Table 2, but for constructed \(\tilde{\pi }^0\tilde{\pi }^0\) pairs. In addition, the Coulomb potentials extracted from the \(\pi ^-\pi ^-\) and \(\pi ^+\pi ^+\) data according to the method described in Sect. 4.2 are given in columns 4 and 8–10
Full size table

Finally, the transverse-mass dependences of constructed \(\tilde{\pi }^{0}\tilde{\pi }^{0}\) radii for different centrality classes are summarized in Fig. 6. The data are fitted (dashed curves) with a function

$$\begin{aligned} R= R_{0} \left( \frac{m_\mathrm{t}}{m_\pi }\right) ^{\alpha }. \end{aligned}$$
(14)

The fit parameters, \(R_0=R(k_\mathrm{t}=0)\) and \(\alpha \), derived from the \(\chi ^2\) minimizations are summarized in Table 5. Here we emphasize that, while \(\alpha \) derived from a fit to \(R_{{\text {side}}}(m_\mathrm{t})\) of \(\tilde{\pi }^0\tilde{\pi }^0\) follows the expectation (i.e. \(\alpha =-0.5\)) of (3+1)D hydrodynamics [42], its absolute value is larger (smaller) by \(\sim 0.1\) when fitting the corresponding dependences of \(\pi ^{-}\pi ^{-}\) (\(\pi ^{+}\pi ^{+}\)).

Fig. 6
figure 6

Source radii as function of the transverse mass for pairs of constructed neutral pions as deduced from interpolations of the charged-pion data according to the procedure of Sect. 4.2. The upper left, upper right, lower left, and lower right panels display the ‘invariant’, ‘out’, ‘side’, and ‘long’ radii, respectively. Red circles, magenta squares, blue triangles, and green inverted triangles represent the centrality classes 0–10%, 10–20%, 20–30%, and 30–40%, respectively. Dashed curves are \(\chi ^2\) fits to the data with a function \(R= R_0 (m_\mathrm{t}/m_\pi )^\alpha \). The fit parameters \(R_0\) and \(\alpha \) are summarized in Table 5

Full size image

4.4 Azimuthally-sensitive HBT analysis

To estimate the geometrical quantities hidden in the azimuthal variation of the correlation function (Eq. (7)), we follow the recipe given bei Wiedemann and Heinz [43, 44] and perform a common fit to our data on squared radii from Eq. (10), using the entire set of fit equations (Eq. (2)) given in Ref. [45], which yields the elements of the spatial correlation tensor as ten angle-independent fit parameters:

$$\begin{aligned} S_{\mu \nu }=\langle \tilde{x}_{\mu }\,\tilde{x}_{\nu }\rangle ,\,\,\tilde{x}_{\mu }=x_{\mu }-\langle x_{\mu } \rangle ,\quad (\mu ,\nu =0,\,1,\,2,\,3). \end{aligned}$$
(15)

Here, \(x_0=t\) is the time component, \(x_1=x\) is parallel to the impact parameter \(\varvec{b}\), and \(x_3=z\) points in beam direction. The Cartesian coordinate system is completed with direction \(x_2=y\) being perpendicular to the reaction plane formed by x and z. The brackets \(\langle \,\rangle \) indicate an average over the emission source.

A typical fit result is displayed in Fig. 7 for \(\pi ^{-}\pi ^{-}\) pairs with average transverse momentum of \(\bar{k}_\mathrm{t}=170~\text{ MeV }/c\) and for centralities of 10–30%. It delivers the uncorrected matrix S (in units of \(\hbox {fm}^2\), with statistical errors attached):

$$\begin{aligned}&S^{\mathrm{meas}} \nonumber \\&\quad = \begin{pmatrix} ~10.12 \pm 0.64 &{}\quad -0.77 \pm 0.21 &{}\quad -0.15 \pm 0.23 &{}\quad ~~~0.09\pm 0.14 \\ -0.77 \pm 0.21 &{}\quad ~21.60 \pm 0.28 &{}\quad -0.07 \pm 0.19 &{}\quad ~~~2.10\pm 0.11 \\ -0.15 \pm 0.23 &{}\quad -0.07 \pm 0.19 &{}\quad ~28.87 \pm 0.34 &{}\quad -0.05\pm 0.11 \\ ~~~0.09 \pm 0.14 &{}\quad 2.10 \pm 0.11 &{}\quad -0.05 \pm 0.11 &{}\quad ~16.12\pm 0.13 \end{pmatrix}. \end{aligned}$$
(16)

As expected from symmetry considerations [46], only the diagonal elements and \(S_{13}\) differ significantly from zero, yielding the following six non-vanishing squared radii

$$\begin{aligned} R^2_{\mathrm{out}} =\,&\,\frac{1}{2}(S_{11}+S_{22})+\frac{1}{2}(S_{22}-S_{11})\cos {(2\varPhi )}\nonumber \\ \,&\,\quad =+\langle \beta ^2_t\rangle S_{00}, \nonumber \\ R^2_{\mathrm{side}}=\,&\,\frac{1}{2}(S_{11}+S_{22})+\frac{1}{2}(S_{22}-S_{11})\cos {(2\varPhi )}, \nonumber \\ R^2_{\mathrm{long}}=\,&\,S_{33}+\langle \beta ^2_l\rangle S_{00}, \nonumber \\ R^2_{\mathrm{outside}}=\,&\,\frac{1}{2}(S_{22}-S_{11})\sin {(2\varPhi )}, \nonumber \\ R^2_{\mathrm{outlong}}=\,&\,S_{13}\cos {(\varPhi )}, \nonumber \\ R^2_{\mathrm{sidelong}}=\,&\,-S_{13}\sin {(\varPhi )}, \end{aligned}$$
(17)

where \(\beta _\mathrm{t}\) and \(\beta _\mathrm{l}\) are the pair velocities in transverse and longitudinal direction, respectively.

The correction of both, the finite reaction-plane resolution and the finite azimuthal bin width, is performed following the method described in refs. [24, 35]:

$$\begin{aligned} R^{2,\mathrm {corr}}_{\mathrm{i,n}} = R^{2,\mathrm {meas}}_{\mathrm{i,n}} \, \dfrac{n \, \varDelta /2}{F_{\mathrm{n}} \sin (n \, \varDelta /2)}, \end{aligned}$$
(18)

where \(R^{2,\mathrm {corr}}_{\mathrm{i,n}}\) (\(i=\) ‘out’, ‘side’, ‘long’, ‘outside’, ‘outlong’, ‘sidelong’,  \(n=0,\,1,\,2\)) are the underlying (“true”) Fourier coefficients [35], \(\varDelta =\pi /4\) is the present \(\varPhi \) interval, and the quantity \(F_{\mathrm{n}}\) represents the n-th event-plane resolution.Footnote 2 The values of \(F_1\) and \(F_2\) for the centrality classes investigated in the present analysis are summarized in Table 1. For the detailed description of the event-plane reconstruction and the centrality dependence of the resolution parameters \(F_\mathrm{n}\) we refer to a forthcoming HADES paper on the n-th (\(n \le 4\)) collective flow observables of protons and light nuclei produced in Au + Au collisions at 1.23\(A~\text{ GeV }\) [47].

Fig. 7
figure 7

Squared radii as function of the pair azimuthal angle relative to the reaction plane, \(\varPhi =\phi _{12}-\phi _\mathrm{RP}\), for \(\pi ^-\pi ^-\) pairs with transverse momenta of \(p_{\mathrm{t},\,12}=300{-}400~\text{ MeV }/c\) and for centralities of 10–30%. The left column from top to bottom shows the ‘out’, ‘side’, and ‘long’ radii. The right column gives the ‘outside’, ‘outlong’, and ‘sidelong’ components. Error bars include only statistical uncertainties. The full curves represent a global fit with Eq. (2) of Ref. [45] to the experimental data points

Full size image
Fig. 8
figure 8

The spatial principal axes (left column),the tilt angle w.r.t. the beam axis in the reaction plane (Eq. (20), top right), the xy-eccentricity (Eq. (22, center right), and the zy-eccentricity (Eq. (23, bottom right) of the Gaussian emission ellipsoid of \(\pi ^{-}\pi ^{-}\) pairs as function of pair transverse momentum for different centrality classes (cf. legend). Error bars include only statistical uncertainties. Dotted lines represent the corresponding initial eccentricity as derived from Glauber simulations

Full size image
Fig. 9
figure 9

The same as Fig. 8, but for \(\pi ^{-}\pi ^{-}\) (black squares) and \(\pi ^{+}\pi ^{+}\) (red circles) pairs in comparison for medium centralities of 10–30%

Full size image
Fig. 10
figure 10

The volume of the region of homogeneity (Eq. (24), full symbols) as derived from the azimuthally-sensitive analysis and the approximate volume (Eq. (25), open symbols) following from the azimuthally-integrated analysis of \(\pi ^{-}\pi ^{-}\) (squares) and \(\pi ^{+}\pi ^{+}\) (circles) pairs as function of pair transverse momentum for different centralities (cf. legend)

Full size image

The matrix S (Eq. (16)) after the correction (Eq. (18)) reads

$$\begin{aligned}&S^{\mathrm{corr}} \nonumber \\&\quad =\begin{pmatrix} ~10.12\pm 0.64 &{}\quad -\,0.91\pm 0.25 &{}\quad -\,0.18\pm 0.27 &{} \quad 0.09\pm 0.14 \\ -\,0.91\pm 0.25 &{} \quad 18.60\pm 0.40 &{} \quad -\,0.13\pm 0.36 &{} \quad 2.49\pm 0.13 \\ -\,0.18\pm 0.27 &{} \quad -\,0.13\pm 0.36 &{} \quad 31.87\pm 0.44 &{} \quad -\,0.06\pm 0.14 \\ ~~~0.09\pm 0.14 &{} \quad 2.49\pm 0.13 &{} \quad -\,0.06\pm 0.14 &{}\quad 16.12\pm 0.13 \end{pmatrix}. \end{aligned}$$
(19)

From this matrix the spatial tilt angle in the reaction plane can be calculated as

$$\begin{aligned} \theta _\mathrm{s}=\frac{1}{2}\tan ^{-1}\left( \frac{2 S^\mathrm {\,corr}_{13}}{S^\mathrm {\,corr}_{33}-S^\mathrm {\,corr}_{11}}\right) =(-32\pm 2)^{\circ }. \end{aligned}$$
(20)

Rotating \(S^{\mathrm{corr}}\) by the angle \(-\theta _\mathrm{s}\) around the y axis, i.e. applying the corresponding rotation matrix \(G_y(\theta _\mathrm{s})\), yields a diagonal tensor

$$\begin{aligned}&S^{\mathrm{diag}} = G^{\dag }_y(\theta _\mathrm{s})\cdot S^{\mathrm{corr}}\cdot G_y(\theta _\mathrm{s}) \nonumber \\&\quad =\begin{pmatrix} ~10.12\pm 0.64 &{}\quad -0.73\pm 0.23 &{}\quad -0.18\pm 0.27 &{}\quad 0.56\pm 0.18 \\ -0.73\pm 0.23 &{}\quad 20.15\pm 0.35 &{} \quad -0.15\pm 0.31 &{}\quad 0.00\pm 0.23 \\ -0.18\pm 0.27 &{} \quad -0.15\pm 0.31 &{}\quad ~31.87\pm 0.44 &{} \quad ~~~0.02\pm 0.22 \\ ~~~0.56\pm 0.18 &{} \quad ~~~0.00\pm 0.23 &{}\quad ~~~0.02\pm 0.22 &{}\quad ~14.58\pm 0.24 \end{pmatrix} \end{aligned}$$
(21)

whose eigenvalues are the temporal and geometrical variances \(\sigma ^2_\mathrm{t}, \sigma ^2_\mathrm{x}, \sigma ^2_\mathrm{y}, \sigma ^2_\mathrm{z}\). The geometrical variances of \(\pi ^{-}\pi ^{-}\) emission as function of transverse-momentum are displayed in the left column of Fig. 8 for all centralities. Figure 8 also shows the transverse-momentum dependences of the tilt angle, \(\theta _\mathrm{s}\), the xy-eccentricity,

$$\begin{aligned} \varepsilon _{xy}=\frac{\sigma _\mathrm{y}^2-\sigma _\mathrm{x}^2}{\sigma _\mathrm{x}^2+\sigma _\mathrm{y}^2}, \end{aligned}$$
(22)

and the zy-eccentricity,

$$\begin{aligned} \varepsilon _{zy}=\frac{\sigma _\mathrm{y}^2-\sigma _\mathrm{z}^2}{\sigma _\mathrm{y}^2+\sigma _\mathrm{z}^2}. \end{aligned}$$
(23)

Note that we included an intermediate centrality class of 25–35% to verify the strong centrality dependence of the tilt angle (upper right panel of Fig. 8 and discussion below). For all transverse momenta and all centralities, the deduced eccentricities represent an almond shape in the plane perpendicular to the beam direction. The xy-eccentricity becomes smallest (almost circular shape) for most central collisions. Almost no centrality dependence appears for the zy-eccentricity. The range of variation with transverse momentum is larger for the zy-eccentricity as compared to the one in the xy-plane. For large momenta, \(\theta _\mathrm{s}\) tends to vanish, provided that the 30% most central event classes (i.e. impact parameters not considerably larger than 50% of the maximum one) are selected, and the final \(\pi ^{-}\pi ^{-}\)xy-eccentricity recovers the corresponding initial eccentricity. No charge-sign difference appears in the transverse-momentum dependence of both, the tilt angle and the eccentricities, while the spatial principal axes differ, as demonstrated in Fig. 9 for medium centralities of 10–30%.

Due to the freedom in the coordinate-system definition, the tilt angle, defined as the angle between the z coordinate (directed along the shortest principal axis in our case) and the beam direction, can be changed from \(\theta _\mathrm{s}\) to \(\theta _\mathrm{s}-90^{\circ }\), while \(\sigma _x\) and \(\sigma _z\) (and accordingly \(\varepsilon _{xy}\) and \(\varepsilon _{xy}\)) interchange. The arrangement of our data (cf. Fig. 8) was chosen such that both dependences, on transverse momentum and on centrality, show smooth trends. Thus, smaller \(\vert \theta _\mathrm{s}\vert \) for more central collisions are ensured, as one would expect from collision geometry. However, at more peripheral collisions and high values of \(p_{\mathrm{t},12}\) a different configuration is conceivably, i.e. requiring \(\theta _\mathrm{s} = 0\) at high transverse momenta (relevant here only for the data points at \(p_{\mathrm{t},12} = 440 \, \mathrm {MeV}/c\) and centrality 25–35%). Within our statistical and systematic uncertainties we are not able to find a final decision, which arrangement is the better one. (Note also that, near \(p_{\mathrm{t},12} = 400 \, \mathrm {MeV}/c\) and for centralities beyond 20%, \(\sigma _x\) and \(\sigma _z\) are about the same size, which disturbs the picture of a well defined tilt angle, eventually represented by large uncertainties of this quantity.)

The volume of the region of homogeneity,

$$\begin{aligned} V_\mathrm{fo}=(2\pi )^{3/2}\sigma _{\mathrm{x}} \sigma _{\mathrm{y}} \sigma _{\mathrm{z}}, \end{aligned}$$
(24)

as derived from the spatial principal axes of the Gaussian emission ellipsoid is presented in Fig. 10 (full symbols) as function of transverse momentum for all centrality classes. For comparison, also the results for the approximate volume,

$$\begin{aligned} V_\mathrm{fo}'=(2\pi )^{3/2}R^2_{\mathrm{side}}R_{\mathrm{long}}, \end{aligned}$$
(25)

following from the azimuthally-integrated analysis (open symbols) is shown. Differences are significant only at small transverse momenta, where the tilt angle vanishes.

Fig. 11
figure 11

The emission-ellipsoid tilt angle, \(\theta _\mathrm{s}\), in the reaction plane for \(\pi ^{-}\pi ^{-}\) (full symbols) and \(\pi ^{+}\pi ^{+}\) (open symbols) pairs as function of the mean impact parameter derived from Glauber MC calculations [21] for different pair transverse momentum classes. Error bars include only statistical uncertainties

Full size image
Fig. 12
figure 12

The \(\pi ^{-}\pi ^{-}\) emission ellipsoid xy-eccentricity (Eq. (22)) as function of the initial nucleon eccentricity derived from Glauber simulations [22] for different pair transverse momentum classes. Error bars include only statistical uncertainties. The dashed line indicates \(\varepsilon _{\mathrm{final}} = \varepsilon _{\mathrm{initial}}\)

Full size image
Fig. 13
figure 13

The temporal components of the spatial correlation tensor S (Eq. (15)) for \(\pi ^{-}\pi ^{-}\) pairs as function of transverse momentum for different centrality classes (cf. legend). The upper left, upper right, lower left, and lower right panels display \(\sigma _\mathrm{t}^2=S_{00}\), \(S_{01}\), \(S_{02}\), and \(S_{03}\), respectively. Error bars include only statistical uncertainties

Full size image

To study the centrality dependence of the tilt angle in more detail, Fig. 11 displays \(\theta _\mathrm{s}\) for \(\pi ^{-}\pi ^{-}\) (full symbols) and \(\pi ^{+}\pi ^{+}\) (open symbols) pairs as function of the impact parameter, b (cf. Table 1), for different average transverse momentum values. While for lower momenta the magnitude of \(\theta _\mathrm{s}\) is proportional to b, this dependence gets weaker with increasing pair transverse momentum until the tilt is close to zero for the highest momentum classes. No charge-sign difference is observed in the centrality dependence of \(\theta _\mathrm{s}\).

Figure 12 explicitly relates the xy-eccentricity for \(\pi ^-\pi ^-\) pairs to the initial eccentricity relative to the participant plane (cf. Eq. (21) of Ref. [41]), as derived from Glauber simulations [22]. For large transverse momenta the source eccentricity derived from the present identical-pion HBT analysis recovers the initial (nucleon) eccentricity. This observation is essentially different from the picture at higher energies, where the \(k_\mathrm{t}\) integrated freeze-out eccentricity is found considerably below the corresponding initial one [41, 48, 49]. The trend of increasing freeze-out eccentricity with increasing transverse momentum is also found by ALICE at LHC, though on a much lower absolute scale [49]. The model of (3+1)D hydrodynamics [42] slightly underestimates the final source eccentricity at LHC. It would be very desirable to discuss in separate future works the present observations at SIS18 energies in the framework of contemporary hydro or transport models.

Figure 13 shows the temporal components of the spatial correlation tensor (cf. Eq. (2) of Ref. [45]). For the emission duration, \(\sigma _\mathrm{t}^2=S_{00}\), the fit yields values significantly deviating from zero, especially for low transverse momenta and for all considered centrality classes except the most central one. The values for the mixed elements, however, are much smaller and are mostly consistent with zero, with the possible exception of pion pairs with small transverse momenta in central collisions. As already mentioned above, the disappearance of \(S_{01}\), \(S_{02}\) and \(S_{03}\) is expected due to symmetry reasons [46]. However, one should keep in mind that the symmetry of the rapidity distribution w.r.t. midrapidity is not perfect in the present fixed-target experiment, even though the rapidity interval was rather restricted for this analysis (\(\vert y - y_{cm} \vert < 0.35\)).

Finally, the geometrical and temporal variances and the tilt angle of the \(\pi ^{-}\pi ^{-}\) (\(\pi ^{+}\pi ^{+}\)) emission source in dependence of centrality and transverse momentum are summarized in Table 6 (7).

Fig. 14
figure 14

The low-energy excitation function of the principal axes (left column), the tilt angle in the reaction plane (Eq. (20), top right), the eccentricities \(\varepsilon _{xy}\) (Eq. (22), center right, open symbols) and \(\varepsilon _{zy}\) (Eq. (23), center right, filled symbols), and the emission duration \(\sigma _\mathrm{t}=\sqrt{S_{00}}\) (Eq. (17), bottom right) as derived from the Gaussian emission ellipsoid of \(\pi ^{-}\pi ^{-}\) pairs. Stars are HADES data for centralities of 10–30% and average transverse momentum of \(\bar{k}_\mathrm{t}=110~\text{ MeV }/c\). The circles are corresponding E895 data for slightly different centralities taken at beam energies of 2, 4, and 6 \(A~\text{ GeV }\) [45]. Error bars include only statistical uncertainties

Full size image

4.5 Comparison with other experiments: excitation functions of source parameters

The energy dependences below \(\sqrt{s_{\mathrm{NN}}}=4~\hbox {GeV}\) of the \(\pi ^{-}\pi ^{-}\) emission ellipsoid principal axes, the tilt angle w.r.t. the beam axis in the reaction plane, the eccentricities, and the emission duration for medium centralities and average transverse momentum of \(\bar{k}_\mathrm{t}=110~\text{ MeV }/c\) are shown in Fig. 14. The displayed HADES/SIS18 data follow the trends of the E895/AGS data [45].

Fig. 15
figure 15

Excitation function of the final \(\pi ^{-}\pi ^{-}\) eccentricities \(\varepsilon _{xy}\) (open symbols) and \(\varepsilon _{zy}\) (filled symbols) for mid-central (10–30%) Au + Au, Pb + Au, or Pb + Pb collisions. Stars are HADES data for three different transverse momentum intervals with average values of \(\bar{k}_\mathrm{t}= 110\), 230, and 310\(~\text{ MeV }/c\). The circles, diamond, triangles, and square are corresponding data by E895 at AGS [45], CERES at SPS [51], STAR at RHIC [41], and ALICE at LHC [49], respectively. Error bars include only statistical uncertainties. Short (long) dashed curves connecting the filled (open) symbols are to guide the eyes. The arrow (to be compared to the open symbols) indicates the initial eccentricity derived from Glauber simulations [22]

Full size image
Fig. 16
figure 16

Excitation function of the squared emission duration, the latter one either directly taken from the corresponding fit parameter of the spatial correlation tensor, \(\sigma _\mathrm{t}^2=S_{00}\) (cf. Eq. (17)), or derived from the difference of the ‘out’ and ‘side’ radii, \((\varDelta \tau )^2 = (R^2_\mathrm{out}-R^2_\mathrm{side})/\langle \beta _t^2\rangle \), for mid-central (10–30%) Au + Au collisions. Stars are HADES \(\pi ^{-}\pi ^{-}\) data of \(\sigma _\mathrm{t}^2\) for three different transverse momentum regions with average values of \(\bar{k}_\mathrm{t}= 110\), 230, and 310\(~\text{ MeV }/c\). The circles are corresponding \(\sigma _\mathrm{t}^2\) data by E895 at AGS [45]. The diamond, triangles, and square represent data of \((\varDelta \tau )^2\) by CERES at SPS [51], STAR at RHIC [41], and ALICE at LHC [49], respectively. Error bars include only statistical uncertainties

Full size image
Fig. 17
figure 17

Excitation function of the source radii \(R_{\mathrm{out}}\) (upper panel), \(R_{\mathrm{side}}\) (central panel), and \(R_{\mathrm{long}}\) (lower panel) for the azimuthally-integrated correlation function of pairs of identical pions with transverse mass of \(m_{t}=330~\text{ MeV }\) in central (0–10%) collisions of Au + Au or Pb + Pb. Squares represent data by ALICE at LHC (\(\pi ^+\pi ^+\)) [52], full triangles STAR at RHIC (\(\pi ^-\pi ^{-}+\pi ^+\pi ^+\)) [41], diamonds are for CERES at SPS (\(\pi ^-\pi ^{-}+\pi ^+\pi ^+\)) [53], open triangles are for NA49 at SPS (\(\pi ^-\pi ^-\)) [54], open circles are \(\pi ^-\pi ^-\) data by E895 at AGS [34, 55], and open (full) crosses involve \(\pi ^-\pi ^-\) (\(\pi ^+\pi ^+\)) data of E866 at AGS [56], respectively. The present data of HADES at SIS18 for pairs of \(\pi ^-\pi ^-\) (\(\pi ^+\pi ^+\)) are given as open (full) stars. Statistical errors are displayed as error bars; if not visible, they are smaller than the corresponding symbols. Note the suppressed zero on the ordinate

Full size image

The excitation function of eccentricities over a wider range of collision energies is presented in Fig. 15. While for the HADES and E895 data the exact equations (22) and (23) are applied (cp. Fig. 14), for \(\sqrt{s_{\mathrm{NN}}}>4~\hbox {GeV}\), where no \(\pi \pi \) emission source data after principal-axis transformation of the spatial correlation tensor are available, the approximations \(\varepsilon _\mathrm{xy}\approx 2R^2_{\mathrm{side,\,2}}/R^2_{\mathrm{side,\,0}}\) and \(\varepsilon _{zy}\approx 1- 2R^2_{\mathrm{long,\,0}}/(R^2_{\mathrm{side,\,0}}+2R^2_{\mathrm{side,\,2}}+R^2_{\mathrm{long,\,0}})\) are used. Here, the quantities \(R^2_{\mathrm{side,\,0}}\) (\(R^2_{\mathrm{long,\,0}}\)) and \(R^2_{\mathrm{side,\,2}}\) are the zeroth- and 2nd-order Fourier coefficients of the azimuthal-angle parameterization of the ‘side’ (‘long’) radius, respectively, where \(2 R^2_{\mathrm{side,\,2}}\) corresponds to the oscillation amplitude of both the squared ‘out’ and ‘side’ radii (cf. Eq. (17) and Fig. 7). These approximations are justified by the small tilt angles found at high collision energies. We observe a remarkable increase of both types of eccentricity at low \(\sqrt{s_\mathrm{NN}}\). It would be very desirable to study the strong increase towards low energy of both eccentricity and tilt angle in the framework of contemporary models. First investigations of both observables by means of several theoretical transport calculations qualitatively reproduced the experimental findings [50]; the quantitative reproduction of their dependences on collision energy, centrality and transverse momentum should allow for further discrimination among the models.

The excitation function of the squared emission duration which is either taken directly from the fit parameter \(S_{00}\), i.e. the temporal variance of the \(\pi \pi \) emission (cf. Eq. (17)), or from the difference of the squared ‘out’ and ‘side’ radii, \((\varDelta \tau )^2=(R_\mathrm{out}^2-R_\mathrm{side}^2)/\langle \beta ^2_\mathrm{t}\rangle \), is shown in Fig. 16. For all data of the various experiments, the emission duration appears small, i.e. typically a few fm/c. With increasing transverse momentum, the emission duration deduced from the HADES data is found to decrease, virtually vanishing for large momenta. Taking the data points at \(k_\mathrm{t}=310~\text{ MeV }/c\), there is a clear energy dependence: A rise towards \(\sqrt{s_\mathrm{NN}}\sim 10{-}20~\hbox {GeV}\) and then a slow drop towards LHC.

Finally, we return to a few results of the azimuthally-integrated HBT analysis, concentrating on central collisions (0–10%). The excitation functions of \(R_{\mathrm{out}}\), \(R_{\mathrm{side}}\), and \(R_{\mathrm{long}}\) for pion pairs produced in most central events are displayed in Fig. 17. All shown radius parameters have been obtained by interpolating the existing measured data points to the same transverse mass of \(m_\mathrm{t}=330~\text{ MeV }\) (\(k_\mathrm{t}=300~\text{ MeV }/c\)) at which the charge differences of the source radii almost vanish (cf. Fig. 4). The statistical errors are properly propagated and quadratically added to the differences between linear and cubic-spline interpolations. Extrapolations were not necessary at this \(m_\mathrm{t}\) value. \(R_{\mathrm{out}}\) and \(R_{\mathrm{side}}\) vary not more than 40% over three orders of magnitude in center-of-mass energy. Only \(R_{\mathrm{long}}\) exhibits a steady increase by about a factor of two when going in energy from SIS18 via AGS, SPS, RHIC to LHC. Note that in the excitation functions shown in Ref. [41] not all, particularly AGS, data points were properly corrected for their \(k_\mathrm{t}\) dependence.

The excitation function of \(R_{\mathrm{out}}^2-R_{\mathrm{side}}^2\) for an average transverse momentum of the pion pairs of \(k_\mathrm{t}=300~\text{ MeV }/c\) in central collisions is shown in Fig. 18. Almost all other measurements below 10 GeV are characterized by large errors and scatter sizeably. The new HADES data show that the difference of the source parameters in the transverse plane almost vanishes at low collision energies. Since this quantity is related to the emission duration via Eq. (17), one would conclude that in the 1\(A~\text{ GeV }\) energy region the observed pions are emitted into free space during a short time span of less than a few fm/c (cp. also Figs. 14 (bottom right) and 16 displaying similar data divided by \(\langle \beta ^2_\mathrm{t}\rangle \) for centralities of 0–10% and 10–30%, respectively, but for different transverse momenta). However, also the opaqueness of the source affects \(R_{\mathrm{out}}^2-R_{\mathrm{side}}^2\) which could cause it to become negative, thus compensating the positive contribution of the emission duration [40]. We also emphasize that, with increasing available energy, this quantity reaches a local maximum at \(\sqrt{s_{\mathrm{NN}}}\sim 20{-}30\) GeV and afterwards decreases towards zero at LHC energies.

Fig. 18
figure 18

The same as Fig. 17, but for the derived quantity \(R_{\mathrm{out}}^2-R_{\mathrm{side}}^2\)

Full size image
Fig. 19
figure 19

The same as Fig. 17, but for the derived approximate volume of the region of homogeneity, Eq. (25)

Full size image
Table 5 Parameters from the fits with Eq. (14) to the data of Fig. 6
Full size table
Table 6 The geometrical and temporal variances (columns 3 to 6) and the tilt angle (column 7) of the \(\pi ^-\pi ^-\) emission source in dependence of centrality and average transverse momentum, \(\bar{k}_\mathrm{t}\). Values in the 1st (2nd) brackets represent the corresponding statistical (systematic) uncertainties in units of the last digit of the respective quantity
Full size table
Table 7 The same as Table 6, but for \(\pi ^+\pi ^+\) pairs
Full size table

The excitation function of the approximate volume of the region of homogeneity, Eq. (25), for central collisions is given in Fig. 19. Here, we chose this approximation, in contrast to Eq. (24), for the sake of comparability with other experiments. Note that this definition of a three-dimensional Gaussian volume does not incorporate \(R_{{\text {out}}}\) since this length is potentially extended due to a finite value of the aforementioned emission duration. For the large transverse momenta selected, the differences of the azimuthally-integrated and azimuthally-sensitive analyses as well as the charge-sign splitting largely vanish, cp. Figure 10 exhibiting the strong transverse-momentum dependence of the volume, especially for \(\pi ^{-}\pi ^{-}\) pairs. From the above HADES data, we estimate a volume of about \(850\,\hbox {fm}^3\) for constructed \(\tilde{\pi }^0\tilde{\pi }^0\) pairs. This volume of homogeneity steadily increases with energy, but appears merely a factor four larger at LHC. Extrapolating the volume to \(k_\mathrm{t}=0\) yields a value of about \(3900\,\hbox {fm}^3\).

Similar excitation functions at other transverse momenta or centralities (with proper interpolation/extrapolation of the transverse-momentum and centrality dependences) can be derived from the source parameters summarized in Tables 2, 3, 6, and 7.

5 Summary and outlook

We presented high-statistics \(\pi ^{-}\pi ^{-}\) and \(\pi ^{+}\pi ^{+}\) HBT data for Au + Au collisions at 1.23\(A~\text{ GeV }\). The three-dimensional Gaussian emission source is studied in dependence on transverse momentum and collision centrality. It is found to follow the trends observed at higher collision energies, extending the corresponding excitation functions towards very low energies. A surprisingly small variation of the space-time extent of the pion emission source over three orders of magnitude of \(\sqrt{s_{\mathrm{NN}}}\) is observed. All source radii increase almost linearly with the number of participants, irrespective of transverse momentum. Substantial differences of the source radii for pairs of negatively and positively charged pions, especially at low transverse momenta, are found, an effect hardly visible at higher collision energies. Correcting for this Coulomb effect, we found for all transverse momenta and centralities a stable hierarchy of the three widths of the emission ellipsoid, revealing an oblate shape in coordinate space with the largest extent oriented perpendicular to the reaction-plane. The corresponding eccentricity is found to recover the initial geometrical eccentricity of the nucleons if sufficiently large transverse momenta of the pions are considered. For low momenta, the shape approaches a circular one. Selecting collisions of the 30% most central event classes, the tilt angle in the reaction plane, \(\theta _\mathrm{s}\), is found to approach zero at high \(p_\mathrm{t}\). For low transverse momenta, \(\vert \theta _\mathrm{s}\vert \) increases and its \(p_{\mathrm{t}}\) dependence is stronger for larger impact parameters. The centrality dependences of the tilt angle and of the eccentricity are found to be independent of pion charge and transverse momentum. Both quantities fit well into a corresponding low-energy excitation function for semi-central collisions [45].

Future femtoscopic studies of other system-energy combinations in the ‘low-energy domain’, e.g. performed with HADES + CBM at FAIR/GSI [57] and MPD at NICA in Dubna [58] or with the STAR fixed-target program [59], may focus also on correlations of other particles (pp, \(\hbox {p}\varLambda \), \(\varLambda \varLambda \), etc.), allowing for deeper insight into their strong interaction [60].