Measurements of multiplicity fluctuations of identified hadrons in inelastic proton-proton interactions at the CERN Super Proton Synchrotron

Measurements of multiplicity fluctuations of identified hadrons produced in inelastic p+p interactions at 31, 40, 80, and 158~\GeVc beam momentum are presented. Three different measures of multiplicity fluctuations are used: the scaled variance $\omega$ and strongly intensive measures $\Sigma$ and $\Delta$. These fluctuation measures involve second and first moments of joint multiplicity distributions. Data analysis is performed using the Identity method which corrects for incomplete particle identification. Strongly intensive quantities are calculated in order to allow for a direct comparison to corresponding results on nucleus-nucleus collisions. The results for different hadron types are shown as a function of collision energy. A comparison with predictions of string-resonance Monte-Carlo models: Epos, Smash and Venus, is also presented.


Introduction
This paper presents experimental results on event-by-event fluctuations of multiplicities of identified particles produced in inelastic proton-proton (p+p) interactions at 31,40,80, and 158 GeV/c ( √ s NN = 7. 6,8.7,12.3,17.3 GeV). The measurements were performed by the multi-purpose NA61/SHINE [1] experiment at the CERN Super Proton Synchrotron (SPS) in 2009. 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 was performed [2].
An interpretation of the experimental results on nucleus-nucleus (A+A) collisions relies to a large extent on a comparison with the corresponding data on p+p and p+A interactions. In addition models of nucleus-nucleus collisions are often tuned based on results on p+p interactions. However, available results on fluctuations of identified hadrons in these reactions are sparse. Moreover, fluctuation measurements cannot be corrected in a model independent manner for partial phase-space acceptance. Thus all measurements of the scan should be performed in the same phase space region. This motivated the NA61/SHINE Collaboration to analyse data on p+p interactions with respect to fluctuations using the same experimental methods, acceptance and measures as used to study nucleus-nucleus collisions.
Fluctuations in A+A collisions are susceptible to two trivial sources: the finite and fluctuating number of produced particles and event-by-event fluctuations of the collision geometry. Suitable statistical tools have to be chosen to extract the fluctuations of interest. In this publication three different event-by-event fluctuation measures are used: the scaled variance ω, the ∆ and Σ measures introduced in Refs. [3,4]. All of them were already successfully utilized by the NA49 experiment at the CERN SPS, see e.g.
Experimental measurements of multiplicity distributions of identified hadrons are challenging because it is very difficult to identify a particle with sufficient precision. In this paper the Identity method [13,14,15,16,17,18,19] is employed to circumvent this problem. The Identity method has already been successfully used in the past by collaborations NA49 [8], NA61/SHINE [20], and ALICE [21,22,23].
The paper is organized as follows. In Sec. 2 intensive and strongly intensive measures of fluctuations used in this analysis are introduced and briefly discussed. The Identity method which allows to take into account the incomplete particle identification is presented in Sec. 3 and Appendix A. The NA61/SHINE set-up and the data reconstruction method are presented in Secs. 4 and 5, respectively. The data analysis procedure is introduced in Secs. 6 and 7. Applied corrections and remaining uncertainties are presented in Sec. 8. Results on the collision energy dependence of multiplicity fluctuations of identified hadrons in inelastic p+p collisions at 31, 40, 80, and 158 GeV/c beam momentum are presented, discussed and compared with model predictions in Sec. 9.
Throughout this paper the rapidity is calculated in the collision center of mass system: y = atanh(β L ), where β L = p L /E is the longitudinal (z) component of the velocity, p L and E are particle longitudinal momentum and energy given in the collision center of mass system. The transverse component of the momentum is denoted as p T and the azimuthal angle φ is the angle between the transverse momentum vector and the horizontal (x) axis. Total momentum in the laboratory system is denoted as p lab and electric charge is denoted as q. The collision energy per nucleon pair in the center of mass system is denoted as √ s NN .
2. Intensive and strongly intensive measures of multiplicity and particle type fluctuations

