Study of the $\Lambda$-$\Lambda$ interaction with femtoscopy correlations in pp and p-Pb collisions at the LHC

This work presents new constraints on the existence and the binding energy of a possible $\Lambda$-$\Lambda$ bound state, the H-dibaryon, derived from $\Lambda$-$\Lambda$ femtoscopic measurements by the ALICE collaboration. The results are obtained from a new measurement using the femtoscopy technique in pp collisions at $\sqrt{s}=13$ TeV and p-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV, combined with previously published results from p-Pb collisions at $\sqrt{s}=7$ TeV. The $\Lambda$-$\Lambda$ scattering parameter space, spanned by the inverse scattering length $f_0^{-1}$ and the effective range $d_0$, is constrained by comparing the measured $\Lambda$-$\Lambda$ correlation function with calculations obtained within the Lednicky model. The data are compatible with hypernuclei results and lattice computations, both predicting a shallow attractive interaction, and permit to test different theoretical approaches describing the $\Lambda$-$\Lambda$ interaction. The region in the $(f_0^{-1},d_0)$ plane which would accommodate a $\Lambda$-$\Lambda$ bound state is substantially restricted compared to previous studies. The binding energy of the possible $\Lambda$-$\Lambda$ bound state is estimated within an effective-range expansion approach and is found to be $B_{\Lambda\Lambda}=3.2^{+1.6}_{-2.4}\mathrm{(stat)}^{+1.8}_{-1.0}\mathrm{(syst)}$ MeV.


Introduction and physics motivation
A detailed characterization of the Λ-Λ interaction is of fundamental interest since it plays a decisive role in the quantitative understanding of the hyperon (Y) appearance in dense neutron-rich matter, in protoneutron and in neutron stars [1]. If hyperons do appear at large densities and their fraction becomes sizeable, the Y-Y interaction is expected to play an important role in the equation of state of the system [2,3]. Even if the hyperon densities in compact objects are negligible, the interplay between the average separations and the Λ-Λ effective range determine the possible onset of phenomena such as fermion superfluidity, and hence influence the transport properties of the system [4][5][6].
The characterization of the Λ-Λ interaction is still an open issue in experimental nuclear physics. The Nagara event, recently measured with the emulsion technique [7,8], reports a clear evidence for a double-Λ hypernucleus 6 ΛΛ He, with a small binding energy between the two Λs of ∆B ΛΛ = 0.67 ± 0.17 MeV. This value was obtained by comparing the binding energy of the two Λs inside the double hypernucleus (B ΛΛ = 6.91 ± 0.16 MeV) with the binding energy of a single Λ in a single-hypernucleus, however, it might be influenced by three-body forces. Nevertheless, this result was used to set a lower limit for the mass of the predicted but so far not observed H-dibaryon, a possible bound state composed of six quarks (uuddss) [9]. Several experimental collaborations have been involved in the search for this state in the decay channels H → Λpπ and H → ΛΛ, in nuclear and elementary (e − e + ) collisions, but no evidence has been found [10][11][12], even though an enhanced Λ-Λ production near threshold was measured by E224 and E522 at KEK-PS [13,14].
Theoretical models constrained to the available nucleon-nucleon and hyperon-nucleon experimental data, assuming either a soft [15][16][17] or a hard [18,19] repulsive core for the Λ-Λ interaction, predict different scattering lengths ( f 0 ) and effective ranges (d 0 ). Throughout this paper the standard sign convention in femtoscopy is used, according to which a positive f 0 corresponds to an attractive interaction, while a negative scattering length corresponds either to a repulsive potential (d 0 < | f 0 |/2) or a bound state (d 0 > | f 0 |/2). It was reported that a small variation of the Λ-Λ repulsive core parametrization leads to inverse scattering lengths within −0.27 fm −1 < f −1 0 < 4 fm −1 and effective ranges up to 16 fm [20]. Other calculations are directly constrained to the Nagara event and result in rather small scattering lengths and moderate effective ranges, like the FG ( f −1 0 = 1.3 fm −1 ; d 0 = 6.59 fm) [21] and the HKMYY ( f −1 0 = 1.74 fm −1 ; d 0 = 6.45 fm) [22] models. It is clear that more experimental data are needed to study the problem in a more quantitative and model-independent way.
An alternative method to study hypernuclei is the investigation of momentum correlations of Λ-Λ pairs produced in hadron-hadron collisions via the femtoscopy technique [23]. The STAR collaboration reported a Λ-Λ scattering length and effective range of f −1 0 = −0.91 ± 0.31 +0.07 −0.56 fm −1 and d 0 = 8.52 ± 2.56 +2.09 −0.74 fm, measured in Au-Au collisions at √ s NN = 200 GeV [24]. These values correspond to a repulsive interaction; however, it was shown that the values and the sign of the scattering parameters strongly depend on the treatment of feed-down contributions from weak decays to the measured correlation. A re-analysis of the data outside the STAR collaboration came to different conclusions [20] and resulted in a shallow attractive interaction.
In a pioneering study [25], the Λ-Λ interaction was studied employing the femtoscopy technique in pp collisions at √ s = 7 TeV. This study demonstrated that the data are consistent with either a bound state or an attractive interaction, however, due to the small data sample no quantitative results were obtained. In this letter, these studies are extended by analyzing final-state momentum correlations in pp collisions at √ s = 13 TeV and p-Pb collisions at √ s NN = 5.02 TeV, recorded by ALICE during LHC Run 2. The small system size in pp and p-Pb gives rise to pronounced correlations from strong final-state interactions due to the small relative distance at which particles are produced. Hence, the large data sets enable a high-precision study of the Λ-Λ strong final-state interaction and provide new experimental constraints on the scattering parameters and the existence of a possible bound state.

