Towards the understanding of the genuine three-body interaction for p$-$p$-$p and p$-$p$-\Lambda$

Three-body nuclear forces play an important role in the structure of nuclei and hypernuclei and are also incorporated in models to describe the dynamics of dense baryonic matter, such as in neutron stars. So far, only indirect measurements anchored to the binding energies of nuclei can be used to constrain the three-nucleon force, and if hyperons are considered, the scarce data on hypernuclei impose only weak constraints on the three-body forces. In this work, we present the first direct measurement of the p$-$p$-$p and p$-$p$-\Lambda$ systems in terms of three-particle correlation functions carried out for pp collisions at $\sqrt{s} = 13$ TeV. Three-particle cumulants are extracted from the correlation functions by applying the Kubo formalism, where the three-particle interaction contribution to these correlations can be isolated after subtracting the known two-body interaction terms. A negative cumulant is found for the p$-$p$-$p system, hinting to the presence of a residual three-body effect while for p$-$p$-\Lambda$ the cumulant is consistent with zero. This measurement demonstrates the accessibility of three-baryon correlations at the LHC.


Introduction
One of the open challenges of nuclear physics is the understanding of many-particle dynamics.Studies of the nuclear structure have unambiguously shown that calculations based only on nucleon-nucleon (N-N) interactions fail to accurately describe many experimental observables, such as nuclear binding energies along the periodic table of elements [1], the position of the neutron drip line for neutron-rich nuclei [2] or the properties of the recently observed four-neutrons resonance [3].A significant improvement in the modelling of nuclear bound objects has been achieved by including three-body forces in theoretical calculations.These three-body forces are implemented in chiral effective field theories [4] and in a number of ab initio many-body methods such as no-core shell model [5], coupled-cluster theory [6,7], self-consistent Green's function theory [8], similarity renormalisation group [9,10], and quantum Monte Carlo [11].Studies conducted on intermediate mass neutron-rich nuclei proved that the sensitivity to the three-body forces increases with the number of neutrons in the system [6].Three-body forces within light and medium-mass nuclei, where the nuclear saturation density corresponds to typical inter-particle distances of 2 fm, contribute about 10 − 15 % to the total interaction strength [12,13].However, at higher densities and shorter inter-particle distances their contribution might increase [2], but no data are available in such a regime and the properties of nuclear matter can be only extrapolated using the available information at saturation densities.The experimental information on the three-body forces involving Λ hyperons is even more scarce since the data available for hypernuclei are much less than the data for nuclei.Recent hypetriton measurements in several colliding systems at RHIC and LHC [14][15][16][17] provide important input to the understanding of N-N-Λ forces and future measurements will resolve the current tensions among the different estimations of the binding energy and life-time.Theoretical works assign to the hypertriton a radius of the order of 5 fm [18] and hence a N-Λ distance of 10 fm [13,19] within this state.Heavier hypernuclei are more compact and their size is comparable to that of normal nuclei so that they represent an optimal test bed for the N-N-Λ interaction [20].However, good fits of the theoretical models to the available hypernuclear data, from 7  Λ Li to 208 Λ Pb [21][22][23][24], require a full understanding of the shell-structure of such bound objects as well as accurate experimental constraints on the spin-dependent N-Λ interaction, in particular for the p-wave and higher partial waves.The lack of precise data as well as the limitations in the microscopic description of the structure of hypernuclei cause large ambiguities on the strength of the N-N-Λ three-body force.Further opportunities are provided by the recently observed 3  Λ n bound state [25] and planned experimental programs focused on neutron-rich hypernuclei [26].Nevertheless, the contribution from three-body forces in bound objects such as nuclei and hypernuclei cannot be separated from the lower-order two-body interactions, hence, complementary experimental methods to investigate three-baryon systems could provide an important contribution to this field.
Neutron-rich and dense baryonic matter constitutes an interesting system also because of its connections to the physics of neutron stars (NS) [27].The structure and composition of the innermost part of NS is not known.Amongst many possible scenarios, some models support the appearance of various hadronic particle species with increasing baryon density inside the star [27,28].The presence of hadronic degrees of freedom and their relative abundances are sensitive to the two-and three-body interaction models which are used to compute the equation of state (EoS) of NS matter.The different hypotheses can be tested by deriving the masses and radii of NS for a specific EoS and comparing them with the corresponding astrophysical observations [28].The suggestion of strange baryons inside NS is motivated by the fact that central densities of NS might become sufficiently large (ρ ≈ 3 − 4ρ 0 , where ρ 0 is the nuclear saturation density) to provide favourable conditions for the onset of strangeness production processes leading to, in particular, the formation of hyperons.The appearance of Λ hyperons in NS matter results in a softening of the EoS which is at variance with astrophysical observations of two solar mass stars [29,30].However, in Ref. [31], it was shown that by adding a strongly repulsive N-N-Λ interaction, tuned to reproduce the separation energies of Λ hyperons in several hypernuclei, a sufficiently stiff EoS can be obtained and even the massive NS observables can be reproduced.This indicates that three-body forces may have a silicon detectors placed radially between 3.9 and 43 cm around the beam vacuum tube.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 ionisation energy loss (dE/dx) in the TPC gas with up to 159 reconstructed space points along the particle trajectory.For high momentum particles, the TPC measurement is combined with information provided by the Time-Of-Flight (TOF) [56] detector system, which 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.
The primary vertex (PV) of the event is reconstructed with the combined track information of the ITS and the TPC, and independently with track segments in the two innermost layers of the ITS.The reconstructed PV of the event is required to have a maximal displacement with respect to the nominal interaction point of 10 cm along the beam axis, in order to ensure a uniform acceptance.Pile-up events with multiple primary vertices are removed following the procedure described in Refs.[39,43,57].This rejects the events with pile-up of collisions occurring in the same or nearby bunch crossings.However, additional clean-up has to be applied on the track selection level to reject particles produced in pile-up collisions in the long TPC readout time.
A total of 1.0 × 10 9 HM events are used for the analysis after event selection.In order to build the three-particle correlation functions of p-p-p and p-p-Λ systems, particle and antiparticle distributions are combined.In the following, p-p-p refers to p-p-p ⊕ p -p -p and p-p-Λ refers to p-p-Λ ⊕ p -p -Λ.The proton and Λ candidates as well as their antiparticles need to be selected.As the particle and antiparticle selections are identical, only the particles are explicitly discussed below.Both particle species are reconstructed using the procedure described in Ref. [57], while the related systematic uncertainties are evaluated by varying the kinematic and topological selection criteria used in the reconstruction.In the following text, the systematic variations are enclosed in parentheses.
The primary protons are selected in the momentum interval 0.5 (0.4, 0.6) < p T < 4.05 GeV/c and |η| < 0.8 (0.77, 0.85).To improve the quality of the tracks a minimum of 80 (70,90) out of the 159 possible spatial points inside the TPC are required.The PID selections are applied by comparing the measured dE/dx and time-of-flight with the expected values for a proton candidate.The agreement is expressed in multiples (n PID σ ) of the detector resolution σ .For protons with p T < 0.75 GeV/c the n PID σ is evaluated only based on the specific energy loss in the TPC, while for p T ≥ 0.75 GeV/c a combined TPC and TOF PID selection is applied n PID σ = n 2 σ ,TPC + n 2 σ ,TOF .The n PID σ of the accepted proton candidates is required to be lower than 3 (2.5, 3.5).To reject particles that are non-primary or come from pile-up collisions, the distance of closest approach (DCA) to the PV of the tracks is required to be less than 0.1 cm in the transverse plane and less than 0.2 cm along the beam axis.The purity of candidates is estimated using Monte Carlo (MC) simulations by taking the ratio of the number of reconstructed true protons produced by the generator and the number of all candidates identified as protons as a function of the reconstructed transverse momentum.The contributions of secondary protons stemming from weak decays of strange baryons and from interactions in the detector material are extracted using MC template fits to the measured distributions of the DCA to the PV [39].The average purity of the identified protons is 98.3% and 86.6% of them are primaries.
The Λ candidates are reconstructed via the weak decay Λ → pπ − (the Λ → pπ + in case of Λ reconstruction).The secondary daughter tracks are selected with similar criteria as for the primary protons regarding |η| and the number of hits in the TPC.However, a less strict PID requirement of n PID σ < 5(4) is used.In addition, the daughter tracks are required to have a DCA to the PV of at least 0.05 (0.06) cm and the DCA between the daughter tracks at the secondary vertex must be smaller than 1.5 (1.2) cm.The cosine of the pointing angle (CPA) between the vector connecting the PV to the decay vertex and the 3-momentum of the Λ candidate is required to be larger than 0.99 (0.995).To reject unphysical secondary vertices, reconstructed with tracks stemming from pile-up of pp collisions occurring in different bunch crossings, the decay tracks are required to possess a hit in the two innermost or the two outermost ITS layers or a matched TOF signal [42].Finally, a selection on the candidate invariant mass (IM) is applied by requiring it to be in a ± 4 MeV/c 2 interval around the nominal Λ mass [58].The primary and secondary contributions to the yield of Λ are extracted employing a similar method as for protons but using the CPA as an observable for the template fits.The Λ hyperons produced in primary interactions contribute to about 58.5% of their total yield.About 19.5% originate from the electromagnetic decays of Σ 0 .The number of Σ 0 particles is related to their ratio to the Λ hyperons, which is fixed to 1/3 based on predictions from the isospin symmetry and a measurement of the corresponding production ratios [59].Further, each of the weak decays of Ξ − and Ξ 0 contributes about 11 % to the yield of Λ hyperons.The purity of Λ and Λ has been extracted by fitting the IM spectra of candidates as a function of the threeparticle kinematic variable Q 3 which is defined in Eq. 5.The fits have been performed in the IM range of 1090 to 1150 MeV/c 2 using a double Gaussian for the Λ signal and a second-order polynomial for the background.The result has been averaged for Q 3 < 1 GeV/c, leading to a combined purity of Λ and Λ of 95.6%.
The systematic uncertainties are evaluated by performing simultaneous variations of the selection criteria for protons and Λ candidates as well as for the corresponding antiparticles.The variations are randomly combined in 44 sets in which at least one of the selection criteria is varied.Such procedure allows to account for the correlations between the systematic uncertainties.Each random set of variations is accepted for the evaluation of the systematic uncertainties only if the yield of the triplets is varied by less than 10% with respect to the standard selection in the kinematic region Q 3 < 0.4 GeV/c.

