Investigation of the p–Σ 0 interaction via femtoscopy in pp collisions

This Letter presents the ﬁrst direct investigation of the p– (cid:2) 0 interaction, using the femtoscopy technique in high-multiplicity pp collisions at √ s = 13 TeV measured by the ALICE detector. The (cid:2) 0 is reconstructed via the decay channel to (cid:4) γ , and the subsequent decay of (cid:4) to p π − . The photon is detected via the conversion in material to e + e − pairs exploiting the capability of the ALICE detector to measure electrons at low transverse momenta. The measured p– (cid:2) 0 correlation indicates a shallow strong interaction. The comparison of the data to several theoretical predictions obtained employing the Correlation Analysis Tool using the Schrödinger Equation (CATS) and the Lednický–Lyuboshits approach shows that the current experimental precision does not yet allow to discriminate between different models, as it is the case for the available scattering and hypernuclei data. Nevertheless, the p– (cid:2) 0 correlation function is found to be sensitive to the strong interaction, and driven by the interplay of the different spin and isospin channels. This pioneering study demonstrates the feasibility of a femtoscopic measurement in the p– (cid:2) 0 channel and with the expected larger data samples in LHC Run 3 and Run 4, the p– (cid:2) 0 interaction will be constrained with high precision. the


Introduction
A quantitative understanding of the hyperon-nucleon interaction in the strangeness S = −1 sector is fundamental to pin down the role of strangeness within low energy quantum chromodynamics and to study the properties of baryonic matter at finite densities. The possible presence of the isoscalar and the isovector ( + , 0 , − ) hyperon states in the inner core of neutron stars (NS) is currently under debate due to the limited knowledge of the interaction of such hyperons with nuclear matter. The inclusion of hyperons in the description of the nuclear matter inside NS typically contains only states, and the on-average attractive nucleon-(N-) interaction leads to rather soft Equations of State (EoS) for NS. These are then unable to provide stability for stars of about two solar masses [1,2]. The hyperons are rarely included in the EoS for NS because of the limited knowledge about the N-strong interaction.
Indeed, while the attractive N-interaction is reasonably well constrained from the available scattering and light hypernuclei data [3][4][5], the nature of the N-interaction lacks conclusive experimental measurements. One of the major complications for experimental studies is the fact that the decay of all states involves neutral decay products [6], thus requiring high-resolution calorimeters.
The main source of experimental constraints on the N-system comes from scattering measurements [7][8][9], analysis of − atoms [10][11][12], and hypernuclei production data [13][14][15][16], although the latter are mainly dominated by large statistical uncertainties and large kaon decay background. Latest hypernuclear results obtained from different nuclear targets point towards an attractive interaction in the isospin I = 1/2 channel of the N-system [13,14], and repulsion in the I = 3/2 channel [15,16]. Hypernuclear measurements, however, are performed at nuclear saturation density and hence in the presence of more than one nucleon, resulting in a substantial model dependence in the interpretation of the experimental data [17].
Additionally, the hyperon-nucleon dynamics are strongly affected by the conversion process N-↔ N-, occurring in the I = 1/2 channel due to the close kinematic threshold between the two systems (about 80 MeV) [18][19][20][21][22]. This coupling is expected to provide an additional attractive contribution in the twobody N-interaction in vacuum [21,22]. Indeed, depending on the strength of the N-↔ N-coupling at the two-body level, the corresponding in-medium hyperon properties are very different. For a strong coupling, this leads to a repulsive single-particle potential U at large densities [21,22]. For the hyperon, the in-medium properties are mostly determined by the overall repulsion in the I = 3/2 component [21,22]. A repulsive component in the hyperon-nucleon interactions could shift the onset for hyperon production to larger densities, above 2-3 times the normal satura- tion density, thus leading to stiffer EoS which are able to describe the experimental constraint of NS.
To this end, different theoretical approaches including chiral effective field theories (χ EFT) [20] and meson-exchange models with hadronic [23] and quark [24] degrees of freedom provide a similar description of the available data by assuming a strong repulsion in the spin singlet S = 0, I = 1/2 and spin triplet S = 1, I = 3/2 and an overall attraction in the remaining channels. Recent ab initio lattice calculations at quark physical masses show a similar dependence on spin-isospin configurations for the central potential term [25]. The strength of the coupled-channel N-↔ Nis strongly model dependent as well. Calculations based on chiral models [20,21] and meson-exchange models [18,26] predict a rather strong or much weaker coupling, respectively. A selfconsistent description of this coupled-channel demands a detailed knowledge of the strong interaction in the N-system.
Recently, the study of two-particle correlations in momentum space measured in ultra-relativistic proton-proton (pp) and proton-nucleus collisions has proven to provide direct access to the interaction between particle pairs in vacuum [27][28][29]. The small size of the colliding systems of about 1 fm results in a pronounced correlation signal from strong final state interactions, which permits the latter to be precisely constrained. These measurements provided additional data in the hyperon sector with an unprecedented precision in the low momentum regime. In this Letter, these studies are extended to the sector. The electromagnetic decay of the 0 is exploited for the first direct measurement of the p-0 interaction in pp collisions. This study paves the way for extending these investigations to the charged states, in particular in light of the larger data samples expected from the LHC Runs 3 and 4.

