1 Introduction and motivation

This paper presents experimental results on two-particle correlations in pseudorapidity and azimuthal angle of charged particles produced in inelastic p + p interactions at 20, 31, 40, 80, and 158 GeV/c. The measurements were performed by themulti-purpose NA61/SHINE [1, 2] experiment at the CERN Super Proton Synchrotron (SPS). They are part of the strong interactions programme devoted to the study of the properties of the onset of deconfinement and search for the critical point of strongly interacting matter. Within this program a two-dimensional scan in collision energy and size of colliding nuclei is in progress. Data on p + p, Be + Be, and Ar + Sc collisions were already recorded and data on previous version: p + Pb and Xe + La collisions will be registered within the coming years. The expected signal of a critical point is a non-monotonic dependence of various fluctuation measures in such a scan; for a recent review see Ref. [3].

Apart from looking for QGP signatures, it is of interest to study specific physical phenomena that happen during and after the collision. The two-particle correlation analysis in pseudorapidity (\(\eta \)) and azimuthal angle (\(\phi \)) allows one to disentangle different correlation sources which can be directly connected with phenomena like jets, collective flow, resonance decays, quantum statistics effects, conservation laws, etc.

Measurements of two-particle correlations in pseudorapidity and azimuthal angle were first published by the ACM collaboration at the Intersecting Storage Rings (ISR) [4]. Two- and three-body decays of resonances (\(\eta \), \(\rho ^0\), \(\omega \)) were found to provide the dominant contributions. Two structures were discovered: an enhancement near \(\Delta \phi = \pi \) (away-side) explained by the two-body decay scenario and another enhancement at \(\Delta \phi \approx 0\) together with an azimuthal ridge (centred at \(\Delta \eta \approx 0\)) consistent with three-body decays.Footnote 1 These features were confirmed at the higher collision energies of RHIC by the PHOBOS [5] collaboration.

At RHIC and the LHC parton scattering processes become important. In addition to high transverse momentum jets, studies of \(\Delta \eta \Delta \phi \) correlations in p + p interactions as well as in collisions of heavy nuclei  [6,7,8,9] found prominent structures explained as arising from the production of minijets, producing a large correlation peak at small opening angles \((\Delta \eta ,\Delta \phi ) \approx (0,0)\) and a broad structure along \(\Delta \eta \) at \(\Delta \phi \approx \pi \) (also referred to as away-side ridge). At SPS energies, however, contributions from hard and semi-hard scattering processes are expected to be much smaller.

Collective flow effects also provide significant contributions to two-particle correlations. In non-central nucleus–nucleus collisions they appear as a modulation in the distribution of \(\Delta \phi \) both at SPS [10] and higher collider energies [11]. Recently a similar behaviour was observed in high-multiplicity p + Pb collisions [12] and even p + p reactions [13] at the LHC. The origin of the effect in p + p reactions is not yet understood. At SPS energies sufficiently high particle multiplicities cannot be produced in p + p collisions. Thus the effects of collective flow will not be discussed further.

In this paper the pseudorapidity variable \(\eta \) is calculated as \(\eta = -\text {ln} \text { tan}(\varTheta /2)\) where \(\text { tan}(\varTheta ) = p_T/p_L\) with \(p_T\) the transverse (xy) and \(p_L\) the longitudinal (z) component of the particle momentum in the collision centre-of-mass system. In Lorentz transformation of \(p_L\) (from NA61/SHINE originally used laboratory system to centre-of-mass system) the pion mass was assumed for all particles. The azimuthal angle \(\phi \) is the angle between the transverse momentum vector and the horizontal (x) axis.

The paper is organised as follows. The correlation function is defined in Sect. 2. The experimental setup is presented in Sect. 3. Data processing and simulation are described in Sect. 4. Data selection and analysis are discussed in Sect. 5. Results are presented in Sect. 6 and compared to model calculations in Sect. 7. A summary in Sect. 8 closes the paper.

2 Two-particle correlations in pseudorapidity and azimuthal angle

Correlations studied in this paper were calculated as a function of the difference in pseudorapidity (\(\eta \)) and azimuthal angle (\(\phi \)) between two particles produced in the same event:

$$\begin{aligned} \Delta \eta = |{\eta }_1 - {\eta }_2|, \quad \Delta \phi = |{\phi }_1 - {\phi }_2|. \end{aligned}$$
(1)

The correlation function \(C(\Delta \eta ,\Delta \phi )\) is defined and calculated as

$$\begin{aligned} C(\Delta \eta ,\Delta \phi )= \frac{N_{\text {mixed}}^{\text {pairs}}}{N_{\text {data}}^{\text {pairs}}} \frac{D(\Delta \eta ,\Delta \phi )}{M(\Delta \eta ,\Delta \phi )}, \end{aligned}$$
(2)

where

$$\begin{aligned} D(\Delta \eta ,\Delta \phi )=\frac{\mathrm{d}^2N_{\text {data}}}{\mathrm{d} \Delta \eta \mathrm{d} \Delta \phi }, \quad M(\Delta \eta ,\Delta \phi )=\frac{\mathrm{d}^2N_{\text {mixed}}}{\mathrm{d} \Delta \eta \mathrm{d} \Delta \phi } \end{aligned}$$

are the distributions of particle pairs from data and mixed events (representing uncorrelated background), respectively. The mixed events were constructed using the mixing procedure described in Sect. 5.1. Distributions \(D(\Delta \eta ,\Delta \phi )\) and \(M(\Delta \eta ,\Delta \phi )\) were obtained by accumulating the number of pairs in intervals of \(\Delta \eta \) and \(\Delta \phi \). For the calculation of \(C(\Delta \eta ,\Delta \phi )\) both distributions were normalised to the number of pairs (\(N^{\text {pairs}}_{\text {data}}\), \(N^{\text {pairs}}_{\text {mixed}}\)) in the given distribution. The measurements are restricted to \(0 \le \Delta \eta \le 3\) and \(0 \le \Delta \phi < \pi \). Finally the correlation functions were mirrored around the point \((\Delta \eta ,\Delta \phi ) = (0,0)\) and plotted in the range \( -\frac{\pi }{2}< \Delta \phi < \frac{3 \pi }{2}\).

In this paper the correlation function \(C(\Delta \eta ,\Delta \phi )\) was calculated for charged primary hadrons produced in strong and electromagnetic processes with transverse momentum \(p_T < 1.5\) GeV/c in the NA61/SHINE acceptance (for details see [14]).

While CMS [15], ATLAS [16] and ALICE [17] used equivalent definitions, other quantities have also been used to measure two-particle correlations. The ACM collaboration [4] at the ISR studied