Three-particle correlation function
The observable of interest in femtoscopy is usually the two-particle momentum correlation function [35,60], which is defined as the probability to simultaneously find two particles with momenta p 1 and p 2 divided by the product of the corresponding single particle probabilities These probabilities are related to the inclusive Lorentz-invariant spectra In the absence of a correlation signal, the value of C(p 1 , p 2 ) is constant and normalised to unity.A similar logic can be followed to construct the three-particle correlation functions as Following [61,62], Eq. 1 can also be written as where S(r * ) is the distribution of the relative distances of particle pairs in the pair rest frame (PRF, denoted by the * ) -the so-called source function.The properties of the source in pp collisions at √ s = 13 TeV have been evaluated in Ref. [57], including the effects of short-lived resonance decays which enlarge the effective source size.The wave function of the particle pair relative motion is denoted by ψ(r * , k * ) where k * = (p * 1 − p * 2 )/2 is the relative momentum.The wave function encapsulates the details of the particle interaction and drives the shape of the correlation function.In case of the threeparticle correlation function, the two-particle source function and the wave function of the particle pair relative motion must be replaced by a three-particle source function and wave function.In this analysis, the measured three-particle correlation functions are not compared to theoretical predictions.The goal here is to extract the three-particle femtoscopic cumulants which provide experimental evidence of the existence, or the absence, of genuine three-particle correlations, as explained in Section 2.3.
The three-particle correlation function can be written as where N s (Q 3 ) and N m (Q 3 ) are the same-event and mixed-event distributions of three particle combinations (triplets) as a function of Q 3 and N is the normalisation parameter.The Lorentz-invariant variable Q 3 is defined in [48] as where q i j is the norm of the four-vector [35] which can be rewritten as Here m i and m j are the particle i and j masses, p µ i and p µ j are the particle four momenta, while q µ i j is the relative four-momentum of the pair i j.In the case of same mass particles, the term where k * i j is the relative momentum of the i j pair in the PRF.The mixed-event sample is obtained using event-mixing techniques, in which the particle triplets of interest are generated by combining single particles stemming from three different events.To maintain the same acceptance effects as in the same event sample, the mixing procedure is conducted only for events with similar z position of the primary vertex and multiplicity [39].Additionally, in order to correct for possible differences in terms of multiplicity distribution between same and mixed events, the yield of the latter is re-weighted in each multiplicity interval to have the same statistical weight as the distribution when particles are from the same event.To account for the two-track merging and splitting effects due to the finite two-track resolution in the same-event sample, a minimum value of the distance between two proton tracks (in case of p-Λ pairs, the proton from Λ decay is considered along with the primary proton) on the azimuthal-polar angles plane ∆η-∆ϕ is applied to both the same-and mixedevent samples.The default selection is ∆η 2 + ∆ϕ 2 ≥ 0.017 2 and a systematic variation of +10 % for the value of the minimum distance is applied in the analysis.The normalisation parameter N is chosen such that the mean value of the correlation function equals unity in a Q 3 region where the effects of FSIs are negligible.The interval Q 3 ∈ (1.0 − 1.2) GeV/c is chosen for all triplets.