Data analysis
The analysis presented in this paper is based on the data samples collected by ALICE [26] during the Run 2 of the LHC (2015)(2016)(2017)(2018) in pp collisions at √ s =13 TeV and p-Pb collisions at √ s NN = 5.02 TeV, combined with the previously analyzed Run 1 data from pp collisions at √ s =7 TeV [25]. The event and particle candidate selection criteria follow closely the procedure applied in the Run 1 analysis [25].
The events are triggered using two V0 detectors, which are small-angle plastic scintillator arrays placed on either side of the collision vertex at pseudorapidities 2.8 < η < 5.1 and −3.7 < η < −1.7 [27]. Minimum bias pp and p-Pb events are triggered by the requirement of coincident signals in both V0 detectors, synchronous with the beam crossing time defined by the LHC clock. The V0 detector is also used to reject background events stemming from the interaction of beam particles with the beam pipe materials or beam-gas interactions. Pile-up events with more than one collision per bunch crossing are rejected by evaluating the presence of secondary event vertices [27]. Charged particles are reconstructed by the Inner Tracking System (ITS) [26] and the Time Projection Chamber (TPC) [28], both immersed in a 0.5 T solenoidal magnetic field directed along the beam axis. A uniform detector coverage is assured by requiring the maximal deviation between the reconstructed primary vertex (PV) and the nominal interaction point to be smaller than 10 cm. The PV can be reconstructed with the combined information of the ITS and TPC, and independently with the Silicon Pixel Detector (SPD -one of the three subdetectors of the ITS). If both methods are available, the difference of the z-coordinate between both vertices is required to be smaller than 5 mm. After applying these selection criteria the remaining number of events is 1.0 × 10 9 for the pp at √ s =13 TeV sample and 6.1 × 10 8 for p-Pb at √ s NN = 5.02 TeV. This corresponds to about 90 % and 84 % of all processed events in pp and p-Pb.
The Λ-Λ interaction is the main focus of the present study. As will be explained in the next section, the p-p correlation function is an essential input for the femtoscopic analysis of Λ-Λ. Therefore, the reconstruction of both protons and Λ particles will be described in the following paragraphs. To increase the statistical significance of the result, the anti-particle pairs are measured as well.
The selection of the proton candidates follows the analysis strategy used in pp collisions at √ s =7 TeV [25]. The particle identification (PID) is determined by the number of standard deviations nσ between the hypothesis for a proton and the experimental measurement of the specific energy loss dE/dx in the TPC or the timing information from the Time-Of-Flight (TOF) detector [29]. The analyzed tracks are selected within the kinematic range 0.5 < p T < 4.05 GeV/c and |η| < 0.8. The PID is performed only with the TPC for tracks with p < 0.75 GeV/c, by requiring |nσ | < 3. To maintain the purity of tracks with p > 0.75 GeV/c, the |nσ | is calculated from combining the TPC and TOF information. The contribution of secondary particles, which stem from electromagnetic and weak decays or the detector material, are a contamination in the signal. The fractions of primary and secondary protons are extracted using Monte Carlo (MC) template fits to the distance of closest approach of the particles to the PV [30]. The MC distributions are generated using Pythia 8.2 [31] for the pp and DPMJET 3.0.5 [32] for the p-Pb case, filtered through the ALICE detector and reconstruction algorithm [26]. The proton purity in pp (p-Pb) is found to be 99 (97)% with a primary fraction of 85 (86)%.
The Λ particles are reconstructed via the decay Λ → pπ − , which has a branching ratio of 63.9% and cτ = 7.89 cm [33]. For the reconstruction of the Λ the charge conjugate decay is employed. The interaction rate of the LHC varied during different periods of the pp running. To maintain a constant purity that is independent of the interaction rate, in addition to the selection criteria used for the analysis of pp collisions at √ s =7 TeV [25], the charged decay tracks must either have a hit in one of the SPD or Silicon Strip Detector (SSD -ITS subdetector) layers or a matched TOF signal. After applying all selection criteria the final Λ and Λ candidates are selected in a 4 MeV/c 2 (∼ 3σ ) mass window around the nominal mass [33]. The fractions of primary and secondary Λ particles are extracted similarly as the protons, while the observable for the template fits is the cosine of the opening angle α between the Λ momentum and the vector pointing from the PV to the Λ decay vertex. The Λ purity in pp (p-Pb) is found to be 97 (94)% with a primary fraction of 59 (50)%.