Data analysis
This Letter presents results obtained from a data sample of pp collisions at √ s = 13 TeV recorded with the ALICE detector [30,31] during the LHC Run 2 (2015-2018). The sample was collected employing a high-multiplicity trigger with the V0 detectors, which consist of two small-angle plastic scintillator arrays located on either side of the collision vertex at pseudorapidities 2.8 < η < 5.1 and −3.7 < η < −1.7 [32]. The high-multiplicity trigger is defined by coincident hits in both V0 detectors synchronous with the LHC bunch crossing and by additionally requiring the sum of the measured signal amplitudes in the V0 to exceed a multiple of the average value in minimum bias collisions. This corresponds, at the analysis level, to the highest multiplicity interval containing the top 0.17% of all inelastic collisions with at least one charged particle in |η| < 1 (referred to as INEL > 0). This data set presents a suitable environment to study correlations due to the enhanced production of strange particles in such collisions [33]. Additionally, the larger charged-particle multiplicity density with respect to the minimum bias sample significantly increases the probability to detect particle pairs. The V0 is also employed to suppress background events, such as the interaction of beam particles with mechanical structures of the beam line, or beam-gas interactions.
In-bunch pile-up events with more than one collision per bunch crossing are rejected by evaluating the presence of additional event vertices [31]. The remaining contamination from pile-up events is on the percent level and does not influence the final results. Charged-particle tracking within the ALICE central barrel is conducted with the Inner Tracking System (ITS) [30] and the Time Projection Chamber (TPC) [34]. The detectors are immersed in a homogeneous 0.5 T solenoidal magnetic field along the beam direction. The ITS consists of six cylindrical layers of high positionresolution silicon detectors placed radially between 3.9 and 43 cm around the beam pipe. The two innermost layers are Silicon Pixel Detectors (SPD) and cover the pseudorapidity range |η| < 2.0 and |η| < 1.4, respectively. The two intermediate layers are composed of Silicon Drift Detectors, and the two outermost layers are made of double-sided Silicon micro-Strip Detectors (SSD), covering |η| < 0.9 and |η| < 1.0, respectively. The TPC consists of a 5 m long, cylindrical gaseous detector with full azimuthal coverage in the pseudorapidity range |η| < 0.9. Particle identification (PID) is conducted via the measurement of the specific ionization energy loss (dE/dx) with up to 159 reconstructed space points along the particle trajectory. The Time-Of-Flight (TOF) [35] detector system is located at a radial distance of 3.7 m from the nominal interaction point and consists of Multigap Resistive Plate Chambers covering the full azimuthal angle in |η| < 0.9. PID is accomplished by measuring the particle's velocity β via the time of flight of the particles in conjunction with their trajectory. The event collision time is provided as a combination of the measurements in the TOF and the T0 detector, two Cherenkov counter arrays placed at forward rapidity [36].
The primary event vertex (PV) is reconstructed with the combined track information of the ITS and the TPC, and independently with SPD tracklets. When both vertex reconstruction methods are available, the difference of the corresponding z coordinates is required to be smaller than 5 mm. A uniform detector coverage is assured by restricting the maximal deviation between the z coordinate of the reconstructed PV and the nominal interaction point to ±10 cm. A total of 1.0 × 10 9 high-multiplicity events are used for the analysis after event selection.
The proton candidates are reconstructed following the analysis methods used for minimum bias pp collisions at √ s = 7 TeV [27] and 13 TeV [28,29], and are selected from the charged-particle tracks reconstructed with the TPC in the kinematic range 0.5 < p T < 4.05 GeV/c and |η| < 0.8. The TPC and TOF PID capabilities are employed to select proton candidates by the deviation n σ between the signal hypothesis for a proton and the experimental measurement, normalized by the detector resolution σ . For candidates with p < 0.75 GeV/c, PID is performed with the TPC only, requiring |n σ | < 3. For larger momenta, the PID information of TPC and TOF are combined. Secondary particles stemming from weak decays or the interaction of primary particles with the detector material contaminate the signal. The corresponding fraction of primary and secondary protons are extracted using Monte Carlo (MC) template fits to the measured distribution of the Distance of Closest Approach (DCA) of the track to the primary vertex [27]. The MC templates are generated using PYTHIA 8.2 [37] and filtered through the ALICE detector [38] and reconstruction algorithm [30]. The resulting purity of protons is found to be 99%, with a primary fraction of 82%.
The 0 is reconstructed via the decay channel 0 → γ with a branching ratio of almost 100% [6]. The decay is characterized by a short life time rendering the decay products indistinguishable from primary particles produced in the initial collision. Due to the small mass difference between the and the 0 of about 77 MeV/c 2 , the γ has typically energies of only few hundreds of MeV. Therefore, it is reconstructed relying on conversions to e + e − pairs in the detector material of the central barrel exploiting the unique capability of the ALICE detector to identify electrons down to transverse momenta of 0.05 GeV/c. For transverse radii R < 180 cm and |η| < 0.9 the material budget corresponds to (11.4 ± 0.5)% of a radiation length X 0 , and accordingly to a conversion probability of (8.6 ± 0.4)% [39]. Details of the photon conversion analysis and the corresponding selection criteria are described in [39,40]. The reconstruction relies on the identification of secondary vertices by forming so-called V 0 decay candidates from two oppositely-charged tracks using a procedure described in detail in [41]. The products of the potential γ conversion are reconstructed with the TPC and the ITS in the kinematic range The Cosine of the Pointing Angle (CPA) between the γ momentum and the vector pointing from the PV to the decay vertex is required to be CPA > 0.999. In addition to the tight CPA selection, the contribution of particles stemming from out-of-bunch pile-up is suppressed by restricting the DCA of the photon to be along the beam direction (DCA z < 0.5 cm). After application of the selection criteria, about 946 × 10 6 γ candidates with a purity of about 95.4% are available for further processing.
The particle candidates are reconstructed via the subsequent decay → pπ − with a branching ratio of 63.9% [6], following the procedures described in [27,28]. For the the charge conjugate decay is exploited, and the same selection criteria are applied. The decay products are reconstructed with the TPC and the ITS within |η| < 0.9. The daughter candidates are identified by a broad PID selection in the TPC |n σ | < 5. The resulting candidate is obtained as the combination of the daughter tracks. The contribution of fake candidates is reduced by requesting a minimum p T > 0.3 GeV/c.
The coarse PID selection of the daughter tracks introduces a residual K 0 S contamination in the sample of the candidates. This contamination is removed by a 1.5σ rejection on the invariant mass assuming a decay into π + π − , where σ corresponds to the width of a Gaussian fitted to the K 0 S signal. Topological selections further enhance the purity of the sample. The radial distance of the decay vertex with respect to the detector centre ranges from 0.2 cm to 100 cm and CPA > 0.999. In addition to the tight CPA selection, particles stemming from out-of-bunch pile-up are rejected using the timing information of the SPD and SSD, and the TOF detector. One of the two daughter tracks is required to have a hit in one of these detectors. After application of the selection criteria, about 188 × 10 6 (178 × 10 6 ) ( ) candidates with a purity of 94.6% (95.3%) are available for further processing.
The 0 ( 0 ) candidates are obtained by combining all ( ) and γ candidates from the same event, where the nominal particle masses [6] are assumed for the daughters. In particular the timing selection on the daughter tracks of the assures that the 0 candidates stem from the right bunch crossing. In case a daughter track is used to construct two γ , , and candidates, or a combination thereof, the one with the smaller CPA is removed from the sample. In order to further optimize the yield and the purity of the sample, only 0 candidates with p T > 1 GeV/c are used.
The resulting invariant mass spectrum is shown in Fig. 1 for two p T intervals. In order to obtain the raw yield, the signal is fitted with a single Gaussian, and the background with a thirdorder polynomial. Due to the deteriorating momentum resolution for low p T tracks, the mean value of the Gaussian M 0 exhibits a slight p T dependence, which is well reproduced in MC simulations. The 0 ( 0 ) candidates for femtoscopy are selected as M 0 (p T ) ± 3 MeV/c 2 . The width of the interval is chosen as a compromise between the candidate counts and purity. In total, about 115 × 10 3 (110 × 10 3 ) 0 ( 0 ) candidates are found at a purity of about 34.6%. Due to the enhanced combinatorial background at low p T , the purity increases from about 20% at the lower p T threshold to its saturation value of about 60% above 5 GeV/c. Only one candidate per event is used, and is randomly selected in the very rare case in which more than one is available. In less than one per mille of the cases when the track of a primary proton is also employed as the daughter track of the γ or the , the corresponding 0 candidate is rejected. Since only strongly decaying resonances feed to the 0 [6], all candidates are considered to be primary particles.