Three-particle femtoscopic cumulants
The measurable three-particle correlation function C(p 1 , p 2 , p 3 ) include all interactions at work in the three-particle system: the two-body interactions among all pairs within the selected triplet and the genuine three-body interaction.To access only the genuine three-body correlations, one can use cumulants.Given random variables X i , the cumulant for a triplet is defined by Kubo [50] as where ⟨X i ⟩ is the expectation value of the variable X i and X i X j , X i X j X k are the two-and three-variable joint moments.The three-particle correlation function, defined in Eq. 4, is the three-particle momentum distribution normalised to the mixed-event distribution.The cumulants method can be applied to the numerator which contains the correlated particles, and then the expression is normalised to the mixedevent distribution.The three-particle femtoscopic cumulant c 3 thus can be defined as where N 3 (p 1 , p 2 , p 3 ) and N 2 p i , p j are the same-event three-and two-particle momentum distributions; N 1 (p i ) is the single-particle momentum distribution; the product terms the mixed event distributions.Thus one can further rewrite the femtoscopic cumulant as This method has been already successfully applied within the ALICE Collaboration to study the possibility of coherent pion production by measuring three-pion femtoscopic cumulants in Ref. [48,49].Theorem I from Ref. [50] enunciates that the three-particle cumulant is zero if the variables X i , X j , ... can be divided into two or more groups that are statistically independent.In case of femtoscopic cumulants, this translates into c 3 (p 1 , p 2 , p 3 ) = 0 in the absence of genuine three-body correlations.Therefore, the measurements of non-vanishing values of c 3 can be used as an experimental confirmation of the existence of genuine three-body effects.
If genuine three-body correlations are not present in the particle triplet, the three-particle correlation function can be expressed using only lower order contributions as follows In Eq. 11, C([p i , p j ], p k ) is built by combining particles i and j from the same event with particle k from another event to obtain the numerator N 2 p i , p j N 1 (p k ) of the correlation function while the denominator (N 1 (p 1 ) N 1 (p 2 ) N 1 (p 3 )) is estimated using three particles from three different events as described in Section 2.2.