Analysis of the correlation function
The method used to investigate the Λ-Λ interaction relies on particle pair correlations measured as a function of k * , defined as the single-particle momentum in the pair rest frame [23]. The observable of interest C( p 1 , p 2 ) is defined as the ratio of the probability of measuring simultaneously two particles with momenta p 1 and p 2 , to the product of the single-particle probabilities: In the absence of correlations, the numerator factorizes and the correlation function becomes unity. The femtoscopy formalism [23] relates the correlation function for a pair of particles, to their effective twoparticle emitting source function S(r) and the two-particle wave function Ψ( k * , r): where r is the relative distance between the points of emission of the two particles. This definition of C(k * ) assumes that the emission source is not dependent on k * , it is spherically symmetric and the emission of all particles is simultaneous. The EPOS transport model [34] predicts an emission source that does not fully satisfy the above assumptions. However, it was verified that the above simplifications result in very mild deviations in the correlation functions, which are negligible for the present analysis.
For a spherical symmetric potential the angular dependence of the wave-function is trivially integrated out. Thus the direction of k * becomes irrelevant on the left-hand side of Eq. 2. Particles with large relative momenta q * = 2k * are not correlated, leading to C(k * → ∞) = 1.
The strong interaction has a typical range of a few femtometers and thus a significant modification of the wave function with respect to its asymptotic form is expected only for r 2 fm. Consequently, for small emission sources the correlation function will be particularly sensitive to the strong interaction potential. Experimentally, a small emission source can be formed in pp and p-Pb collisions [25,35].
In the current analysis, it is assumed that the emission profile is Gaussian and that the p-p and Λ-Λ systems are characterized by a common source size r 0 = r p-p = r Λ-Λ , which is determined by fitting the p-p correlation function and then used for the investigation of the Λ-Λ interaction.
Two different frameworks are available for the computation of C(k * ). The first tool used in this analysis is the "Correlation Analysis Tool using the Schrödinger equation" (CATS) [35]. Here, a local potential V (r) is used as the input to a numerical evaluation of the wave function and the corresponding correlation function. CATS delivers an exact solution and this tool is used to model the p-p correlation using a Coulomb and an Argonne v 18 potential [36] for the strong interaction. The known p-p interaction allows the source size r 0 to be extracted from the fit to the measured correlation function.
The second tool is the Lednický model [37], which assumes a Gaussian emission source and evaluates the wave function in the effective-range expansion. In this approach, the interaction is parameterized in terms of the scattering length f 0 and the effective range d 0 . This approach produces a very accurate approximation for C(k * ) in case d 0 r 0 , while for smaller values of r 0 the approximate solution may become unstable, in particular for negative values of f 0 [25]. However, it is known that the Lednický model can be used to model the p-Λ correlation function even for a source size of r 0 = 1.2 fm, with a deviation from the exact solution of less than 4% [35]. It is therefore expected that this model can successfully be used to study the Λ-Λ interaction, even in small collision systems. Nevertheless, the validity of the approximation will be further verified in the next section. of the individual components of the p-p, p-Λ, p-Ξ − and Λ-Λ correlation functions. The sub-indexes are used to indicate the mother particle in case of feed-down. Only the non-flat feed-down (residual) contributions are listed individually, while all other contributions are listed as "flat residuals (res.)". All misidentified (fake) pairs are assumed to be uncorrelated, thus resulting in a flat correlation signal. Experimentally, the correlation function is defined as where N same (k * ) and N mixed (k * ) are the same and mixed event distributions, while N is a normalization constant determined by the condition that particle pairs with large relative momenta are not correlated. In small collision systems C exp (k * ) often has a long-range tail due to momentum conservation, and a related approximately linear non-femtoscopic background extending to low k * [25]. The latter is incorporated by including a linear term in the fit function.
To increase the statistical significance of C exp (k * ) the particle-particle (PP) and antiparticle-antiparticle (PP) correlations are combined using their weighted mean C exp (k * ) = N PP C exp,PP (k * ) ⊕ N PP C exp,PP (k * ), with the normalization performed in the range 240 < k * < 340 MeV/c, which is unaffected by femtoscopic correlations. It was verified that N PP C exp,PP (k * ) = N PP C exp,PP (k * ) within the statistical uncertainties.
The systematic uncertainties of the experimental correlation function are evaluated by varying the selection criteria of the proton and Λ candidates within 20%, following the procedure used for the analysis of the pp collisions at √ s =7 TeV [25]. Nevertheless, by performing a Barlow test [38], the systematic uncertainties were found to be insignificant compared to the statistical uncertainties.
Momentum resolution effects modify the correlation function by at most 10% and are accounted for by correcting the theoretical correlation function [25]. The measured experimental correlation function contains not only the correlation signal of interest, but additionally accumulates residual contributions from feed-down particles. These are considered in the theoretical description of the correlation by using the linear decomposition of the total correlation function into where the sum runs over all contributions, the λ parameters are the weight factors for the different contributions to the total correlation and i = 0 corresponds to the primary correlation. The λ coefficients are determined in a data-driven approach by performing Monte Carlo template fits to the data, using Pythia and DPMJET in pp and p-Pb collisions, respectively. The values obtained are summarized in Table 1. The individual contributions C i (k * ) are modeled either using CATS or the Lednický model. The non-genuine (i = 0) contributions include additional kinematic effects which lead to a smearing of their corresponding correlation functions [39]. As the correlation strength of these residuals is strongly damped one can assume that C i =0 (k * ) ≈ 1 [40]. The only significant contribution is p-Λ→p-p, where the p-Λ interaction is modeled using the scattering parameters from a next-to-leading order (NLO) χEFT calculation [41] and the corresponding correlation function is computed using the Lednický model. The remaining residuals are considered flat, apart from p-Ξ − →p-Λ, p-Σ 0 →p-Λ and p-Ξ(1530) − →p-Ξ − , where the interaction can be modeled. For the p-Ξ − interaction a recent lattice QCD potential, from the HAL QCD collaboration [42,43], is used. The p-Σ 0 is modeled as in [44], while p-Ξ(1530) − is evaluated by taking only the Coulomb interaction into account.
After all corrections have been applied to C tot (k * ), the final fit function is obtained by multiplying it with a linear baseline (a + bk * ) describing the normalization and non-femtoscopy background [25] C fit (k * ) = (a + bk * )C tot (k * ).   Two further systematic variations are performed for the p-p correlation. The first concerns the possible effect of non-femtoscopy contributions to the correlation functions, which can be modeled by a linear baseline (see Eq. 5) with the inclusion of b as a free fit parameter. The final systematic variation is to model the p-Λ feed-down contribution by using a leading-order (LO) [41,45] computation to model the interaction. The effect of the latter is negligible, as the transformation to the p-p system smears the differences observed in the pure p-Λ correlation function out.
To investigate the Λ-Λ interaction the source sizes are fixed to the above results and the Λ-Λ correlations from all three data sets are fitted simultaneously in order to extract the scattering parameters. The correlation functions show a slight non-flat behaviour at large k * , especially for the pp collisions at √ s = 13 TeV (right panel in Fig. 1). Thus the fit is performed by allowing a non-zero slope parameter b (see Eq. 5). The fit range is extended to k * < 460 MeV/c in order to better constrain the linear baseline. Due to the small primary λ parameters (see Table 1) the Λ-Λ correlation signal is quite weak and the fit shows a slight systematic enhancement compared to the expected C tot (k * ) due to quantum statistics only, suggestive of an attractive interaction. However, the current statistical uncertainties do not allow the Λ-Λ scattering parameters to be extracted from the fit. Therefore, an alternative approach to study the Λ-Λ interaction will be presented in the next section. Systematic uncertainties related to the Λ-Λ emission source may arise from several different effects, which are discussed in the rest of this section.
Previous studies have revealed that the emission source can be elongated along some of the spatial directions and have a multiplicity or m T dependence [46,47]. In the present analysis it is assumed that the correlation function can be modeled by an effective Gaussian source. The validity of this statement is verified by a simple toy Monte Carlo, in which a data-driven multiplicity dependence is introduced into the source function and the resulting theoretical p-p correlation function computed with CATS. The deviations between this result and a correlation function obtained with an effective Gaussian source profile are negligible.
Possible differences in the effective emitting sources of p-p and Λ-Λ pairs due to the strong decays of broad resonances is evaluated via simulations and estimated to have at most a 5% effect on the effective source size r 0 . This is taken into account by including an additional systematic uncertainty on the r Λ-Λ value extracted from the fit to the p-p correlation.