$$\begin{aligned} C^\mathrm{II}(\Delta \eta ,\Delta \phi )= & {} \langle (n-1)(\rho _n^\mathrm{II}(\eta _1,\phi _1,\eta _2,\phi _2) \\&- \rho _n^{\text {I}}(\eta _1,\phi _1)\rho _n^{I}(\eta _2,\phi _2) \rangle , \end{aligned}$$

where n is the multiplicity in the event, \(\rho _n^\mathrm{II}\) the pair density, \(\rho _n^{\text {I}}\) the inclusive density and the averaging runs over the events.

PHOBOS [5] employed the observable

$$\begin{aligned} R(\Delta \eta ,\Delta \phi )=\langle (n-1)(\rho _n^\mathrm{II}(\Delta \eta ,\Delta \phi ) / \rho _n^{\text {mixed}}(\Delta \eta ,\Delta \phi ) - 1) \rangle \end{aligned}$$

with \(\rho _n^\mathrm{II}\) and \(\rho _n^{\text {mixed}}\) the pair density in data and mixed events, respectively, both normalised to unity.

The STAR collaboration [6] used the quantity

$$\begin{aligned}&\Delta \rho / \sqrt{\rho _\mathrm{mixed}}\\&\quad = (\rho (\Delta \eta ,\Delta \phi ) - \rho _\mathrm{mixed}(\Delta \eta ,\Delta \phi ) ) / \sqrt{\rho _\mathrm{mixed}(\Delta \eta ,\Delta \phi )}, \end{aligned}$$

where \(\rho \) and \(\rho _\mathrm{mixed}\) are the normalised densities of data and mixed pairs, respectively.

3 Experimental setup, beams, target, triggers

The NA61/SHINE experiment uses a large acceptance hadron spectrometer located in the CERN North Area [2]. The schematic layout of the NA61/SHINE detector is shown in Fig. 1.

Fig. 1
figure 1

The schematic layout of the NA61/SHINE experiment at the CERN SPS (horizontal cut, not to scale)

A set of scintillation and Cherenkov counters as well as beam position detectors (BPDs) upstream of the spectrometer provide timing reference, identification and position measurements of incoming beam particles. The trigger scintillator counter S4 placed downstream of the target is used to select events with collisions in the target area by the absence of a charged particle hit.

Secondary beams of positively charged hadrons at 20, 31, 40, 80 and 158 GeV/c are produced from 400 GeV/c protons extracted from the SPS accelerator. Particles of the secondary hadron beam are identified by two Cherenkov counters, a CEDAR [18] (either CEDAR-W or CEDAR-N) and a threshold counter (THC). The CEDAR counter, using a coincidence of six out of the eight photo-multipliers placed radially along the Cherenkov ring, provides positive identification of protons, while the THC, operated at pressure lower than the proton threshold, is used in anti-coincidence in the trigger logic. Due to their limited range of operation two different CEDAR counters were used, namely for beams at 20, 31, and 40 GeV/c the CEDAR-W counter and for beams at 80 and 158 GeV/c the CEDAR-N counter. The threshold counter was used for the lower three beam momenta. A selection based on signals from the Cherenkov counters allowed one to identify beam protons with a purity of about 99%. A consistent value for the purity was found by bending the beam into the TPCs with the full magnetic field and using identification based on its specific ionisation energy loss \(\mathrm{d}E/\mathrm{d}x\)  [19].

The main tracking devices of the spectrometer are four large volume Time Projection Chambers (TPCs). Two of them, the vertex TPCs (VTPC-1 and VTPC-2), are located in the magnetic fields of two super-conducting dipole magnets with a maximum combined bending power of 9 Tm, which corresponds to about 1.5 and 1.1 T fields in the upstream and downstream magnets, respectively. In order to optimise the acceptance of the detector, the field in both magnets was adjusted proportionally to the beam momentum.

Two large main TPCs (MTPC-L and MTPC-R) are positioned downstream of the magnets symmetrically to the beam line. The fifth small TPC (GAP-TPC) is placed between VTPC-1 and VTPC-2 directly on the beam line. It closes the gap between the beam axis and the sensitive volumes of the other TPCs. The TPCs are filled with Ar and CO\(_2\) gas mixtures. Particle identification in the TPCs is based on measurements of the specific energy loss (\(\mathrm{d}E/\mathrm{d}x\)) in the chamber gas and is augmented by mass measurements using time-of-flight (ToF) detectors.

The p + p data, which is the topic of this paper, was taken by colliding the proton beam with a liquid hydrogen target (LHT), a 20 cm long cylinder positioned about 80 cm upstream of VTPC-1.

4 Data processing, simulation and detector performance

Detector parameters were optimised by a data-based calibration procedure which also took into account their time dependence. Small adjustments were determined in consecutive steps for:

  1. (i)

    detector geometry, TPC drift velocities and distortions due to the magnetic field inhomogeneities in the corners of the VTPCs,

  2. (ii)

    magnetic field setting,

  3. (iii)

    specific energy loss measurements.

Each step involved reconstruction of the data required to optimise a given set of calibration constants and time dependent corrections followed by verification procedures. Details of the procedure and quality assessment are presented in Ref. [20]. The resulting performance in the measurements of quantities relevant for this paper is discussed below.

The main steps of the data reconstruction procedure are:

  1. (i)

    cluster finding in the TPC raw data, calculation of the cluster centre-of-gravity and total charge,

  2. (ii)

    reconstruction of local track segments in each TPC separately,

  3. (iii)

    matching of track segments into global tracks,

  4. (iv)

    track fitting through the magnetic field and determination of track parameters at the first measured TPC cluster,

  5. (v)

    determination of the interaction vertex using the beam trajectory (x and y coordinates) fitted in the BPDs and the trajectories of tracks reconstructed in the TPCs (z coordinate),

  6. (vi)

    refitting the particle trajectory using the interaction vertex as an additional point and determining the particle momentum at the interaction vertex.

A simulation of the NA61/SHINE detector response is used to correct the reconstructed data. For this purpose inelastic p + p collisions generated with the Epos 1.99 [21] model were used to obtain the corrections for reconstruction inefficiency and trigger bias of the NA61/SHINE detector and obtain final corrected results.

The simulation consists of the following steps (see Ref. [22] for more details):

  1. (i)

    generation of inelastic p + p interactions using the Epos model,

  2. (ii)

    propagation of outgoing particles through the detector material using the GEANT 3.21 package [23] which takes into account the magnetic field as well as relevant physics processes, such as particle interactions and decays,

  3. (iii)

    simulation of the detector response using dedicated NA61/SHINE packages which introduce distortions corresponding to all corrections applied to the real data,

  4. (iv)

    simulation of the interaction trigger selection by checking whether a charged particle hits the S4 counter, see Sect. 3,

  5. (v)

    storage of the simulated events in a file which has the same format as the raw data,

  6. (vi)

    reconstruction of the simulated events with the same reconstruction chain as used for the real data and

  7. (vii)

    matching of reconstructed with simulated tracks based on the cluster positions.

It should be underlined that only inelastic p + p interactions in the hydrogen within the target cell were simulated and reconstructed. Thus the simulation-based corrections can be applied only for inelastic events. The contribution of elastic events is removed by the event selection cuts (see Sect. 5.1).

5 Data selection and analysis

This section presents the procedures used for data analysis [24]. These consist of the following steps: applying event and particle selections, obtaining uncorrected experimental results, evaluation of correction factors based on simulations, and finally calculation of the corrected results.

These steps are described in the subsections below.

5.1 Event and particle selection

5.1.1 Event selection criteria

The events selected for the analysis had to satisfy the conditions:

  1. (i)

    event is recorded with the interaction trigger,

  2. (ii)

    the beam particle trajectory is measured in at least one of BPD-1 or BPD-2 and in the BPD-3 detector,

  3. (iii)

    no off-time beam particle is detected within \(\pm 1\) \(\upmu \)s around the trigger particle,

  4. (iv)

    there is at least one track reconstructed in the TPCs and fitted to the interaction vertex,

  5. (v)

    the vertex z position (fitted using the beam and TPC tracks) is not farther away than 10 cm from the centre of the LHT,

  6. (vi)

    events with a positively charged track with absolute momentum (p) close to the beam momentum (\(p > p_{\text {beam}} - 3\) GeV/c) are rejected to exclude elastic interactions.

The data statistics for analysed reactions is shown in Table 1.

5.1.2 Track selection criteria and acceptance

The tracks selected for the analysis had to satisfy the conditions:

  1. (i)

    track fit converged,

  2. (ii)

    the total number of reconstructed points on the track should be at least 30,

  3. (iii)

    the sum of the number of reconstructed points in VTPC-1 and VTPC-2 should be at least 15 or the number of reconstructed points in the GAP-TPC should be at least five,

  4. (iv)

    the distance between the track extrapolated to the interaction plane and the interaction point (impact parameter) should be smaller than 4 cm in the horizontal (bending) plane and 2 cm in the vertical (drift) plane,

  5. (v)

    transverse momentum of particle track should be lower than 1.5 GeV/c,

  6. (vi)

    tracks with \(\mathrm{d}E/\mathrm{d}x\) and total momentum values characteristic for electrons are rejected (see Ref. [25] for the details of this cut).

Model simulations were performed in \(4\pi \) acceptance, thus the NA61/SHINE detector acceptance filter needs to be applied before comparisons with data. The detector acceptance was defined as a three-dimensional matrix (p, \(p_T\), \(\phi \)) of 1 and 0 depending on whether the bin was or was not populated by particles reconstructed and accepted in the events (see Ref. [14]).

Table 1 Event statistics before and after event cuts for p + p data. See text for detailed description of the cuts
Fig. 2
figure 2

Example distributions of number of pairs of all charged particles in (\(\Delta \eta \), \(\Delta \phi \)) from the data (left) and mixed (right) events in p + p interactions in 158 GeV/c

5.2 Raw correlation function

Uncorrected two-particle correlation functions were evaluated using Eq. 2 for events and tracks selected according to the criteria stated above. This was done by calculating differences between each accepted track in pseudorapidity and azimuthal angle with every other track in the same event. The differences were evaluated in the centre-of-mass (CMS) frame and accumulated in two-dimensional bins of (\(\Delta \eta \), \(\Delta \phi \)). To increase statistics, the \(\Delta \phi \) range was folded, i.e. for \(\Delta \phi \) larger than \(\pi \) its value was recalculated as \(2\pi - \Delta \phi \). This is allowed due to the symmetry of the p + p reactions. The transformation from laboratory (LAB) to CMS frame was performed assuming the pion mass for all produced particles.

The uncorrelated reference was constructed by mixing particles between events with two main constraints: (a) the multiplicity distribution of mixed data had to be exactly the same as the original; (b) mixed events could not contain two particles from the same original (data) event. Example distributions of number of pairs of charged particles in bins of (\(\Delta \eta \), \(\Delta \phi \)) before and after mixing are shown in Fig. 2. The overall shape of the pair distribution in Fig. 2 (left) results from the shape of the inclusive single particle distribution whereas the hump structure is caused by the acceptance limitations of the detector. As shown in Fig. 2 (right) these features are reproduced in the mixed pair distribution, which does not contain effects of true correlations by construction. Thus in the correlation function defined as the ratio in Eq. 2 these background effects cancel and only the true correlation effects remain.

Examples of raw correlation functions, uncorrected for trigger bias, track selection cuts and reconstruction efficiency are presented in the left column of Fig. 3 for all charge pairs.Footnote 2

Fig. 3
figure 3

Examples of uncorrected correlation functions from data (left plots), correction factors (centre plots) and corrected two-particle correlation functions \(C(\Delta \eta ,\Delta \phi )\) (right plots) for all charge pairs. Results are for p + p interactions at 20 GeV/c (top) and 158 GeV/c (bottom)

5.3 Corrections and uncertainties

5.3.1 Correction for reconstruction inefficiency and trigger bias

In order to correct the results for biases due to trigger and off-line event and track selection, detection efficiency, contribution of weak decays and secondary interaction products, the same analysis was also performed on simulated data. These effects, in general, may change the correlation function. In particular, weak decays of \(\Lambda \) and \(K^0_S\) hadrons lead to production of positively and negatively charged hadrons which are correlated by weak decay kinematics. The EPOS [21] model was used for event generation since it provides a good description of the NA61/SHINE results on the yields of both non-strange and strange particles in p + p collisions [25,26,27]. Correction factors \(\text {Corr}(\Delta \eta ,\Delta \phi )\) were calculated bin-by-bin as the ratio of the correlation functions for simulated events from the Epos  [21] model (“pure”) and the same events after processing through Geant [23] detector simulation and reconstruction (“rec”), and filtered using event and track selection cuts.

$$\begin{aligned} \text {Corr}(\Delta \eta ,\Delta \phi ) = \frac{C^{\text {sim}}_{\text {pure}}(\Delta \eta ,\Delta \phi )}{C^{\text {sim}}_{\text {rec}}(\Delta \eta ,\Delta \phi )}, \end{aligned}$$
(3)

where \(C^{\text {sim}}_{\text {pure}}\) and \(C^{\text {sim}}_{\text {rec}}\) are the correlation function \(C(\Delta \eta ,\Delta \phi )\) obtained for simulated events before and after detector simulation and reconstruction. For both “pure” and “rec” events the NA61/SHINE acceptance filter was applied.

Corrected results \(C(\Delta \eta ,\Delta \phi )\) were obtained by multiplying the uncorrected correlation function by the corresponding correction, namely:

$$\begin{aligned} C(\Delta \eta ,\Delta \phi ) = C^{\text {raw}}(\Delta \eta ,\Delta \phi ) \cdot \text {Corr}(\Delta \eta ,\Delta \phi ). \end{aligned}$$
(4)
Table 2 Event (top) and track (bottom) selection cuts. The standard cuts (centre) are used to obtain the final results, whereas the loose (left) and tight (right) cuts are employed to estimate systematic uncertainties (see Sects. 5.1.1 and 5.1.2, respectively)
Fig. 4
figure 4

Two-particle correlation function \(C(\Delta \eta ,\Delta \phi )\) for all charge pairs in inelastic p + p interactions at 20–158 GeV/c. The correlation function is mirrored around \((\Delta \eta ,\Delta \phi )=(0,0)\)

Fig. 5
figure 5

Two-particle correlation function \(C(\Delta \eta ,\Delta \phi )\) for unlike-sign pairs in inelastic p + p interactions at 20–158 GeV/c. The correlation function is mirrored around \((\Delta \eta ,\Delta \phi )=(0,0)\)

Fig. 6
figure 6

Two-particle correlation function \(C(\Delta \eta ,\Delta \phi )\) for positive charge pairs in inelastic p + p interactions at 20–158 GeV/c. The correlation function is mirrored around \((\Delta \eta ,\Delta \phi )=(0,0)\)

Fig. 7
figure 7

Two-particle correlation function \(C(\Delta \eta ,\Delta \phi )\) for negative charge pairs in inelastic p + p interactions at 20–158 GeV/c. The correlation function is mirrored around \((\Delta \eta ,\Delta \phi )=(0,0)\)

Correction factors \(\text {Corr}(\Delta \eta ,\Delta \phi )\) and corrected correlation functions \(C(\Delta \eta ,\Delta \phi )\) for all charge pairs in p + p interactions at 20 and 158 GeV/c are shown in Fig. 3 as examples. The values of the correction factors range between 0.9 and 1.1.

Fig. 8
figure 8

Two-particle correlation function \(C(\Delta \eta )\) obtained from projection of \(C(\Delta \eta ,\Delta \phi )\) onto the \(\Delta \eta \) axis for two subranges of \(\Delta \phi \). Left column shows the results for all charge pairs, middle column unlike-sign pairs, right column positive charge pairs. Vertical bars denote statistical and shaded regions denote systematic uncertainties. Predictions of the Epos model are shown by solid curves and the UrQMD model by dotted curves. Legend applies to all panels

Fig. 9
figure 9

Two-particle correlation function \(C(\Delta \phi )\) obtained from projection of \(C(\Delta \eta ,\Delta \phi )\) onto the \(\Delta \phi \) axis for two subranges of \(\Delta \eta \). Left column shows the results for all charge pairs, middle column unlike-sign pairs, right column positive charge pairs. Vertical bars denote statistical and shaded regions denote systematic uncertainties. Predictions of the Epos model are shown by solid curves and the UrQMD model by dotted curves. Legend applies to all panels

5.3.2 Statistical uncertainties

Statistical uncertainties of the correlation function are calculated as

$$\begin{aligned} \sigma ^2(C) = \sqrt{\left[ \text {Corr} \cdot \sigma (C^{\text {raw}}) \right] ^2 + \left[ C^{\text {raw}} \cdot \sigma (\text {Corr}) \right] ^2} \end{aligned}$$
(5)

in each \((\Delta \eta ,\Delta \phi )\) bin. Statistical uncertainties are found to be approximately independent of \(\Delta \phi \) but increase significantly with increasing \(\Delta \eta \). Depending on the beam momentum and charge combination the statistical uncertainties of the C function in individual (\(\Delta \eta ,\Delta \phi \)) bins are about 5% for \(\Delta \eta \approx 0\) (all charge pairs at top beam momentum) to more than 20% for \(\Delta \eta \approx 3\) (negatively charge pairs at the lowest beam momentum).

5.3.3 Estimation of systematic uncertainties

In order to estimate systematic uncertainties the data were analysed with loose and tight event and track selection cuts. By changing the cuts, one changes the magnitude of the corrections due to various biasing effects. In case of simulation perfectly reproducing the data, corrected results should be independent of the cuts. A dependence on cuts is due to imperfections of the simulation and is used as an estimate of the systematic uncertainty. For example, systematic uncertainty caused by weakly decaying particles is estimated by varying \(B_x\) and \(B_y\) cuts. The standard set of cuts was presented in Sect. 5.1 and is tabulated in Table 2 together with loose and tight cuts.

Results for both sets of cuts were subtracted bin-by-bin (loose–tight). Since the differences in all bins follow Gaussian distributions with mean close to 0, the systematic uncertainties were estimated as the standard deviation of the distribution. The procedure was performed for all charge combinations (all charge, unlike-sign, positively and negatively charge pairs) and at 20, 40 and 158 GeV/c. The systematic uncertainties are generally below 1%. The exception are correlation functions at 20 GeV/c at higher \(\Delta \eta \), where the systematic uncertainties are about 5%.

6 Results and discussion

6.1 Two-particle correlation function \(C(\Delta \eta ,\Delta \phi )\)

The corrected correlation functions for all charge pair combinations (all charge pairs, unlike-sign pairs, positively and negatively charge pairs) are presented in Figs. 4, 5, 6 and 7, respectively. Their values lie in the range between 0.8 and 1.4. Vanishing two-particle correlations result in \(C = 1\).

Fig. 10
figure 10

NA61/SHINE detector acceptance effects studied with the Epos  [21] model. Left plots show correlation functions in almost complete acceptance (only \(p_T < 1.5\) GeV/c cut was applied), middle plots show correlation functions in NA61/SHINE acceptance, right plots ratios of these two. Top row shows results for 20 GeV/c, bottom row for 158 GeV/c

Fig. 11
figure 11

The ratio of correlation functions \(C(\Delta \eta ,\Delta \phi )\) for data and models. The two left columns show results for unlike-sign pairs and the two right columns results for positive charge pairs. Each row shows results for one beam momentum

Statistical and systematic uncertainties were calculated using the procedures described in Sects. 5.3.2 and 5.3.3. The numerical values of the correlation functions and of statistical and systematic uncertainties are available in Ref. [28].

The main features of the results are:

  1. (i)

    A maximum at \((\Delta \eta ,\Delta \phi ) \approx (0,\pi )\), most prominent for all charge and unlike-sign pairs and significantly weaker for like-sign pairs. The most probable explanation is the contribution from resonance decays, mostly from abundantly produced \(\rho ^0 \rightarrow \pi ^+ + \pi ^-\). The weaker maximum for positive charge pairs can be attributed to e.g. \(\Delta ^{++}\) resonance decay. No such maximum is observed for negative charge pairs consistent with the fact that there are almost no resonances decaying into two or more negative charge pairs.

  2. (ii)

    A \(\cos (\Delta \phi )\) modulation appearing as a minimum near \(\Delta \phi = 0\) and maximum near \(\Delta \phi = \pi \) at all values of \(\Delta \eta \) for all combinations of charges. Stronger in all charge and unlike-sign pairs, weaker but still visible in positive charge and barely noticeable for negative charge pairs. The structure is probably due to momentum conservation.

  3. (iii)

    A Gaussian-like enhancement around \(\Delta \eta = 0\) along the full \(\Delta \phi \) range, clearly visible for all charge and unlike-sign pairs, however, significantly weaker, but still noticeable in like-sign pairs. This feature may be connected with string fragmentation or flux-tube fragmentation [29].

  4. (iv)

    A hill around \((\Delta \eta ,\Delta \phi ) = (0,0)\) for like-sign pairs. For positive charged pairs it grows with increasing beam momentum, but it is independent of beam momentum for negative charge pairs. Since the products of \(\gamma \) decay were rejected during the analysis, the hill is probably caused by Bose–Einstein statistics.

  5. (v)

    Clearly, there is no jet-like peak at \((\Delta \eta ,\Delta \phi ) = (0,0)\). This indicates that the contribution of hard scattering effects is small as expected from their small probability at SPS energies (omitting the cut \(p_T < 1.5\) GeV/c was tried and did not change the results).

6.2 Projections onto \(C(\Delta \eta )\) and \(C(\Delta \phi )\)

To make the results easier to compare quantitatively, slices of the two-dimensional correlation function are also presented. The \(C(\Delta \eta )\) function was obtained by projecting \(C(\Delta \eta ,\Delta \phi )\) onto \(\Delta \eta \) in the \(\Delta \phi \) intervals: \(0< \Delta \phi < \pi /4\) and \(3\pi /4< \Delta \phi < \pi \). The results for 20, 40, and 158 GeV/c beam momenta and all charge, unlike-sign and positive charged pairs, together with uncertainties as well as predictions of the Epos and UrQMD models are presented in Fig. 8. Note that due to the projection the statistical uncertainties of the projected correlation function are considerably smaller than for the individual two-dimensional bins.

The two-particle correlation function \(C(\Delta \phi )\) was obtained from projection of \(C(\Delta \eta ,\Delta \phi )\) onto the \(\Delta \phi \) axis in two \(\Delta \eta \) intervals: \(0< \Delta \eta < 1\) and \(2< \Delta \eta < 3\). The results are presented in Fig. 9. The tendencies shown in Figs. 8 and 9 are reproduced well by the Epos model.

6.3 Acceptance effects

Acceptance effects were studied using events generated with the Epos model. Namely, correlation functions \(C(\Delta \eta ,\Delta \phi )\) were produced for Epos events in almost complete acceptance (full range of pseudorapidity and azimuthal angle but \(p_T < 1.5\) GeV/c) and the NA61/SHINE detector acceptance (applying the particle population matrix [14]). The results (Fig. 10) show that the correlation functions are qualitatively similar (only the enhancement close to \(\Delta \phi \approx \pi \) is 3–7% stronger with the near complete acceptance). The right panels of Fig. 10 display the ratios of the two correlation functions. Thus, the Epos model suggests that the effects of the limitation of the NA61/SHINE acceptance are not large. However, it is important to point out that correlations cannot be corrected for incomplete phase space acceptance in a model independent way. Therefore it is essential to apply the NA61/SHINE acceptance filter to model calculations before comparing to the measurements.

7 Comparison with models

7.1 Comparison with the Epos and UrQMD models

In order to compare the experimental results and predictions of the Epos  1.99 [21] and UrQMD 3.4 [30, 31] models the ratio of correlation functions obtained from data and models was calculated. The same acceptance was used for the NA61/SHINE data and for the Epos and UrQMD models (see Ref. [14]). The results for unlike-sign and positive charge pairs are shown in Fig. 11. In general, Epos  1.99 reproduces the experimental results from p + p interactions well (see first and third column in Fig. 11). However, the model does not reproduce the peak of about 20% at \((\Delta \eta ,\Delta \phi )=(0,0)\) seen for positive charge pairs (third column in Fig. 11) at 40, 80 and 158 GeV/c. This demonstrates an excess of the number of real data over model pairs in that region. Note that the Epos model does not include either Bose–Einstein correlations or Coulomb interactions which are expected to produce this kind of correlation. Moreover, a weak enhancement of the ratio of about 10% is seen near \(\Delta \eta \approx 0\) and \(\Delta \phi \approx \pi \text { and } \frac{3}{4}\pi \) (first column in Fig. 11). This suggests a slight underestimate of the correlation in the valleys of the \(\cos (\Delta \phi )\) modulation in Epos.

The ratio of the experimental results and the predictions of the UrQMD 3.4 model is shown in the second and fourth column of Fig. 11. The UrQMD model reproduces data less well than Epos. The most visible discrepancy is a \(\Delta \eta \)-wide suppression of the ratio in unlike-sign pairs (second column in Fig. 11) which suggests that UrQMD generates 25–30% stronger long-range near-side (\(\Delta \phi \approx 0\)) correlations than existing in the data. Another inconsistency between data and model results can be seen in positive charge correlations (fourth column in Fig. 11) where a suppression of up to 25% is present on the away-side (\(\Delta \phi \approx \pi \)), strongest at the three lowest beam momenta, as well as an enhancement of more than 10% at large \(\Delta \eta \) on the same side in azimuthal angle.

8 Summary

Two-particle correlations of charged particles in azimuthal angle and pseudorapidity were measured in inelastic p + p collisions at 20–158 GeV/c by the NA61/SHINE experiment at the CERN SPS. The results show structures which can be connected to phenomena such as resonance decays, momentum conservation and Bose–Einstein correlations. A comparison with the Epos and UrQMD models was performed. The predictions of the Epos model are close to the NA61/SHINE data whereas the results from UrQMD show larger discrepancies.