Projector method
An alternative method to isolate the genuine three-body contribution to the measured three-particle correlation functions is the projector method [51].This method makes use of the subtraction rule provided by the Kubo's cumulant decomposition (Eq.10) but, instead of evaluating them with the data-driven approach based on event mixing described above, it calculates C([p i , p j ], p k ) using the measured or the calculated two-particle correlation function and the projection of the third non-interacting (spectator) particle.The method is described in Ref. [51].Given the three-particle correlation function, C(Q 3 ), and the two-body correlation functions, C(k * i j ), the projector method provides a kinematic transformation from the relative momentum k * i j of the interacting pairs i j to the Q 3 of the three-body system (i − j) − k under study.The transformation is given by the following integral in the momentum space where the indices i j denote the interacting pair and the projector function W i j is equal to [51] Three-baryon interaction

ALICE Collaboration
The constants α, β and γ depend on the particle masses 1 .The integral in Eq. 12 can be evaluated us- The upper panels show the comparison of the two-particle correlations projected on three-particle phase space obtained using the data-driven approach based on event mixing (green points) and the projector method (grey band).The resulting correlation functions are shown for (p-p)-p (panel a), (p-p)-Λ (panel b) and p-(p-Λ) (panel c) cases.The error bars and the boxes represent the statistical and systematic uncertainties, respectively.The grey band includes systematic and statistical uncertainties summed in quadrature.The lower panels show the deviations between the data-driven approach and the projector method, expressed in terms of n σ .ing the measured p-p and p-Λ correlation functions from Refs.[41,57,63].The resulting correlation functions are compared to the ones obtained by employing the data-driven method (Eq.11) and shown in Fig. 1.Panel a) shows the (p-p)-p correlation function, the green points are obtained using the datadriven approach and the grey band is obtained with the projector method.The statistical and systematic uncertainties are shown separately for the data driven method, while the width of the grey band represents the sum in quadrature of the statistical and systematic uncertainties for the projector method.The statistical uncertainties of all the measured correlation functions have been estimated using a bootstrap [63] method by sampling same-and mixed-event counts from Poisson distributions.The statistical uncertainties shown correspond to the central 68% confidence interval and are consistent with the uncertainties obtained employing the standard error propagation method.The systematic uncertainties are estimated by varying the selection criteria of the particle candidates as described in Section 2.1.Panels b) and c) show the same comparison for the (p-p)-Λ and the p-(p-Λ) correlation functions.The number of events used for mixing to obtain (p-p)-p and (p-p)-Λ correlation functions is 30.The numerator of p-(p-Λ) correlation function requires p-Λ pairs in same event sample, which are less abundant than p-p pairs.Thus to obtain good statistical precision, the number of events used for mixing (to account for the third uncorrelated particle) must be increased to 100 in case of p-(p-Λ) triplets.
The results from the data-driven and the projector method are in good agreement between each other.The number of deviations n σ in each bin are shown in the bottom panels of Fig. 1, where σ is the combined statistical and systematic uncertainty for both the experimental data and the projector.The agreement in the region Q 3 < 0.8 GeV/c has been evaluated by performing a χ 2 test.The χ 2 is calculated combining the n σ values of each bin.Finally, the p-value from the χ 2 -distribution is computed and the global n σ values are extracted.The latter amount to 0.167, 0.0006 and 2.75 for (p-p)-p, (p-p)-Λ and p-(p-Λ), respectively.The data-driven method requires the usage of the third particle in the triplet from the mixedevent data sample and consequently the statistical uncertainty depends on the number of events used for mixing, while the projector method does not have this limitation.Thus, the latter significantly reduces the total uncertainty in the evaluation of the two-particle correlation effect on the three-particle correlation functions.For this reason, the projector method is used to calculate the three-particle cumulants for the p-p-p and p-p-Λ triplets.
The total two-particle contribution to the three-particle correlation function is obtained by substituting all terms on the right-hand side of Eq. 11 with the corresponding kinematic transformation, i.e.
where the indices refer to the label of the correlated pairs.In the case of p-p-p we have and in the case of p-p-Λ we have The resulting total lower-order contributions to the three-particle correlation functions (Eqs.15 and 16) are shown in Fig. 2. The agreement between the data-driven approach and the projector method predictions translate into n σ = 0.167 and n σ = 0.0014 for the p-p-p and p-p-Λ lower-order contributions, respectively.Comparison of the total two-particle contribution to the three-particle correlation functions obtained using the data-driven approach (green points) and the projector method (grey band).The resulting correlation functions are shown for p-p-p (left panel) and p-p-Λ (right panel).The error bars and the boxes represent the statistical and systematic uncertainties, respectively.The grey band includes systematic and statistical uncertainties summed in quadrature.