Intensive quantities
Measures of multiplicities and fluctuations are called intensive when they are independent of the volume (V) of systems modelled by the ideal Boltzmann grand canonical ensemble (IB-GCE). In contrast, extensive quantities (for example mean multiplicity or variance of the multiplicity distribution) are proportional to the system volume within IB-GCE. One can also extend the notion of intensive and extensive quantities to the Wounded Nucleon Model (WNM) [24], where the intensive quantities are those which are independent of the number of wounded nucleons (W), and extensive those which are proportional to the number of wounded nucleons. Here it is assumed that the number of wounded nucleons is the same for all collisions. The ratio of two extensive quantities is an intensive quantity [3]. Therefore, the ratio of mean multiplicities N a and N b , as well as the scaled variance of the multiplicity distribution ω[a] ≡ ( N 2 a − N a 2 )/ N a , are intensive measures. As a matter of fact, due to its intensity property, the scaled variance of the multiplicity distribution ω[a] is widely used to quantify multiplicity fluctuations in high-energy heavy-ion experiments. The scaled variance takes the value ω[a] = 0 for N a = const. and ω[a] = 1 for a Poisson distribution of N a .

Strongly intensive quantities
In nucleus-nucleus collisions the volume of the produced matter (or number of wounded nucleons) cannot be fixed -it changes from one event to another. The quantities, which within the IB-GCE (or WNM) model are independent of V (or W) fluctuations are called strongly intensive quantities [25,3]. The ratio of mean multiplicities is both an intensive and a strongly intensive quantity, whereas the scaled variance is an intensive but not strongly intensive quantity.
Strongly intensive quantities ∆ and Σ used in this paper are defined as [4]: and where N a and N b stand for multiplicities of particles of type a and b, respectively. First and second pure moments, N a , N b , and N 2 In addition, the second mixed moment, N a N b , is needed to calculate Σ[a, b].
With the normalization of ∆ and Σ used here [4], the quantities ∆[a, b] and Σ[a, b] are dimensionless and have a common scale required for a quantitative comparison of fluctuations of different, in general dimensional, extensive quantities. The values of ∆ and Σ are equal to zero in the absence of event-byevent fluctuations (N a = const., N b = const.) and equal to one for fluctuations given by the model of independent particle production (Independent Particle Model) [4]. The model assumes that particle types are attributed to particles independent of each other. Positive correlations between particle types, for example π + and π − coming in pairs from ρ 0 decays, lead to ∆ and Σ values below one. Anti-correlations between particle types, for example due to conservation laws, for example energy conservation leads to anti-correlation of multiplicities of different hadron types, may increase ∆ and Σ above one. For detailed discussion see Refs. [26,27].

Identity method
Experimental measurement of a joint multiplicity distribution of identified hadrons is challenging. Typical tracking detectors, like time projection chambers used by NA61/SHINE, allow for a precise measurement of momenta of charged particles and sign of their electric charges. In order to be able to distinguish between different particle types (e.g. a particle type a being e + , π + , K + or p) a determination of particle mass is necessary. This is done indirectly by measuring for each particle a value of the specific energy loss dE/dx in the tracking detectors, the distribution of which depends on mass, momentum and charge. The resolution of dE/dx measurements is not sufficient for particle-by-particle identification without a radical reduction of considered statistics. Probabilities to register particles of different types with the same value of dE/dx may be comparable. Consequently, it is impossible to identify particles individually with reasonable confidence for fluctuation analysis. The Identity method [13,14,15,16,17,18] is a tool to measure moments of multiplicity distribution of identified particles, which circumvents the experimental issue of incomplete particle identification.
The method employs the fitted inclusive dE/dx distribution functions of particles of type a, ρ a (dE/dx) in momentum bins. Each event has a set of measured dE/dx values corresponding to each track in the event. For each track in an event the probability w a of being a particle of type a is calculated: where Next, an event variable W a (a smeared multiplicity of particle a in the event) is defined as: where N is the number of measured particles in the event. The Identity method unfolds moments of the true multiplicity distributions from moments of the smeared multiplicity distribution P(W a ) using a response matrix calculated from the measured ρ a (dE/dx) distributions [14].