Analysis of the correlation function
The experimental definition of the two-particle correlation function, for both p-p and p-0 pairs, is given by [44], with the same (N same ) and mixed (N mixed ) event distributions of k * and a normalization constant N . The relative momentum of the pair k * is defined as k * = 1 2 × |p * 1 − p * 2 |, where p * 1 and p * 2 are the momenta of the two particles in the pair rest frame, denoted by the * . The normalization is evaluated in k * ∈ [240, 340] MeV/c for p-p and in k * ∈ [250, 400] MeV/c for p-0 pairs, where effects of final state interactions are absent and hence the correlation function approaches unity.
The trajectories of the p-p and p-p pairs at low k * are almost collinear, and might therefore be affected by detector effects like track splitting and merging [45]. Accordingly, the reconstruction efficiency for pairs in the same and mixed event might differ. To this end, a close-pair rejection criterion is employed removing p-p and p-p pairs fulfilling η 2 + ϕ * 2 < 0.01, where the azimuthal coordinate ϕ * considers the track curvature in the magnetic field.
A total number of 1.7 × 10 6 (1.3 × 10 6 ) p-p (p-p ) and 587 (539) p-0 (p-0 ) pairs contribute to the respective correlation function in the region k * < 200 MeV/c. To enhance the statistical significance of the results, the correlation functions of baryonbaryon and antibaryon-antibaryon pairs are combined. Therefore, in the following p-0 denotes the combination p-0 ⊕ p-0 , and correspondingly for p-p.
The systematic uncertainties of the experimental correlation function are evaluated by simultaneously varying all proton, , γ , and 0 single-particle selection criteria by up to 20% around the nominal values. Only variations that modify the pair yield by less than 10% (20%) for p-0 (p-p) with respect to the default choice are considered, and the 0 purity by less than 5%. The impact of statistical fluctuations is reduced by evaluating the systematic uncertainties in intervals of 100 MeV/c (20 MeV/c) in k * for p-0 (p-p). The resulting systematic uncertainties are parametrized by an exponential function and interpolated to obtain the final pointby-point uncertainties. At the respectively lowest k * , the total systematic uncertainties are of the order of 2.5% for both p-p and p-0 .
Using the femtoscopy formalism [44], the correlation function can be related to the source function S(r * ) and the two-particle wave function ( r * , k * ) incorporating the interaction, where r * refers to the relative distance between the two particles. As demonstrated in [27][28][29]46] the correlation function becomes particularly sensitive to the strong interaction for small emission sources formed in pp and p-Pb collisions. For this study, a spherically symmetric emitting source is assumed, with a Gaussian shaped core density profile parametrized by a radius r 0 , which is obtained from a fit to the p-p correlation function, similarly as in [28,29]. Following the premise of a common emission source the such extracted radius is then used as an input to fit the p-0 correlation function. Possible modifications of the source profile due to the influence of strongly decaying resonances [47-49] are considered in the evaluation of the systematic uncertainties associated with the fitting procedure. The genuine p-p correlation function is modeled using the Correlation Analysis Tool using the Schrödinger equation (CATS) [46], which allows one to use either a local potential V (r) or directly the two-particle wave function, and additionally any source distribution as input to compute the correlation function. For the p-p correlation function the strong Argonne v 18 potential [50] in the S, P , and D waves is used as an input to CATS.
The theoretical correlation function for p-0 is modeled employing two different approaches. On the one hand, in CATS the correlation function is computed from the isospin-averaged wave functions obtained within a coupled-channel formalism. On the other hand, the Lednický-Lyuboshits approach [51] relies on the effective-range expansion using scattering parameters as input to evaluate the correlation function. The coupling of the n-+ system to p-0 considering the different thresholds is explicitly included by means of a coupled-channel approach, while the effect Table 1 Weight parameters for the individual components of the measured correlation function. Contributions from feed-down contain the mother particle listed as a subindex. Non-flat contributions are listed individually. The experimental data are compared with the modeled correlation function considering the finite experimental momentum resolution [27]. In addition to the genuine correlation function of interest, the measured correlation function also contains residual correlations due to protons coming from weak decays of other particles, such as and + (feed-down), and misidentifications. These effects are included by modeling the total correlation function as a decomposition, where the sum runs over all contributions. Their relative contribution is given by the λ parameters computed in a data-driven way from single-particle properties such as the purity and feed-down fractions [27], and are summarized in Table 1.
Apart from the genuine p-p correlation function, a significant contribution comes from the decay of particles feeding to the proton pair. The residual p-correlation function is modeled using either the Usmani potential [53], chiral effective field theory calculations at Leading (LO) [54], or Next-To-Leading order (NLO) [20]. The resulting correlation function is transformed into the momentum basis of the p-p pair by applying the corresponding decay matrices [55]. All other contributions are assumed to be C (k * ) ∼ 1. Due to the challenging reconstruction of the 0 , the experimental purity of the 0 sample is rather low, and additionally exhibits a strong dependence on the transverse momentum p T as demonstrated in Fig. 1. The average p T of the 0 candidates used to construct the correlation function at k * < 200 MeV/c, however, is lower than the p T of all inclusive 0 candidates. Considering this effect, the 0 purity employed to compute the λ parameters is found to be 27.4%. Accordingly, the main contribution to the p-0 correlation function stems from the combinatorial background appearing in the invariant mass spectrum around the 0 peak, which in the following is referred to as ( γ ). The shape of the p-( γ ) correlation function is extracted from the sidebands of the invariant mass selection, and is found to be independent of the choice of mass window. The non-flat behavior is mainly determined by residual p-correlations which are smeared by an uncorrelated γ , and defines the baseline of the measurement of the p-0 correlation function. The shape is parametrized with a Gaussian distribution and weighted by its λ parameter. All other contributions stemming from misidentified protons or from feeddown are assumed to be C (k * ) ∼ 1.
The total correlation function including all corrections is then multiplied by a polynomial baseline C non-femto (k * ), to account for the normalization and non-femtoscopic background effects stemming e.g. from momentum and energy conservation [27]. The p-p correlation function is fitted in the range The baseline is modeled as a polynomial of zeroth, first, and second order. Additionally, as discussed above, all three models for the p-residual correlation function are employed, and the input to the λ parameters is modified by ±20% while maintaining a constant sum of the primary and secondary fractions. The p-p correlation function is shown in Fig. 2, where the width of the bands corresponds to one standard deviation of the total systematic uncertainty of the fit. The inset shows a zoom of the p-p correlation function at intermediate k * , where the effect of repulsion becomes apparent. The femtoscopic fit yields a radius of r 0 = 1.249 ± 0.008 (stat) +0.024 −0.021 (syst) fm.
Analyses of π -π and K-K correlation functions at ultrarelativis-  [58]. The main resonances feeding to pions, ρ and ω, are significantly longer-lived than those feeding to protons ( ) and 0 ( (1405)). Hence, it is not surprising that the source distribution for π -π deviates from a Gaussian. These conclusions are underlined when fitting the p-p correlation function with a Lévy-stable source distribution [59,60].
Leaving both the femtoscopic radius and the stability parameter α for the fit to determine, the Gaussian source shape (α = 2) is recovered. Employing a Cauchy-type source distribution (α = 1), the data cannot be satisfactorily described. Therefore, the premise of a Gaussian source holds for baryon-baryon pairs. Accordingly, a Gaussian source with femtoscopic radius r 0 is used to fit the p-0 correlation function. The parameters of the linear baseline are obtained from a fit to the p-( γ ) correlation function in k * ∈ [250, 600] MeV/c, where it is consistent and kinematically comparable with p-0 , however featuring significantly smaller uncertainties. The experimental p-0 correlation function is then fitted in the range k * < 550 MeV/c, and varied during the fitting procedure within k * ∈ [500, 600] MeV/c to determine the systematic uncertainty. Additionally, the input to the λ parameters is modified by ±20% while maintaining a constant sum of the primary and secondary fractions. The parameters of the base-  EFT [20], NSC97f [26] and ESC16 [23], and using the Lednický-Lyuboshits approach [51,52] for fss2 [24]. The width of the bands corresponds to one standard deviation of the systematic uncertainty of the fit. The absolute correlated uncertainty due to the modeling of line are varied within 1σ of their uncertainties considering their correlation, including the case of a constant baseline. Finally, the femtoscopic radius is varied according to its uncertainties. Possible variations of the p-0 source due to contributions of m T scaling and strong decays are incorporated by decreasing r 0 by 15%, similarly as in [28,29]. The corresponding resonance yields are taken from the statistical hadronization model within the canonical approach [58].
All correlation functions resulting from the above mentioned variations of the selection criteria are fitted during the procedure, additionally considering variations of the mass window to extract the p-( γ ) baseline. The width of the bands in Fig. 3 corresponds to one standard deviation of the total systematic uncertainty of the fit. The absolute correlated uncertainty due to the modeling of the p-( γ ) baseline correlation function is shown separately at the bottom of the figure.