Results
In order to extract the Λ-Λ scattering parameters, the correlation functions measured in pp collisions at √ s =7, 13 TeV as well as in p-Pb collisions at √ s NN = 5.02 TeV are fitted simultaneously. The right panel in Fig. 1 shows the Λ-Λ correlation function obtained in pp collisions at √ s = 13 TeV together with the result from the fit.   One option is to use a local potential and obtain C(k * ) based on the exact solution from CATS, with the source size fixed to the value obtained from the fit to the p-p correlations. Many of the existing model predictions are summarized in [20] and the corresponding potentials V (r) are parametrized in a local form using a double-Gaussian function. The correlation function depends on the nature of the underlying interaction and Fig. 2 shows the experimental Λ-Λ correlations measured in pp collisions at √ s = 13 TeV (left panel) and p-Pb collisions at √ s NN = 5.02 TeV (right panel) together with the correlation functions obtained for different meson-exchange interaction potentials employing CATS. Models with a strongly attractive interaction ( f −1 0 1 and positive), like the Ehime [17] potential, result in a large enhancement of the correlation function at low momenta which overshoots the data significantly both in pp and p-Pb collisions. The same is valid for potentials corresponding to a shallow bound state ( f −1 0 → 0 and negative), e.g. NF44 [19].
The other tested potentials correspond either to a bound state or a shallow attractive ( f −1 0 1) nonbinding interaction. However, those two very different scenarios result in similar correlations and are difficult to separate. This is evident from Fig. 2 as all of the ESC08 [48], HKMYY [22] and Nijmegen ND46 [18] models produce comparable results and are compatible with the experimental data, even though their scattering parameters are different. In particular, ND46 predicts a bound state, while the ESC08 and HKMYY models describe a shallow attractive potential and the latter is consistent with hypernuclei data [7,8].
The Lednický model can be used to compute C(k * ) for any f −1 0 and d 0 . Thus a scan over the scattering parameters can be preformed and the agreement to the experimental data can be quantified. The Lednický model breaks down for source sizes smaller than the effective range, especially when dealing with repulsive interactions [25], as it produces unphysical negative correlation functions. As there are no realistic models predicting such an interaction, this study is not affected. Nevertheless, all models described in [20] are explicitly tested by comparing the correlation functions obtained using the exact solution provided by CATS with the approximate solution evaluated using the Lednický model. The deviations are on the percent level and are neglected.
Another assumption, which the Lednický model is based on, is a Gaussian profile of the source. The EPOS [34] transport model predicts a non-Gaussian emission profile [35], and the effects of short lived resonances are included. This source was adopted in CATS, by tuning its width such as to describe the p-p correlation function, and the predicted C(k * ) for all of the ND and NF models, shown in Fig. 3, were compared to the Λ-Λ correlation function in pp collisions at √ s = 13 TeV. The deviations in χ 2 compared to the case of a Gaussian source are within the uncertainty, justifying the use of a Gaussian source.
To quantify the uncertainties of f −1 0 and d 0 , and estimate the confidence level of each parameter set, a Monte Carlo method is used. In the current work the approach described in [49] is followed, which is closely related to the Bootstrap method. The strategy is to use the Lednický model to perform a scan over the parameter space spanned by f −1 0 ∈ [−2, 5] fm −1 and d 0 ∈ [0, 18] fm and refit the Λ-Λ correlation using Eq. 5 when fixing the scattering parameters to a specific value ( f −1 0 , d 0 ) i . The corresponding χ 2 i is evaluated by taking all data sets (pp at √ s = 7 and 13 TeV and p-Pb at √ s NN = 5.02 TeV) into account. The different scattering parameters can be compared by finding the lowest (best) χ 2 best and evaluating ∆χ 2 i = χ 2 i − χ 2 best for each parameter set. This observable, and the associated ( f −1 0 , d 0 ) i , can be directly linked to the confidence level [49]. This can be achieved either by assuming normally distributed uncertainties of ( f −1 0 , d 0 ), or invoking a more sophisticated Monte Carlo study, like the Bootstrap method. The latter is used in the current analysis.
The resulting exclusion plot is presented in Fig. 3, where the color code corresponds to the confidence level nσ for a specific choice of scattering parameters. In the computation only the statistical uncer-  tainties are taken into account, as the systematic uncertainties are negligible according to the Barlow criterion [38]. The predicted scattering parameters of all discussed potentials are highlighted with different markers and the phase space region in which the Lednický model produces an unphysical correlation is specified by the black hatched area. In this region the effective range expansion breaks down and the Lednický equation leads to a negative correlation function. While the STAR result [24] is located in this region, all theoretical models exclude the possibility of a repulsive Λ-Λ interaction with large effective range. Moreover a re-analysis of the STAR data [20] demonstrated that a more realistic treatment of the residual correlations leads to an inversion of the sign of the scattering length, that corresponds to an attractive potential. The imposed limit on the scattering length is f −1 0 > 0.8 fm −1 [20]. This result can be tested within the current work, and Fig. 3 demonstrates that the ALICE data can extend those constraints. In particular the region corresponding to a strongly attractive or a very weakly binding short-range interaction (small | f −1 0 | and small d 0 ) is excluded by the data, while a shallow attractive potential (large f −1 0 ) is in very good agreement with the experimental results obtained from this analysis. A Λ-Λ bound state would correspond to negative f −1 0 and small d 0 values. The present data are compatible with such a scenario, but the available phase space is strongly constrained. The HKMYY [22], FG [21] and HAL QCD [50] values are of particular interest, as the first two models are tuned to describe the modern hypernuclei data, while the latter is the latest state-of-the-art lattice computation from the HAL QCD collaboration. The lattice results are preliminary and predict the scattering parameters f −1 0 = 1.45 ± 0.25 fm −1 and d 0 = 5.16 ± 0.82 fm [50]. All three models are compatible with the ALICE data, providing further support for a shallow attractive Λ-Λ interaction potential. A possible bound state is investigated within the effective-range expansion by computing the corresponding binding energy from the relation [51,52] This relation is only valid for bound states, which are characterized by negative f −1 0 values. Further, the binding energy has to be a real number, thus the expression 1 + 2d 0 f −1 0 has to be positive, which implies that at least one of the parameters f −1 0 or d 0 has to be small in absolute value. With these restrictions