Decomposition of the three-particle cumulants
The experimental determination of the correlation function is mainly distorted by two distinct impurities in the candidate sample: misidentified particles and feed-down particles originating from weakly decaying particles.This introduces additional contributions to the correlation function of interest.These contributions are either assumed to be flat or, when the interaction is known, they are explicitly considered as discussed in Ref. [39].The contributions to the correlation function stemming from decaying particles or impurities of the sample are weighted with the so-called λ parameters.By adopting this technique the residual correlations can be included in the final description of the experimental correlation function of two particles as where the i j ̸ = 00 denote all possible impurity and feed-down contributions and the i j = 00 is the correctly identified primary particle contribution.These λ parameters are obtained employing single particle properties such as the purity and feed-down probability.The underlying mathematical formalism is outlined in Ref. [39].This mechanism has been extended to the three-particle case and the genuine three-particle cumulants can be obtained by subtracting the impurity and feed-down contributions from the measured cumulants.The full mathematical derivation is presented in Appendix C. The final expression of the genuine three-particle cumulants is where X, Y and Z represent three generic particle species, the index 0 refers to correctly identified primary particle and the indexes i, j, k refer to misidentified or to secondary particles of a generic particle species.As shown in Appendix C, the specific weights λ depend on the purity and feed-down fraction of the single particles and are found to be equal to λ X 0 Y 0 Z 0 (ppp) = 0.618 and λ X 0 Y 0 Z 0 (ppΛ) = 0.405 for the p-p-p and p-p-Λ cumulants, respectively.Only 60% (40%) of the p-p-p (p-p-Λ) triplets correspond to correctly identified primary particles.
In the following, the results for the p-p-p cumulants will be corrected according to the evaluated λ parameters assuming that all the three-particle contributions stemming from feed-down and impurities are flat in the momentum space.This assumption is supported by the observation that the measured p-p-Λ cumulants are consistent with zero within uncertainties (see Fig. 4 and the discussion in Section 3).The correction is not applied to the p-p-Λ cumulants because the shape of the feed-down contribution is not known and also because the statistical uncertainties are too large to provide any sensitivity to the three particle correlations.