Experimental setup
The NA61/SHINE experimental facility [1] consists of a large acceptance hadron spectrometer located in the H2 beam line of the CERN North Area. The schematic layout of the NA61/SHINE detector is shown in Fig. 1.
The results presented in this paper were obtained using measurement from the Time Projection Chambers (TPC), the Beam Position Detectors (BPD) and the beam and trigger counters. These detector components as well as the proton beam and the liquid hydrogen target (LHT) are briefly described below. Further information can be found in Refs. [1,28,29].
Secondary beams of positively charged hadrons at 31, 40, 80, and 158 GeV/c were produced from 400 GeV/c protons extracted from the SPS onto a beryllium target. A selection based on signals from a set of detectors along the H2 beam-line (Cerenkov detectors CEDAR, scintillation counters S, THC and BPDs (see inset in Fig. 1)) allowed to identify beam protons with a purity of about 99%. A coincidence of these signals provided the beam trigger T beam . For data taking on p+p interactions a liquid hydrogen target of 20.29 cm length (2.8% interaction length) and 3 cm diameter was placed 88.4 cm upstream of the first TPC (see Fig. 1). The interaction trigger T int was provided by the anti-coincidence of the incoming proton beam and a scintillation counter S4 (T int = T beam ∧ S4). The S4 counter with 2 cm diameter, was placed between the magnets along the beam trajectory at about 3.7 m from the target, see Fig. 1.
The main tracking devices of the spectrometer are four large volume TPCs. 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 T and 1.1 T fields in the upstream and downstream magnets, respectively. In order to optimize the acceptance of the detector, the fields in both magnets were 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. Simultaneous measurements of dE/dx and p lab allow to extract information on particle mass, which is used to identify charged particles.
Behind the MTPCs there are three Time-of-Flight (ToF) detectors.

Data reconstruction and simulation
The event vertex and the produced particle tracks were reconstructed using the standard NA61/SHINE software [28]. Detector parameters were optimized by a data-based calibration procedure which also took into account their time dependence, for details see Refs. [28,30].
A simulation of the NA61/SHINE detector response was used to correct the reconstructed data. Several Monte Carlo models were compared with the NA61/SHINE results on p+p, p+C and π+C interactions [28,29,31,32]. Based on these comparisons and taking into account continuous support the Epos1.99 model [33,34] was selected for Monte Carlo simulations and calculation of corrections. In order to estimate systematic uncertainties simulations were also performed using Venus4.12, which was previously used by the NA49 Collaboration at the CERN SPS energies [35,9]. Generated and reconstructed tracks were matched based on the number of common points along their path. Possible differences due to the different identification procedures followed in the MC simulations and the real data are addressed in Ref. [29] and Sec. 8.3.
Since the contribution of elastic events is removed by the event selection (see Sec. 6), only inelastic p+p interactions in the hydrogen of the target cell were simulated and reconstructed. Thus the MC based corrections (see Sec. 8.1) can be applied only for inelastic events.

Event and track selection
The final results presented in this paper refer to identified hadrons produced in inelastic p+p interactions by strong interaction processes and in electromagnetic decays of produced hadrons. Such hadrons are referred to as primary hadrons. The event and track selection cuts described below are selected in order to minimize unavoidable biases in measured data with respect to the final results.

Event selection
An event was recorded if there was no beam particle detected downstream of the target (S4 counter in Fig. 1). Reconstructed events with off-time beam proton detected within a time window of ±1.5 µs around the trigger proton were removed. The trajectory of the trigger proton was required to be measured in at least three planes out of four of BPD-1 and BPD-2 and in BPD-3. Events without the fitted primary interaction vertex were removed. The fitted interaction vertex z-coordinated was requested to be within ±20 cm from the LHT target center. Events with a single positively charged track with absolute momentum close to the beam momentum (for details see Ref. [28]) are removed in order to eliminate remaining elastic scattering reactions. This requirement removes less then 2% of all events in data at 31-80 GeV/c beam momenta. The same cut applied to the Epos simulated and reconstructed data removes less then 5% of all inelastic events. The cut was not applied to the 158 GeV/c data-set. The fraction of inelastic events removed by the event selection cuts estimated based on the Epos simulation ranges between 15% at 31 GeV/c and 22% at 158 GeV/c.

Track selection
In order to select good quality tracks of primary charged hadrons and to reduce the contamination of tracks from secondary interactions, weak decays and off-time interactions, the following track selection criteria were applied: (i) Track momentum fit at the interaction vertex should have converged.
(ii) Total number of reconstructed points used to fit the track trajectory should be greater than 30.
(iii) Sum of the number of reconstructed points in VTPC-1 and VTPC-2 should be greater than 15 or the number of reconstructed points in the GAP TPC should be greater than 4. . This defines the analysis acceptance given in Ref.
[36] in a form of three dimensional maps in particle momentum space. Examples of the analysis acceptance are shown in Fig. 4.
The event and track statistics after applying the selection criteria are summarized in Table 1.