Results
The experimental p-0 ⊕ p-0 correlation function is shown in Fig. 3. The k * value of the data points is chosen according to the k * of the same event distribution N same (k * ) in the corresponding interval. Therefore, due to the low number of counts in the first bin, the data point is shifted with respect to the bin centre. Since the uncertainties of the data are sizable, a direct determination of scattering parameters via a femtoscopic fit is not feasible. Instead, the data are directly compared with the various models of the interaction. These include, on the one hand, meson-exchange models, such as fss2 [24] and two versions of soft-core Nijmegen models (ESC16 [23], NSC97f [61]), and on the other hand results of χ EFT at Next-to-Leading Order (NLO) [20]. The correlation function is modeled using the Lednický-Lyuboshits approach [51] considering the couplings of the p-0 system to p-and n-+ [52] with scattering parameters extracted from the fss2 model. For the case of ESC16, NSC97f and χ EFT, the wave function of the p-0 system, including the couplings, is used as an input to CATS to compute the correlation function. The degree of consistency of the data with the discussed models is expressed by the number of standard deviations n σ , computed in the range k * < 150 MeV/c from the p-value Table 2 Degree of consistency of the different models with the experimental correlation function.  Table 2 is computed as one standard deviation of the corresponding distribution. The data are within (0.2−0.8)σ consistent with the p-( γ ) baseline, indicating the presence of an overall shallow strong potential in the p-0 channel. The main source of uncertainty of the modeling of the correlation function is the parametrization of the p-( γ ) baseline due the sizeable statistical uncertainties of the latter. All employed models for the N-interaction potential succeed in reproducing the scattering data in the S = −1 sector [7].
Due to the available experimental constraints, the overall description of the p-interaction yields a consistent description. On the other hand, the corresponding p-0 correlation functions differ significantly among each other. This demonstrates that femtoscopic measurements can discriminate and constrain models, and therefore represent a unique probe to study the N-interaction. Both fss2 and χ EFT exhibit an overall repulsion in N-at intermediate k * , which mainly occurs in the spin singlet S = 0, I = 1/2 and spin triplet S = 1, I = 3/2 components [20,24]. In the low momentum region, below roughly 50 MeV/c, both models yield attraction, which is reflected in the profile of the correlation function. The Nijmegen models, on the other hand, are characterized by a rather constant attraction over the whole range of k * . In particular at low relative momenta, however, the behavior of the two models deviates significantly. The shape of the correlation function of the most recent Nijmegen model, ESC16, differs significantly from that of the other calculations. This is mainly due to the fact that the occurrence of bound states in the strangeness sector (S = −1, −2, −3) is not allowed in the model [23]. This leads to a repulsive core in all the N-channels, which can well be observed in Fig. 3 as the non-monotonic behavior at small relative momenta. In contrast to all other discussed models, NSC97f yields attraction in the spin triplet S = 1, I = 3/2 channel [61]. Accordingly, the corresponding correlation function demonstrates the strongest attraction at low momenta. The rather large differences among the modeled p-0 correlation functions demonstrate that the shape of the latter is very sensitive to details of the strong interaction, and driven by the interplay of the different spin and isospin channels. This shows the strength of femtoscopic measurements, in particular in the Nchannel. The underlying two-body N-interaction obtained within these models, however, translates into significantly different values for the in-medium single-particle potential U when included in many-body calculations. Both the fss2 quark-model, along with χ EFT, deliver similar results at nuclear saturation density, leading to an overall repulsive U of around [10][11][12][13][14][15][16][17]21,24]. This is in agreement with evidence from relativistic mean field calculations fitting experimental data of − atoms [12] and the experimental absence of bound states in hypernuclei [16]. On the contrary, both Nijmegen models yield a slightly attractive single-particle potential, ranging from ≈ −16 MeV for NSC97f to ≈ −3 MeV for ESC16. As already mentioned, however, the interpretation of hypernuclear measurements introduces a significant model dependence. This concerns not only the extraction of the experimental results, relying for instance on the framework of the distorted-wave impulse approximation [17], but also the extrapolation of theoretical calculations to finite density via e.g. the G-matrix approach [62,63].