Results
The measured three-particle correlation functions for p-p-p and p-p-Λ triplets are shown in Fig. 3 on the left and right panels, respectively.The number of events used for mixing for both cases is 30.The total number of same event triplets that are present at the range Q 3 < 0.8 GeV/c are 17840 for p-p-p, 10980 for p -p -p , 9191 for p-p-Λ and 5886 for p -p -Λ.The green symbols represent the data points with their statistical and systematic uncertainties, while the grey bands correspond to the lowerorder two-body interaction contributions obtained using the projector method already shown in Fig. 2.
The non-femtoscopic contributions to the measured correlation functions, evaluated using Monte Carlo simulations, are found to be negligibly small (see Appendix A for a detailed discussion).The green points show the experimental results, the error bars and the boxes represent the statistical and systematic uncertainties, respectively.The grey bands represent the expectations for the lower-order two-particle correlations obtained using the projector method and the band width is obtained including systematic and statistical uncertainties summed in quadrature.
In the low Q 3 region, the measured correlation functions deviate from the projected lower order contributions obtained using only two-particle correlations.The genuine three-body effects are then isolated by evaluating the cumulants The lower-order contribution C two-body (Q 3 ) obtained with the projector method is used.The results for p-p-p and p-p-Λ triplets are shown in Fig. 4 on the left and right panels, respectively.The p-p-p cumu-  3. The p-p-p cumulant in the left panel is further corrected for the feed-down contributions from decaying particles and represents thus, the cumulant for the correctly identified primary protons (see Section 2.5 for details).The dashed lines correspond to the assumption that there are no genuine three-body correlations c 3 (Q 3 ) = 0.The red open circles in left panel represent the cumulant for p-p-p triplets (for more details see the main text).lant, already corrected for the feed-down contributions, is negative for 0.16 < Q 3 < 0.22 GeV/c, while the large statistical uncertainty in the lowest Q 3 interval prevents a conclusion on the sign for Q 3 < 0.16 GeV/c.The agreement between the measured cumulant and the assumption that there are no genuine three-body effects is evaluated using the χ 2 test in the region Q 3 < 0.4 GeV/c, where the two-body interactions are prominent.There is no theoretical or experimental knowledge on the exact Q 3 range where three-body effects become relevant, however they are expected to contribute at lower or same Q 3 values as the two-body interactions.For this reason, the region of two-body forces was chosen.The obtained p-value corresponds to 6.7 standard deviations.If the cumulant is obtained using data-driven method to estimate lower order contributions, it results in 6.0 standard deviations.This result hints to the presence of an effect beyond the two-body interactions in the p-p-p system that could be either due to Pauli blocking (Fermi-Dirac quantum statistics) at the three particle level [64] or to the contribution of a three-body nuclear repulsive interactions.Long-range Coulomb interactions may also lead to significant contributions [65].More quantitative conclusions on the interpretation of the non-zero cumulant require more sophisticated calculations for the three-body system.The present analysis demonstrates the experimental accessibility of the three-baryon cumulant in the data sample of pp collisions at √ s = 13 TeV recorded by ALICE.In addition to the p-p-p system, the mixed-charge p-p-p case has also been studied (see Appendix B for details).In the case of p-p-p , the two-body p-p interaction contains both elastic and inelastic components, previously measured in an independent analysis [46]; the p-p interaction is well known from scattering data [66,67] and verified by correlation measurements [39,41,57]; and the genuine three-body strong interaction for the p-p-p triplet should be negligible.The extracted cumulant for p-p-p is shown in the left panel of Fig. 4 by the red open circles.Since the number of the mixedcharge triplets is a factor four higher than the one of the same-charge triplets, it is possible to extend the measurement of the three-particle correlation to lower Q 3 values.The correction for the feed-down contributions has been applied and the statistical and systematic uncertainties are shown.The p-p-p cumulant evaluated using the projector method (data-driven approach) agrees with the assumption of only two-body correlations present in the system within 2.1 (2.2) standard deviations in the region Q 3 < 0.4 GeV/c and within 0.9 (0.9) standard deviations in the region Q 3 < 0.2 GeV/c, suggesting that genuine three-body effects are not statistically significant.This result as well demonstrates that the measured p-p-p cumulant deviation from zero is not due to detector effects.
In the case of p-p-Λ, a positive cumulant is measured at Q 3 < 0.16 GeV/c.The p-value obtained from the χ 2 test in the region Q 3 < 0.4 GeV/c corresponds to a deviation of 0.8 σ from the assumption that no genuine three-body correlations are present.Using data-driven method to estimate lower order contributions results in 0.8 σ as well.A similar value is found by repeating the significance test with the Fisher method, meaning that the measured cumulant is compatible with the assumption of no genuine three-body effects within the uncertainties.The current measurement does not allow to draw any firm conclusion yet on the three-body interaction in the p-p-Λ system, but since in this case only two of the particles are identical and charged, a non zero cumulant can be directly linked to the presence of a strong three-body interaction.It is estimated that employing a three-baryon event filtering during the upcoming Run 3 data taking should increase the number of triplets by a factor up to 500 for the target integrated luminosity of 200 nb −1 at √ s = 13.6 TeV [68].This opens up the possibility of measuring precisely the three-body correlations for both the p-p-p and p-p-Λ systems.