Identity analysis
In order to calculate moments of multiplicity distributions of identified hadrons corrected for incomplete particle identification the analysis was performed using the Identity method. The analysis consists of three steps: (i) parametrization of inclusive dE/dx spectra, (ii) calculation of smeared multiplicity distributions and their moments, (iii) correction of smeared moments for incomplete particle identification using the dE/dx response matrix.
The Identity analysis steps are briefly described below and in App. A.
For each particle its specific energy loss dE/dx is calculated as the truncated mean (smallest 50%) of cluster charges measured along the track trajectory. As an example, dE/dx measured in p+p interactions at 80 GeV/c, for positively and negatively charged particles, as a function of q · p lab is presented in Fig. 2 Figure 2: Example distribution of charged particles in the dE/dxq · p lab plane in p+p interactions at 80 GeV/c. Expectations for the dependence of the mean dE/dx on p lab for the considered particle types are shown by the curves calculated based on the Bethe-Bloch function.
The parametrization of dE/dx spectra of e + , e − , π + , π − , K + , K − , p, andp were obtained by fitting the dE/dx distributions separately for positively and negatively charged particles in bins of p lab and transverse momentum p T .
The fitted function was defined as a sum of the dE/dx distributions of e + , π + , K + , and p for positively charged particles. The sum of the contributions of the corresponding antiparticles was used for negatively charged particles.
The details of this fitting procedure can be found in Refs. [29,37,38]. In contrast to the spectra analysis [29] separate fits were performed in order to extend acceptance by adding particles with negative p x /q, where p x is x-component of the particle total momentum in the laboratory system. Systematic uncertainties arising from the fitting procedure are estimated in Sec. 8. In order to ensure similar particle numbers in each bin, 20 logarithmic bins were chosen in p lab in the range 1 − 100 GeV/c. Furthermore, the data were binned in 20 equal p T intervals in the range 0 − 2 GeV/c.
The dE/dx spectrum for a given particle type was parametrized by the sum of asymmetric Gaussians with widths σ a,l depending on the particle type a and the number of points l measured in the TPCs. The peak position of the dE/dx distribution for particle type a is denoted as x a . The contribution of a reconstructed particle track to the fit function reads: where x is the dE/dx of the particle, n l is the number of tracks with number of points l and Y a is the amplitude of the contribution of particles of type a. The second sum is the weighted average of the lineshapes from the different numbers of measured points (proportional to track-length) in the sample. The quantity σ l is written as: where the width parameter σ 0 is assumed to be common for all particle types and bins. A 1/ √ l dependence on number of points is assumed following Ref. [39]. The asymmetry parameter δ is introduced in order to take into account a possible small asymmetry of the truncated mean distribution resulting from a strong asymmetry of the Landau energy loss distribution. Examples of fits for p+p interactions at 31 and 158 GeV/c are shown in Fig. 3.
In order to ensure good fit quality, only bins with number of tracks greater than 500 were used for further analysis. The Bethe-Bloch curves for different particle types cross each other at low values of the total momentum. Thus, the proposed technique is not sufficient for particle identification at low p lab and bins with p lab < 4.3 GeV/c were excluded from this analysis based solely on dE/dx. The requirement of at least 500 tracks with good quality dE/dx measurement in each p lab , p T bin reduces the acceptance available for the analysis. Due to different multiplicities the acceptance is different for positively and negatively charged particles. Moreover, it also changes with beam momentum. Thus, the largest acceptance was found for positively charged hadrons at 158 GeV/c and the smallest at 31 GeV/c for negatively charged hadrons. The acceptance used in this analysis is given separately for negatively and positively charged particles by a set of publicly available acceptance tables [36]. The corresponding rapidity and transverse momentum acceptances at 31 and 158 GeV/c are shown in Fig. 4.
The parametrization of inclusive dE/dx spectra of identified particles is first used to calculate probabilities w a and, then, W a . The first moments of the multiplicity distributions for complete particle identification, N a are equal to the corresponding first moments of the smeared distributions:   The Identity method was quantitatively tested by numerous simulations, see for example Refs. [15,17].

Corrections and uncertainties
This section briefly describes the corrections for biases and presents methods to calculate statistical and systematic uncertainties.