Summary
In this Letter, new data on p-p and Λ-Λ correlations in pp collisions at √ s = 13 TeV and p-Pb collisions at √ s NN = 5.02 TeV are presented. Together with the results from a pioneering study on two-baryon correlations in pp at √ s = 7 TeV, these data allow for a detailed study of the Λ-Λ interaction with unprecedented precision.
Each data set is analyzed separately by extracting the p-p and Λ-Λ correlation functions. The former are used to constrain the size of the source r 0 , which is assumed to be the same for p-p and Λ-Λ pairs. The Λ-Λ interaction is then investigated by testing the combined compatibility of all data sets to different model predictions and scattering parameters. The HKMYY and FG models, which are tuned to hypernuclei data, and the lattice calculations performed by the HAL QCD collaboration predict a shallow attractive interaction potential. The ALICE data manifest very good agreement with these predictions. Nevertheless, the data is also compatible with the existence of a bound state, given a binding energy of B ΛΛ = 3.2 +1.6 −2.4 (stat) +1.8 −1.0 (syst) MeV. The Run 3 of the LHC is expected to further increase the statistical significance of the Λ-Λ correlation function and allow the scattering parameters to be constraint even more precisely in the future.

Acknowledgements
The ALICE collaboration is grateful to the HAL QCD collaboration for providing lattice results regarding the Λ-Λ interaction. We are particularly thankful to Prof. Tetsuo Hatsuda and Prof. Kenji Sasaki for the precious suggestions and stimulating discussions.
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