Conclusions
In this article, the first femtoscopic study of the p-p-p and p-p-Λ systems measured in high-multiplicity pp collisions at √ s = 13 TeV with the ALICE detector has been presented.In the chosen colliding system, hadrons are emitted at average relative distances of about 1 fm providing a unique environment to test three-body interactions at scales shorter than inter-particle ones in nuclei.The data collected during the LHC Run 2 enabled the measurement of the p-p-p and p-p-Λ correlation functions in the low Q 3 region down to 0.1 GeV/c, giving access to the region where the effects of the hadronic two-and threebody interactions are more pronounced.The genuine three-particle correlations have been isolated using the Kubo's cumulant expansion method.The lower-order two-body contributions have been estimated employing both a data-driven event mixing technique and a newly developed projector method.The two approaches have been compared and found to be in good agreement between each other, providing the first validation of the projector method using the data.The extracted p-p-p and p-p-Λ cumulants deviate from zero in the low Q 3 region.In the case of p-p-p, a negative three-particle cumulant is measured.The p-value extracted from the χ 2 test corresponds to a deviation of 6.7σ from the assumption of only two-body correlations present in the system for Q 3 < 0.4 GeV/c.The obtained result provides an experimental hint to the presence of an effect beyond pairwise interactions in the p-p-p system.The observed deviation could be due to genuine three-body effects arising from: Pauli blocking, shortrange strong interactions, or long-range Coulomb interactions.Refined three-body system calculations are required to give a solid interpretation of the measurement.The mixed-charge p-p-p cumulant has also been measured as a benchmark and the result is consistent with the assumption that only twobody correlations are present in the system showing that the effect observed for the p-p-p system is a genuine one.In the case of p-p-Λ, a positive cumulant is observed at low Q 3 .The deviation from zero at Q 3 < 0.4 GeV/c is 0.8 σ , which suggests no significant deviation from the assumption that only two-body correlations are present in the system within the current uncertainties.For this system, where one particle is uncharged and only two particles are identical, genuine three-particle correlations can be directly linked to the three-body strong interaction.The upcoming LHC Run 3 data taking will provide significantly larger samples of measured triplets, allowing more quantitative conclusions to be drawn for many-body dynamics.
The analysis presented in this article represents a first important step towards the direct measurement of the three-body interaction among baryons, demonstrating that genuine three-particle effects can be studied using three-particle correlation functions as experimental observables.

A Monte Carlo studies
One of the benchmarks of the presented analysis consists in verifying that the three-particle correlations obtained from PYTHIA 8 [69] (Monash 2013 Tune) simulations do not show any significant deviation from unity.Indeed, no FSIs -either two-or three-particle -are included in the simulated and reconstructed Monte Carlo data using the PYTHIA 8 event generator for pp collisions at √ s = 13 TeV, thus Monte Carlo can neither be used to estimate the lower-order contributions (two-body correlations) nor genuine three-body effects.Figure A.1 shows the comparison between the measured and simulated correlation functions as a function of Q 3 , where the simulation includes a dedicated high-multiplicity selection to mimic the V0 high-multiplicity trigger in the real data.The green symbols represent the experimental data while the black symbols refer to the simulation.The simulated correlation functions are consistent with unity for the entire Q 3 < 0.8 GeV/c range, showing that there are no effects caused by the track reconstruction in the detector, as well as no sign of mini-jets contribution (see more details in Ref. [46]) and that the energy and momentum conservation effects are eventually present at larger values of Q 3 that are not relevant for the studies carried out in this work.Also the simulations of the lower-order contributions display the same behaviour.
The comparison between the correlation functions obtained from the measurements (green) and from the PYTHIA 8 event generator with a dedicated high-multiplicity selection to mimic the V0 high-multiplicity trigger (black).

B Mixed-charge correlation studies
An additional benchmark for the p-p-p cumulant result and the measured deviation from zero has been considered by studying the p-p-p triplets.Identical event and track selection criteria and systematic variations as those employed for the p-p-p analysis have been used.The rejection in the ∆η-∆ϕ plane (see Section 2.2 for the details) has been applied only for same-charge pairs in the triplet.1.This is consistent with the expectation from Eqs. 12 and 13, since such correlation functions depend only on the correlation of the two particles in the same event and on the mass of the uncorrelated particle, which are identical in the (p-p)-p and (p-p)-p systems.The (p-p)-p correlation function reflects the interplay of FSIs and non-femtoscopic correlations measured in the study of p-p pairs [46].The grey band in panel c) of Fig. B.1 is obtained by using as input of Eq. 12 the correlation function C(k * ) of p-p pairs emitted in p-p-p triplets with Q 3 < 1 GeV/c.Such selection is performed in order to use p-p pairs produced in similar shape events as those where low Q 3 triplets are found.The requirement is necessary for p-p correlations due to the mini-jet contribution [46] which is instead not present in p-Λ and p-p (see e.g.[70]).The systematics induced by this additional selection are considered in the uncertainties.
The p-p-p cumulant is obtained using the projector method and it is shown in Fig. 4 by red circular symbols.The result is compatible with zero within the statistical and systematic uncertainties, demonstrating that strong-interaction as well as Coulomb effects on three-particle level are not statistically significant in p-p-p .