Corrections for event and track losses and contribution of unwanted tracks
The first and second moments of multiplicity distributions corrected for incomplete particle identification were also corrected for: (i) loss of inelastic events due to the on-line and off-line event selection, (ii) loss of particles due to the detector inefficiency and track selection, (iii) contribution of particles from weak decays and secondary interactions (feed-down).
A simulation of the NA61/SHINE detector response was used to correct the data for the above mentioned biases. Corrections were calculated for moments of identified hadron multiplicity distributions. Events simulated with the Epos model were reconstructed with the standard NA61/SHINE software as described in Sec. 5. The multiplicative correction factors C (k) a and C ab , where a and b denote the particle type (a, b = π +/− , K +/− , p,p, e +/− ; and a b), are defined as: where: (i) (N k a ) MC gen -moment k of particle type a generated by the model, (ii) (N k a ) MC sel -moment k of particle type a generated by the model with the detector response simulation, reconstruction and selection, (iii) (N ab ) MC gen/sel -mixed second moment of particle types a and b generated by the model (gen) and with the detector response simulation, reconstruction and selection (sel).
This way of implementing the corrections was tested using simulations based on the Epos and Venus models, for details see Sec. 8.3.2. Multi-dimensional unfolding of distributions which would be the best approach for correcting experimental biases is too complex to be implemented [40]. This is why the Identity method -the unfolding of moments -was used to correct for the main bias -the incomplete particle identification. Then the unfolded moments were corrected for remaining biases.
The correction factors for first, second and mixed moments of identified hadrons are shown in Figs. 5, 6 and 7. Note that a single particle reconstruction inefficiency is typically lower than 10%, whereas the feed-down contribution is about 10% [28,29].

Statistical uncertainties
The sub-sample method was used to calculate statistical uncertainties of final results. All selected events were grouped into M = 30 non-overlapping sub-samples of events. Then a given fluctuation measure Q (for example Σ[π + , p]) was calculated for each sub-sample separately, and the variance of its distribution, Var[Q], was obtained. The statistical uncertainty of Q for all selected events was estimated as √ Var[Q]/M. The dE/dx parametrization requires a minimum number of tracks in a given momentum bin, thus the acceptance in which the dE/dx parametrization can be obtained is larger for all selected events than for sub-samples of events. In order to have the maximum acceptance the same dE/dx parametrization obtained using all events was used in the sub-sample analysis. It was checked using the bootstrap method [41,42] that the above approximation leads only to a small underestimation of statistical uncertainties.

Systematic uncertainties
Systematic uncertainties originate from imperfectness of the detector response and systematic uncertainties in the modelling of physics processes implemented in the models. The total systematic uncertainties were calculated by adding detector-related (see Sec. 8.3.1) and model-related (see Sec. 8.3.2) contributions in quadrature.

Detector related effects
These uncertainties were studied by applying standard (see Sec. 6) and loose cuts defined by: (i) reducing the rejection window for events with off-time beam particle to < 0.5 µs, For each choice the complete analysis was repeated including the dE/dx fitting and recalculating the corrections. Observed differences between results for the standard and loose cuts are related to imperfectness of the reconstruction procedure and to the acceptance of events with additional tracks from off-time particles. The corresponding systematic uncertainty was calculated as the difference between the results for the standard and loose cuts.
An additional possible source of uncertainty is imperfectness of the dE/dx parametrization. Here the largest uncertainty comes from uncertainties of the parameters of the kaon dE/dx distribution. The kaon distribution significantly overlaps with the proton and pion distributions. In the most difficult low momentum range the dE/dx fits were cross-checked using the time-of-flight information and found to be in agreement at the level of single particle spectra (see Ref. [29]).
In this analysis, as it considers second order moments, two additional tests were performed. First, fits of dE/dx distributions with fixed asymmetry parameter and without any constraint on asymmetry were used to estimate the resulting possible biases of fluctuation measures.  Table 2: Comparison of mean multiplicity in the analysis acceptance to mean multiplicity of identified hadrons in the full phase-space (only statistical uncertainty indicated) [29].
all, it is within or comparable to the systematic uncertainty. The only exceptions are the scaled variance of protons at 158 GeV/c (10% which normally is 5%) and pions at 31 GeV/c (15% which normally is 8%) as well as ∆ of pions and protons at 158 GeV/c (17% compared to 11%).