Summary
This Letter presents the first direct investigation of the p-0 interaction in high-multiplicity pp collisions at √ s = 13 TeV, hence proving the feasibility of femtoscopic studies in the N-sector.
The p-0 correlation function is consistent with the p-( γ ) baseline, and therefore the measurement indicates the presence of an overall shallow strong potential. The data are compared with stateof-the-art descriptions of the interaction, including chiral effective field theory and meson-exchange models. Due to the scarce experimental constraints in the N-sector, the modeled correlation functions differ significantly among each other. The shape of the modeled correlation functions appears to be very sensitive to details of the strong interaction, and is driven by the interplay of the different spin and isospin channels. This proves that femtoscopic measurements in high-energy pp collisions provide a direct study of the genuine two-body N-strong interaction. The presented femtoscopic data cannot discriminate between different models, which is also the case for the available scattering and hypernuclei data. , will therefore shed light on the N-sector and provide constraints on the models describing the interaction.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements
The ALICE Collaboration is grateful to J. Haidenbauer and T. Rijken for valuable discussions and for providing the theoretical results for the p-0 interaction.
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector: A. I. Alikhanyan National Science Laboratory ( [50] R.B. Wiringa, V.G.J. Stoks, R. Schiavilla, An accurate nucleon-nucleon potential with charge independence breaking, Phys. Rev. C 51 (1995) 38-51.