C Feed-down contributions
The measured femtoscopic correlations do not originate only from correctly identified primary particles but they include as well contributions from misidentified particles and feed-down particles originating from weakly decaying hadrons.In case of two-body femtoscopy, the decomposition method explained in Section 2.5 and the Eq. 17 can be used to account for such impurities and feed-down effects.This method can be extended to the three-particle case.The total data sample X contains the particles which stem from feed-down X F , misidentified particles X M and the correctly identified primary particles X 0 .Both feed-down and misidentified particles can originate from different channels and the contributions can be expressed as where N F and N M are the number of feed-down and misidentification contributions.The purity is the fraction of correctly identified particles to the total number of particles in the data sample and can be defined as The correctly identified particles can stem from the decays of particles and for this purpose the channel fraction f (X i ) is defined as The fraction of specific channel in the whole data sample then can be written as The correlation function for three particles can be written as where N and M are the yields of XY Z triplet in same and mixed events, respectively.Using the identities introduced before one can write The correlation function then can be rewritten as where C i, j,k (XY Z) is the correlation function of the i, j, k-th channel of origin of the particles X,Y, Z and the λ i, j,k (XY Z) is the weight for such correlation.This parameter depends only on the mixed event sample and can be related to previously introduced single particle quantities, channel fraction and purity, as follows To study the lower-order contributions in the measured three-particle correlation function, one needs to define the decomposition for the measurement of two correlated particles.In such case, the origin of the third particle in the numerator does not matter as it is uncorrelated.Equation C.9 becomes C(XY, Z) = ∑ i, j N (X i Y j , Z) M(XY, Z) N (X i Y j , Z) M (X i Y j , Z) C i, j (XY,Z) M (X i Y j , Z) M(XY, Z) Here C i, j (XY, Z) denotes the correlation function where the two correlated particles X and Y are from origin i and j respectively and Z is from any origin.Here λ i, j (XY, Z) is (C.12) This notation is valid only if one wants to study the (X −Y ) − Z correlation.To obtain the cumulant of the primary particles which were identified correctly, one has to subtract the lower-order correlations, such as (X −Y ) − Z, from the three-particle correlation.For this purpose, the Eq.C.11 must be rewritten.As previously explained, in case of (X −Y ) − Z correlation, the origin of the third uncorrelated particle is not important, which means that C(XY, Z l ) = C(XY, Z m ), where Z l and Z m have different origin.As previously shown, the fraction of particles Z l in the whole data sample is P(Z l ), which can be as well expressed with the λ parameter of one particle λ l (Z).Using the property 1 = ∑ k λ k (Z) of the λ parameters in Eq.C.12, one can write (C.16) Written in such a way, one can already see that the cumulant of the measured correlation function can be expressed as the sum of the correctly identified primary particle cumulant and the cumulant of the rest of the contributions as follows (C.17) The terms inside the brackets are almost a cumulant expression, except the +2 term is missing, but one can add and subtract 2 to obtain λ i, j,k (XY Z)c(X i Y j Z k ).

(C.18)
This is the final result -the cumulant calculated using the measured correlation functions consists of the three correctly identified primary particle cumulant and the cumulant which consist of the rest of possible contributions.In such case, the correctly identified particle cumulant is λ i, j,k (XY Z)c(X i Y j Z k ) .

Figure 2 :
Figure2: Comparison of the total two-particle contribution to the three-particle correlation functions obtained using the data-driven approach (green points) and the projector method (grey band).The resulting correlation functions are shown for p-p-p (left panel) and p-p-Λ (right panel).The error bars and the boxes represent the statistical and systematic uncertainties, respectively.The grey band includes systematic and statistical uncertainties summed in quadrature.

Figure 3 :
Figure3: Measured p-p-p (left panel) and p-p-Λ (right panel) three-particle correlation functions.The green points show the experimental results, the error bars and the boxes represent the statistical and systematic uncertainties, respectively.The grey bands represent the expectations for the lower-order two-particle correlations obtained using the projector method and the band width is obtained including systematic and statistical uncertainties summed in quadrature.

Figure 4 :
Figure4: Three-particle cumulants for p-p-p (left panel, blue square symbols) and p-p-Λ (right panel) triplets obtained by subtracting the lower-order contributions from the measured three-particle correlation functions shown in Fig.3.The p-p-p cumulant in the left panel is further corrected for the feed-down contributions from decaying particles and represents thus, the cumulant for the correctly identified primary protons (see Section 2.5 for details).The dashed lines correspond to the assumption that there are no genuine three-body correlations c 3 (Q 3 ) = 0.The red open circles in left panel represent the cumulant for p-p-p triplets (for more details see the main text).
The obtained correlation functions for p-p-p , (p-p)-p and (p-p)-p triplets with the corresponding statistical and systematic uncertainties are shown in panels a), b) and c) of Fig. B.1, respectively.The grey bands are obtained using the projector method following the procedure described in Section 2.4.The (p-p)-p correlation functions obtained with the data-driven approach and the projector method (panel b) in Fig. B.1 are in agreement with the (p-p)-p results (panel a) in Fig.

Figure B. 1 :
Figure B.1: Panel a) shows the correlation function for p-p-p triplets (green data points) and the total lower order contributions (grey band).Panels b) and c) show the (p-p)-p and (p-p)-p lower-order contributions to the measured p-p-p correlation function.The error bars and the boxes represent the statistical and systematic uncertainties, respectively.The grey band includes systematic and statistical uncertainties summed in quadrature obtained from the projector method.