Model-related effects
Systematic uncertainty originating from imperfectness of the Epos model used to calculate corrections in describing p+p interactions is discussed here. The uncertainty was estimated using simulations performed within the Epos and Venus models. The simulated Epos data were corrected using corrections obtained based on the Venus model and compared to the unbiased Epos results. Then the same procedure was repeated swapping Epos and Venus. In both cases the correction improves agreement between the obtained results and the true ones. The differences between the unbiased and simulated-corrected results were added to the systematic uncertainty. They are on average about 20% (Epos data) and 25% (Venus data) of the total systematic uncertainty. Note that the models show similar agreement with results on p+p interactions at the CERN SPS energies.

Results, discussion and comparison with models
In this section final experimental results are presented and discussed as well as compared with predictions of selected string-hadronic models.

Results
The final results presented in this section refer to identified hadrons produced in inelastic p+p interactions by strong interaction processes and in electromagnetic decays of produced hadrons. They were obtained within the kinematic acceptances given in Ref.
[36] and illustrated in Fig. 4. Note, that the kinematic acceptances for positively and negatively charged hadrons are different. Mean multiplicities of pions, kaons and anti-protons in the acceptance region of the fluctuation analysis are plotted in Fig. 8 and compared to corresponding mean multiplicities measured in the full phase-space in Table 2.
Pions are the most abundantly produced particles and are the majority of accepted charged hadrons in all analyzed reactions. With decreasing beam momentum the contribution of protons increases and the small contributions of kaons and protons decrease. Almost all charged hadrons are pions, except that protons are the majority of positively charged hadrons at the lowest beam momentum, 31 GeV/c. The changes of particle type composition with charge of selected hadrons and beam momentum are related to different thresholds for production of pions, kaons and anti-protons. The mean proton multiplicity in the models in full phase-space is about one (0.3-0.4 in the acceptance) and approximately independent of beam momentum. This is because final state protons are strongly correlated with two initial state protons via baryon number conservation. Figure 9 shows the collision energy dependence of the scaled variance of pion, kaon and proton multiplicity distributions. Note, the intensive fluctuation measure ω is one for a Poisson distribution and zero in the case of a constant multiplicity for all collisions. The scaled variance quantifies the width of the multiplicity distribution relatively to the width of the Poisson distribution with the same mean multiplicity. The results for all charged, positively charged and negatively charged hadrons are presented separately. One observes: (i) the scaled variance for pions increases with the collision energy. The increase is the strongest for all charged pions. This is probably related to the well established KNO scaling of the charged hadron multiplicity distributions in inelastic p+p interactions with the scaled variance being proportional 10 15 [     Figure 11: The collision energy dependence of ∆ for different particle type combinations in inelastic p+p interactions. Results for all charged, positively and negatively charged hadrons are presented separately. The solid, dashed and dotted lines show predictions of Epos, Smash and Venus models, respectively. Statistical uncertainties are smaller than the symbol size and systematic uncertainty is indicated by red bands.
to mean multiplicity [43,44,45]. Global and local (resonance decays) electric charge conservation correlates multiplicities of positively and negatively charged pions and thus the effect is the most pronounced for all charged hadrons.
(ii) The dependence of ω on beam momentum for kaons is qualitatively similar to the one for pions but weaker. This is probably related to the significantly smaller mean multiplicity of kaons than pions. One notes that the scaled variance of the multiplicity distribution approaches one when the mean multiplicity decreases to zero. The latter effect is likely responsible for ω[K − ] and ω[p] being close to one, as the mean multiplicity of K − and p in the acceptance is below 0.1 and 0.03, respectively.
(iii) The scaled variance of protons is about 0.8 and depends weakly on the beam momentum. The netbaryon (baryon -anti-baryon) multiplicity in full phase-space is exactly two. This is because the initial baryon number is two and baryon number is conserved. Thus the scaled variance of the netbaryon multiplicity distribution is zero. Anti-baryon production at the SPS energies is small and thus the net-baryon multiplicity is close to the baryon multiplicity. The baryons are predominately protons and neutrons. Thus the proton fluctuations are expected to be mostly due to the fluctuation of the proton to neutron ratio and fluctuations caused by the limited phase space acceptance of protons.
Figures 10 shows the results on Σ for pion-proton, pion-kaon and proton-kaon multiplicities measured separately for all charged, positively charged, and negatively charged hadrons produced in inelastic p+p collisions at 31-158 GeV/c beam momentum. The Σ measure assumes the value one in the Independent Particle Production Model which postulates that particle types are attributed to particles independently of each other. This implies that Σ unlike ω is insensitive to the details of the particle multiplicity distribution. One observes: (i) For all and positively charged pions-proton combinations Σ is significantly below one (approximately 0.8) and depends weakly on the beam momentum. This is likely due to a large fraction of pion-proton pairs coming from decays of baryonic resonances [26,27]. Corresponding results for Pb+Pb collisions were reported in Refs. [46,47,48].
(ii) Σ for all charged pion-kaon combinations increases significantly with the beam momentum and is about 1.2 at 158 GeV/c. The origin of this behaviour is unclear.
(iii) For the remaining cases Σ is somewhat below or close to one suggesting a small contribution of hadrons from resonance decays. Figure 11 shows the results for ∆ of identified hadrons calculated separately for all charged, positively charged, and negatively charged hadrons produced in inelastic p+p collisions at beam momenta from 31 to 158 GeV/c. Results for ∆[π, p] at 31 and 40 GeV/c are not shown since they have large statistical uncertainties. This is because for these reactions N[π] ≈ N[p] (see Fig. 8) and thus C ∆ ≈ 0 (see Eq. 1.) The general properties of ∆ are similar to the properties of Σ discussed above. Unlike Σ, ∆ does not include a correlation term between multiplicities of two hadron types, see Eqs. 1 and 2. One observes: (i) ∆ for all and positively charged pions and protons is below one. This is qualitatively similar to Σ and thus likely to be caused by resonance decays.
(ii) ∆[(p +p), K] increases with the collision energy from about one to two. The origin of this dependence is unclear.
These models define the baseline for heavy ion collisions from which any critical phenomena are expected to emerge. However, the models should first be tuned to the experimental data on p+p interactions presented here. In p+p interactions at CERN SPS energies one expects none of the high matter density phenomena usually studied and searched for in nucleus-nucleus collisions. Any deviations from independent particle production are considered to be caused by well established effects discussed in Sec. 9.1. In general, the Epos and Venus models reproduce the results reasonably well. However, none of the models agrees with all features of the presented results. For a number a number of observables qualitative disagreement is observed for Smash and Venus.

Summary and outlook
In this paper experimental results on multiplicity fluctuations of identified hadrons produced in inelastic p+p interactions at 31, 40, 80, and 158 GeV/c beam momentum are presented. Results were corrected for incomplete particle identification using a data-based procedure -the Identity method. Remaining biases were corrected for utilising full physics and detector response simulations. The sub-sample method was used to calculate statistical uncertainties whereas systematic uncertainties were estimated by changing event and track selection cuts as well as models.
Results on the scaled variance of multiplicity fluctuations ω[a] of pions, kaons and protons for all charged, positively charged and negatively charged hadrons are presented. Moreover results on the strongly intensive measures of multiplicities fluctuations of two hadron types Σ[a, b] and ∆[a, b] are shown. These were obtained for pion-proton, pion-kaon and proton-kaon particle type combinations for all charged as well as separately for positively and negatively charged pair combinations. The results are presented as a function of the collision energy and discussed in the context of KNO-scaling, conservation laws and resonance decays.
Finally, the NA61/SHINE measurements are compared with string-resonance models Smash, Epos and Venus. In general, the Epos and Venus models reproduce the results reasonably well. However, none of the models agree with all presented results. For some observables even qualitative disagreement is observed for the Smash and Venus models. Thus, before the models can provide the baseline for heavy ion collisions in the search for critical phenomena, the models need to be tuned to the experimental data on p+p interactions presented in this paper.

A. Details on the Identity Method
The parametrization of inclusive dE/dx spectra of identified particles is first used to calculate probabilities w a (see Sec. 3). Distributions of w a for p+p interactions at 158 GeV/c are shown in Fig. 12 for positively and negatively charged particles, separately.
In the second step smeared multiplicities of identified particles W a (see Sec. 3) are calculated for each selected event and their distributions are obtained. Examples of smeared multiplicity distributions for p+p interactions at 158 GeV/c are shown in Fig. 13 for positively and negatively charged particles, separately.
Finally, first and second moments of smeared multiplicity distributions are calculated for positively and negatively charged particles, separately.