# Boosting quantum Monte Carlo and alleviating sign problem by Gutzwiller
projection

Wei-Xuan Chang  Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
University of Chinese Academy of Sciences, Beijing 100049, China  Zi-Xiang Li
<EMAIL_ADDRESS>Beijing National Laboratory for Condensed Matter Physics
and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
University of Chinese Academy of Sciences, Beijing 100049, China

###### Abstract

Here we develop a new scheme of projective quantum Monte-Carlo (QMC)
simulation combining unbiased zero-temperature (projective) determinant QMC
and variational Monte-Carlo based on Gutzwiller projection wave function,
dubbed as “Gutzwiller projection QMC”. The numerical results demonstrate that
employment of Gutzwiller projection trial wave function with minimum energy
strongly speed up the convergence of computational results, thus tremendously
reducing computational time in the simulation. More remarkably, we present an
example that sign problem is enormously alleviated in the Gutzwiller
projection QMC, especially in the regime where sign problem is severe. Hence,
we believe that Gutzwiller projection QMC paves a new route to improving the
efficiency, and alleviating sign problem in QMC simulation on interacting
fermionic systems.

Introduction: Demystifying quantum many-body physics in strongly correlated
systems is of central importance in modern condensed matter physics.
Developing efficient numerical approaches to solve quantum many-body systems
in more than one dimension is particularly crucial. Among various quantum
many-body numerical algorithms, quantum Monte-Carlo (QMC) plays a vital role
because it is unbiased and approximation-free[1, 2, 3, 4, 5, 6]. However, QMC
suffers from the notorious sign problem[7, 8, 9, 10], which strongly hinders
application of QMC simulation to many strongly correlated systems, for
instance Hubbard model at generic filling. Hence, solving or alleviating sign
problem in quantum many-body models potentially featuring intriguing physics
will definitely lead to a substantial progress in understanding strongly
correlated physics in quantum many-body systems[11, 12, 13, 14, 15, 16, 17,
18, 19, 20]

On the other side, for quantum Monte-Carlo simulation on interacting fermionic
systems, the computational complexity is generally cubic in linear system
size, which strongly hampers the applications to the fermionic systems with
large system sizes[6]. Despite the recent progress on the development of
improved algorithms for fermionic QMC[21, 22, 23], in general the QMC
simulations on interacting fermionic systems are much more costly compared
with the ones on spin or bosonic systems. Hence, although in the simulation
free from sign problem, it is extremely difficult to access the accurate
properties close to the thermodynamic limit in studying some important issues,
for instance quantum criticality[24, 25, 26] and competing ordering[27].
Designing powerful algorithms for interacting fermionic models with high
efficiency is thus immensely desirable.

Figure 1:  The contour plot of energy versus variational parameters in
Gutzwiller projection wave function for $U=3.5$(a) and $U=4.0$(b). The
parameters $(g,M_{\rm{N}})$ with minimum energy are indicated by the white
circles. The optimal parameters are $g=0.55$ and $M_{\rm{N}}=0.07$ for
$U=3.5$, and $g=0.6$ and $M_{\rm{N}}=0.21$ for $U=4.0$.

To this end, we propose a new scheme of QMC simulation dubbed as “Gutzwiller
projection QMC” to speed up the simulation, and more remarkably, alleviate
sign problem in interacting fermionic models. The basic idea of the approach
is the combination of variational Monte-Carlo based on the Gutzwiller
projection variational wave function[28, 29, 30, 31] and intrinsically
unbiased projection QMC. We implement a mean-field wave function under
Gutzwiller projection as the trial wave function, and employ standard
procedure of projection QMC to access the ground-state of an interacting
Hamiltonian without involving any uncontrolled approximations. In the
framework of Hubbard-Strotonovich transformation utilized in the PQMC, the
optimal Gutzwiller variational wave function with minimum energy is achieved
efficiently. Compared with the conventional PQMC simulation which implements a
slater-determinent trial wave function, much smaller projection parameter is
required to guarantee the convergence of results in Gutzwiller projection QMC,
hence immensely reducing the computational time. More crucially, the
systematic calculations in specific models reveal that sign problem is greatly
mitigated in the Gutzwiller projection QMC.

Figure 2: The results of simulation on spinful Hubbard model at half filling
for Gutzwiller QMC and conventional PQMC with slater-determinant trial wave
function. The results of ground-state energy versus projective parameter
$\Theta$ for (a) $U=3.5$, (b) $U=3.75$ and (c) $U=4$. The results of AFM
structure factor $S_{AFM}$ versus projective parameter $\Theta$ for (d)
$U=3.5$, (e) $U=3.75$ and (f) $U=4$. GP, SD(AFM), SD(NI) denote the results of
simulation with Gutzwiller projective wave function, AFM mean-field slater-
determinant wave function and non-interacting slater-determinant wave function
as the trial wave functions, respectively.

Method:  In this section, we briefly illustrate the methods of Gutzwiller QMC,
combining variational Monte-Carlo based on Gutzwiller projection wave function
and unbiased projection QMC. To avoid complexity, we consider the typical wave
function with on-site Gutzwiller projection to illustrate our strategy:
$\left|\psi_{G}\right\rangle=e^{-g\sum_{i}n_{i\uparrow}n_{i\downarrow}}\left|\psi_{M}\right\rangle$,
where $\left|\psi_{M}\right\rangle$ is a slater-determant wave function
featuring mean-field ordering, and $e^{-gn_{i\uparrow}n_{i\downarrow}}$ is
Gutzwiller projection with projective parameter $g$. The first step is
minimizing the expectation value of Hamiltonian under the given Gutzwiller
projection variational wave function, which yields the optimal parameters of
projective parameter $g$ and mean-field order parameter in slater-determinant
wave function $\left|\psi_{M}\right\rangle$. Here we perform H-S
transformation on the Gutzwiller projection term:

$\displaystyle
e^{-gn_{i\uparrow}n_{i\downarrow}}=\frac{1}{2}e^{-g/4}\sum_{s_{i}=\pm
1}e^{\lambda s_{i}(n_{i\uparrow}+n_{i\downarrow})}$ (1)

where $\cosh\lambda=e^{\frac{g}{2}}$ and $s_{i}=\pm 1$ is the auxiliary field
defined on each site $i$. Then the expectation value
$\langle\hat{O}\rangle=\frac{\langle\psi_{G}|\hat{O}|\psi_{G}\rangle}{\langle\psi_{G}|\psi_{G}\rangle}$
is straightforwardly achieved utilizing the standard procedure of determinant
QMC. Notice that this approach of yielding the expectation value of
observables based on Gutzwiller projection wave function is significantly
faster than the conventional approach in variational Monte-Carlo. Then we
search the optimal parameters in Gutzwiller projection wave function by
minimizing the expectation value of energy.

Upon obtaining the optimal Gutzwiller projection wave function, we compute the
expectation value of observable as
$\langle\hat{O}\rangle=\frac{\langle\psi_{T}|e^{-\Theta\hat{H}}\hat{O}e^{-\Theta\hat{H}}|\psi_{T}\rangle}{\langle\psi^{o}_{G}|e^{-2\Theta\hat{H}}|\psi^{o}_{G}\rangle}$
with the trial wave function
$|\psi_{T}\rangle=e^{-g\sum_{i}n_{i\uparrow}n_{i\downarrow}}\left|\psi_{M}\right\rangle$.
We employ the standard procedure of Trotter decomposition and H-S
transformation in PQMC, and decouple the Gutzwiller projection using Eq. (1).
Then the ground-state expectation values of observable are easily accessed in
the scheme of conventional PQMC.

The honeycomb Hubbard model: To demonstrate the efficiency of Gutzwiller QMC,
we first apply the approach to spin-1/2 Hubbard model[32] on honeycomb
lattice:

$\displaystyle H=-t\sum_{\left\langle
ij\right\rangle,\sigma}(c^{\dagger}_{i\sigma}c_{j\sigma}+h.c.)+U\sum_{i}n_{i\uparrow}n_{i\downarrow},$
(2)

where $c^{\dagger}_{i\sigma}$ creates an electron on site $i$ with spin
polarization $\sigma=\uparrow$/$\downarrow$, $t$ is the nearest-neighbor (NN)
hopping, and $U$ is the amplitude of onsite Hubbard repulsion. Hereafter we
set $t=1$ as unit of energy. We focus the study at half filling where sign
problem is circumvent upon choosing appropriate HS transformation channel. In
the non-interacting limit namely $U=0$, the model at half filling features
Dirac fermions with Fermi energy located at Dirac point. With increasing
Hubbard interaction amplitude, a quantum phase transition between Dirac
semimetal and antiferromagnetic Mott insulator occurs
$U\\!=\\!U_{c}\\!\approx\\!3.85$ and the transition belongs to chiral-
Heisenberg universality class[33, 34, 35, 36, 37].

We implement Gutzwiller QMC algorithm to study the ground-state properties of
Eq. (2) at half filling. Because AFM is the dominant ordering in the model, it
is natural to choose AFM mean-field wave function with Gutzwiller projection
as the trial wave function in the simulation:

$\displaystyle\left|\psi_{T}\right\rangle=e^{-g\sum_{i}n_{i\uparrow}n_{i\downarrow}}\left|\psi_{N}\right\rangle$
(3)

where $g$ is the parameter of Gutzwiller projection, and $\psi_{\rm{N}}$ is
mean-field wave function featuring Neel AFM ordering, more explicitly, the
ground-state wave function of the Hamiltonian
$H_{N}=H_{0}+M_{N}\sum_{i}(-1)^{\delta_{i}}(n_{i\uparrow}-n_{i\downarrow})$,
where $H_{0}$ is the non-interacting part of Eq. (2), $M_{\rm{N}}$ is Neel AFM
order parameter, and $\delta_{i}=\pm 1$ if site $i$ belongs to A(B)
sublattice. Employing the procedure introduced in the section of Method, it is
straightforward to access the expectation value of Eq. (2) in terms of wave
function under the choice of $g$ and $M_{\rm{N}}$. Fig. 1 depicts the
expectation value of energy with varying Neel order parameter $M_{\rm{N}}$ and
Gutzwiller projection parameter $g$, for several choices of Hubbard
interaction $U=3.5$ and $U=4.0$ located in the DSM phase and AFM ordered
phase, respectively. We obtain the optimal values of parameters $g$ and
$M_{\rm N}$ by minimizing the expectation value of Eq. (2), and achieve the
Gutzwiller variational wave function utilized in the Gutzwiller QMC
simulation.

We compare the convergence of results against projection parameter $\Theta$
using distinct choices of trial wave function. The results of ground-state
energy and AFM structure factor (The definition are shown in the Supplementary
Materials) are evaluated at varying $\Theta$, as depicted in Fig. 2, which
unambiguously show accurate ground-state energy is achieved at much smaller
value of $\Theta$ in Guzwiller QMC, compared with conventional DQMC simulation
which implements slater-determinant wave function as trial wave function. In
the conventional QMC simulation, two different choices of slater-determinant
trial wave function are employed, including the ground-state wave function of
non-interacting part in Eq. (2) and AFM mean-field wave function in the
absence of Gutzwiller projection. The results of ground-state energy and AFM
structure factor $S_{\rm{AFM}}$ exhibit much slower convergence against
projection parameter $\Theta$ in both cases. More surprisingly, within same
number of Monte-Carlo sampling, the statistic errors of observable are
significantly reduced in the Gutzwiller QMC, compared with the conventional
DQMC with slater-determinant trial wave function. The corresponding results of
statistic error are included in the Supplementary Materials.

Figure 3: The results of observable in spinless t-V model at half filling for
Gutzwiller QMC and conventional PQMC with slater-determinant trial wave
function. The results of CDW structure factors $S_{\rm{CDW}}$ for (a) V=1.6
and (c) V=1.35. The results of ground-state energy $E_{0}$ for (b) V=1.6 and
(d) V=1.35. SD denotes the results of conventional PQMC with employment of
non-interacting slater-determinant trial wave function, and GP denotes the
results of Gutzwiller QMC.  Figure 4: The results of sign problem in spinless
t-V model at half filling for Gutzwiller QMC and conventional PQMC with
slater-determinant trial wave function. (a) The result of sign problem versus
projective parameter $\Theta$ for $L=15$ and $V=1.6$. (b) The results of sign
problem versus projective parameter $\Theta$ for $L=15$ and $V=1.35$. (c) The
results of sign problem versus linear system size $L$ for $V=1.6$. SD denotes
the results of conventional PQMC with the employment of non-interacting
slater-determinant trial wave function, and GP denotes the results of
Gutzwiller QMC.

Repulsive spinless honeycomb model: For Honeycomb Hubbard model, we
unambiguously show that utilizing Gutzwiller trial wave function strongly
expedite the convergence of results against projection parameter $\Theta$,
resulting in a tremendous improvement of the efficiency in the QMC simulation.
In this section, we apply Gutzwiller QMC to the spinless t-V model, namely
spinless fermions model with nearest-neighbor (NN) interaction, on honeycomb
lattice. The Hamiltonian of the model reads:

$\displaystyle H=-t\sum_{\left\langle
ij\right\rangle}(c^{\dagger}_{i}c_{j}+h.c.)+V\sum_{\left\langle
ij\right\rangle}(n_{i}-\frac{1}{2})(n_{j}-\frac{1}{2})$ (4)

where $c_{i}$ is the annihilation operator of fermion on site $i$, $t$ is the
NN hopping amplitude, and $V>0$ denotes the density repulsive interaction
between NN sites. We focus on the half filling of the model. The quantum phase
diagram of the model at half filling has been extensively investigated in
recent years, which features a quantum phase transition from DSM to a charge-
density wave (CDW) insulating phase with increasing NN density interaction
strength[38, 39]. The appearance of sign problem depends on the schemes of HS
transformation. The simulation is sign problematic if the NN density
interaction is decoupled in the hopping channel. Recently, it is shown that
although sign problem exists, the behaviour of average sign exhibit two
distinct behaviours in the weak and strong coupling regimes[40]. In the weak
coupling regime, the model is asymptotic sign-free, namely average sign
asymptotically increases to one as system size increases, whereas in the
strong coupling regime the average sign exhibit exponentially decaying scaling
consistent with conventional scaling behaviour of sign problem in QMC
simulation. Here, we decouple the interaction in the density channel and
perform Gutzwiller QMC simulation on spinless t-V model at half filling,
aiming at investigating the effect of Gutzwiller trial wave function on the
efficiency of simulation, and more importantly, the behaviour of sign problem.

In the simulation of spinless honeycomb t-V model, we choose CDW mean-field
wave function with Gutzwiller projection on the NN bond as the trial wave
function:

$\displaystyle\left|\psi_{T}\right\rangle=e^{-g\sum_{\left\langle
ij\right\rangle}n_{i}n_{j}}\left|\psi_{C}\right\rangle$ (5)

where $g$ is the parameter of Gutzwiller projection on NN bonds. $\psi_{C}$ is
a mean-field wave function with CDW ordering, generated as the ground-state
wave function of the mean-field Hamiltonian $H_{\rm C}=H_{0}+\Delta_{\rm
C}\sum_{i}(-1)^{\delta_{i}}n_{i}$, where $H_{0}$ is the non-interacting part
in Hamiltonian Eq. (4), $\Delta_{\rm C}$ is CDW order parameter and
$\delta_{i}=\pm 1$ if site $i$ belongs to A(B) sublattice. Similar to spinful
Honeycomb Hubbard model, we access the optimal variational parameters $g$ and
$\Delta_{\rm C}$ by minimizing expectation value of $H$ defined in Eq. (4). We
perform QMC simulation with the corresponding Gutzwiller trial wave function.
The results of ground-state energy and CDW structure factors, with the
detailed definitions included in Supplementary Materials, for several values
of $V$ are presented in Fig.3, explicitly demonstrating that smaller value of
$\Theta$ is sufficient to access the accurate ground-state results of energy
and CDW structure factors in Gutzwiller QMC. The improvement of efficiency is
particularly pronounced in the CDW ordered phase ($V>V_{c}\sim 1.35$).

Then we investigate the behaviour of sign problem in the Gutzwiller QMC
simulation. For simulation of Eq. (4) with HS transformation in density
channel, as aforementioned, the model is intrinsically sign-problematic only
as interaction is strong $V>V^{*}\sim 1.2$, whereas the model displays
asymptotic sign-free behaviour in the weak coupling regime, namely average
sign increases and approaches to one as system size increases. Hence, we focus
on the strong coupling regime where sign problem is severe. We calculate the
average sign as a function of projection parameter $\Theta$ for $V=1.35$ and
$V=1.6$, as depicted in Fig. 4(a) and Fig. 4(b), respectively. In comparison,
we present the results employing conventional DQMC with slater-determinant
trial function generated as the ground state of non-interacting part in Eq.
(4). Intriguingly, the results of average sign is obviously increased in the
simulation with Gutzwiller trial wave function, compared with the ones in
conventional DQMC with slater-determinant trial wave function, hence sign
problem is significantly alleviated in Gutzwiller QMC. Furthermore, we plot
the average sign versus linear system size in the simulation with Gutzwiller
QMC and conventional QMC for $V=1.6$ (shown in Fig. 4(c)), which unequivocally
demonstrates that sign problem in spinless honeycomb t-V model at half filling
is tremendously mitigated with the employment of Gutzwiller trial wave
function.

Discussions and concluding remarks: In summary, we develop a new scheme of
zero-temperature (projective) DQMC boosted by Gutzwiller projection. To
demonstrate the remarkable efficiency, we apply the approach to two typical
quantum many-body models. In spinful Honeycomb Hubbard model, the results
clearly show that important observables including ground-state energy and AFM
structure factors exhibit much faster convergence against projection parameter
$\Theta$, thus resulting in tremendous reduction of needed computational time
to achieve certain accuracy. In spinless Honeycomb t-V model, the convergence
of observable is also achieved with smaller projection parameter in Gutzwiller
QMC, compared with the conventional algorithm of projective DQMC. More
appealingly, the sign problem of the model is significantly alleviated in
Gutzwiller QMC simulation if the HS transformation is performed in sign-
problematic channel, especially in the regime where sign problem is severe. To
conclude, Gutzwiller QMC offers a potential route to speeding up the
simulation on the ground-state properties on the interacting fermionic models,
and alleviating sign problem in the simulation on the sign-problematic models.

Acknowledgement: We would like to thank Yi Zhou and Hong Yao for helpful
discussions. This work is supported in part by the start-up grant of IOP-CAS.

## References

  * Scalapino and Sugar [1981] D. J. Scalapino and R. L. Sugar, Phys. Rev. Lett. 46, 519 (1981).
  * Hirsch  _et al._ [1981] J. E. Hirsch, D. J. Scalapino, R. L. Sugar, and R. Blankenbecler, Phys. Rev. Lett. 47, 1628 (1981).
  * White  _et al._ [1989] S. R. White, D. J. Scalapino, R. L. Sugar, and N. E. Bickers, Phys. Rev. Lett. 63, 1523 (1989).
  * Zhang  _et al._ [1995] S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. Lett. 74, 3652 (1995).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Assaad and Evertz [2008] F. F. Assaad and H. G. Evertz, Computational Many-Particle Physics, page 277-356, springer (2008).
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94, 170201 (2005).
  * Li and Yao [2019] Z.-X. Li and H. Yao, Annual Review of Condensed Matter Physics 10, 337 (2019).
  * Loh  _et al._ [1990] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B 41, 9301 (1990).
  * Mondaini  _et al._ [2022] R. Mondaini, S. Tarat, and R. T. Scalettar, Science 375, 418 (2022).
  * Li  _et al._ [2015a] Z.-X. Li, Y.-F. Jiang, and H. Yao, Phys. Rev. B 91, 241117 (2015a).
  * Li  _et al._ [2016] Z.-X. Li, Y.-F. Jiang, and H. Yao, Phys. Rev. Lett. 117, 267002 (2016).
  * Umrigar  _et al._ [2007] C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, Phys. Rev. Lett. 98, 110201 (2007).
  * Hangleiter  _et al._ [2020] D. Hangleiter, I. Roth, D. Nagaj, and J. Eisert, Science Advances 6, eabb8341 (2020).
  * Wan  _et al._ [2020] Z.-Q. Wan, S.-X. Zhang, and H. Yao, arXiv preprint arXiv:2010.01141 (2020).
  * Zhang  _et al._ [2022] X. Zhang, G. Pan, X. Y. Xu, and Z. Y. Meng, Phys. Rev. B 106, 035121 (2022).
  * Ouyang and Xu [2021] Y. Ouyang and X. Y. Xu, Phys. Rev. B 104, L241104 (2021).
  * Wang  _et al._ [2015a] L. Wang, Y.-H. Liu, M. Iazzi, M. Troyer, and G. Harcos, Phys. Rev. Lett. 115, 250601 (2015a).
  * Wei  _et al._ [2016] Z. C. Wei, C. Wu, Y. Li, S. Zhang, and T. Xiang, Phys. Rev. Lett. 116, 250601 (2016).
  * Levy and Clark [2021] R. Levy and B. K. Clark, Phys. Rev. Lett. 126, 216401 (2021).
  * He  _et al._ [2019] Y.-Y. He, H. Shi, and S. Zhang, Phys. Rev. Lett. 123, 136402 (2019).
  * Cohen-Stead  _et al._ [2022] B. Cohen-Stead, O. Bradley, C. Miles, G. Batrouni, R. Scalettar, and K. Barros, Phys. Rev. E 105, 065302 (2022).
  * Wang  _et al._ [2015b] L. Wang, M. Iazzi, P. Corboz, and M. Troyer, Phys. Rev. B 91, 235151 (2015b).
  * Berg  _et al._ [2019] E. Berg, S. Lederer, Y. Schattner, and S. Trebst, Annual Review of Condensed Matter Physics 10, 63 (2019).
  * Berg  _et al._ [2012] E. Berg, M. A. Metlitski, and S. Sachdev, Science 338, 1606 (2012).
  * Xu  _et al._ [2019] X. Y. Xu, Z. H. Liu, G. Pan, Y. Qi, K. Sun, and Z. Y. Meng, Journal of Physics: Condensed Matter 31, 463001 (2019).
  * Fradkin  _et al._ [2015] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Rev. Mod. Phys. 87, 457 (2015).
  * Gutzwiller [1963] M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
  * Rice and Ueda [1986] T. M. Rice and K. Ueda, Phys. Rev. B 34, 6420 (1986).
  * Rokhsar and Kotliar [1991] D. S. Rokhsar and B. G. Kotliar, Phys. Rev. B 44, 10328 (1991).
  * Jin  _et al._ [2021] H.-K. Jin, H.-H. Tu, and Y. Zhou, Phys. Rev. B 104, L020409 (2021).
  * Arovas  _et al._ [2022] D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, Annual Review of Condensed Matter Physics 13, 239 (2022), https://doi.org/10.1146/annurev-conmatphys-031620-102024 .
  * Otsuka  _et al._ [2016] Y. Otsuka, S. Yunoki, and S. Sorella, Phys. Rev. X 6, 011029 (2016).
  * Janssen and Herbut [2014] L. Janssen and I. F. Herbut, Phys. Rev. B 89, 205403 (2014).
  * Ladovrechis  _et al._ [2023] K. Ladovrechis, S. Ray, T. Meng, and L. Janssen, Phys. Rev. B 107, 035151 (2023).
  * Zerf  _et al._ [2017] N. Zerf, L. N. Mihaila, P. Marquard, I. F. Herbut, and M. M. Scherer, Phys. Rev. D 96, 096010 (2017).
  * Parisen Toldin  _et al._ [2015] F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Phys. Rev. B 91, 165108 (2015).
  * Li  _et al._ [2015b] Z.-X. Li, Y.-F. Jiang, and H. Yao, New Journal of Physics 17, 085003 (2015b).
  * Wang  _et al._ [2014] L. Wang, P. Corboz, and M. Troyer, New Journal of Physics 16, 103008 (2014).
  * Li  _et al._ [2022] Z.-X. Li, Z.-Q. Wan, and H. Yao, arXiv preprint arXiv:2211.00663 (2022).


HIN-21-006

HIN-21-006

[cern]The CMS Collaboration

#  and () two-particle femtoscopic correlations in PbPb collisions at
$\sqrtsNN=5.02\TeV$

###### Abstract

Two-particle correlations are presented for , , and strange hadrons as a
function of relative momentum in lead-lead collisions at a nucleon-nucleon
center-of-mass energy of 5.02. The dataset corresponds to an integrated
luminosity of 0.607and was collected using the CMS detector at the CERN LHC.
These correlations are sensitive to quantum statistics and to final-state
interactions between the particles. The source size extracted from the
$\PKzS\PKzS$ correlations is found to decrease from 4 to 1fm in going from
central to peripheral collisions. Strong interaction scattering parameters (,
scattering length and effective range) are determined from the $\PGL\PKzS$ and
$\PGL\PGL$ (including their charge conjugates) correlations using the
Lednický–Lyuboshitz model and are compared to theoretical and other
experimental results.

##  0.1 Introduction

Two-particle correlations in relative momentum, so-called femtoscopic
correlations, arising from relativistic heavy ion collisions provide a
powerful tool for studying both the quark-gluon plasma (QGP) that is created
in the collisions, and the subsequent interactions of the emitted particles
[1]. All two-particle correlations are affected by final-state interaction
(FSI) effects, and correlations of identical particles are also sensitive to
the constraints of quantum statistics (QS). The correlations among the neutral
, , and particles, collectively referred to as $\textrm{V}^{0}$ particles, are
of special interest. First, they can be used to determine the space-time
extent of the QGP. In addition, information can be extracted about the strong-
interaction scattering parameters, , the scattering length and the effective
range, that is impossible to obtain from currently achievable scattering
experiments [2, 3, 4, 5, 6]. Because of their relatively heavy mass and the
absence of a Coulomb interaction, femtoscopy based on particles supplements
the more commonly studied pion and charged kaon pairs [7]. The results from
$\PGL\PGL$ correlation studies can help constrain baryon-baryon and, more
specifically, hyperon-hyperon interaction models that are used, for example,
in modeling the composition of neutron stars [8, 9, 10].

Regarding the scattering parameters, of particular interest is establishing
whether the interaction between two particles allows for the existence of the
H-dibaryon, a bound state with quantum numbers $I=0$, $J^{P}=0^{+}$, $S=-2$.
In 1977, R. L. Jaffe predicted the existence of such a six-quark
($\PQu\PQu\PQd\PQd\PQs\PQs$) state having a mass about 81below the threshold
of twice the mass by considering the strong attraction resulting from color
magnetic interactions [11]. Although a double hypernucleus,
${}^{\;\;\;\;6}_{\PGL\PGL}{\mathrm{He}}$, was subsequently observed in the
NAGARA event from the E313 hybrid emulsion experiment at KEK [12, 13], the
observed $\PGL\PGL$ binding energy was not consistent with the conjectured
H-dibaryon [14]. A study of $\PGL\PGL$ correlations may provide additional
information on whether the baryon-baryon interaction can lead to the formation
of the conjectured H-dibaryon.

Recently, the ALICE Collaboration reported on $\PGL\PK$ correlations in lead-
lead (PbPb) collisions at a center-of-mass energy per nucleon pair of
$\sqrtsNN=2.76\TeV$ [15]. According to their findings, the strong force is
repulsive in $\PGL\PKp$ interactions, yet attractive in $\PGL\PKm$
interactions. For the $\PGL\PKzS$ pairs, the uncertainty of the ALICE results
does not permit a definite conclusion on whether the associated strong
interaction is repulsive or attractive. A more precise measurement of
$\PGL\PKzS$ correlations should improve our understanding of the strong
interaction in baryon-meson systems.

This Letter presents $\PKzS\PKzS$, $\PGL\PKzS$, and $\PGL\PGL$ femtoscopic
correlations as a function of relative momentum in PbPb collisions at
$\sqrtsNN=5.02\TeV$, using data recorded by the CMS experiment during the 2018
LHC run. The $\PKzS\PKzS$ correlations are measured in six centrality
intervals within the 0–60% range, where centrality refers to the percentage of
the total inelastic hadronic nucleus-nucleus cross section [16], and 0%
corresponds to the maximum overlap of the colliding nuclei. The $\PKzS\PKzS$,
$\PGL\PKzS$, and $\PGL\PGL$ correlations are measured in an integrated
centrality range (0–80%), with the $\PGL\PGL$ femtoscopic correlation measured
in PbPb collisions at the LHC for the first time. The source size and strong
interaction parameters are determined using the Lednický–Lyuboshitz (LL) model
[17]. Unless otherwise indicated, all measurements include the charge
conjugate states, so $\PGL\PKzS$ and $\PGL\PGL$ include $\PAGL\PKzS$ and
$\PAGL\PAGL$, respectively. Tabulated results are provided in the HEPData
record for this analysis [18].

##  0.2 Experimental setup and data sample

The central feature of the CMS apparatus is a superconducting solenoid of 6m
internal diameter, providing a magnetic field of 3.8T. Within the solenoid
volume there is a silicon pixel and strip tracker, a lead tungstate crystal
electromagnetic calorimeter, and a brass and scintillator hadron calorimeter,
each composed of a barrel and two endcap sections. The silicon pixel detector
[19] is composed of 1856 silicon pixel modules distributed in four 54 cm long
barrel layers at radii of 2.9–16.0cm plus three pairs of endcap disks covering
radii of 4.5–16.1cm at longitudinal distances of 31–51cm from the origin. The
15 148 silicon strip module are arranged in 10 barrel layers at radii of
20–116cm plus 3 pairs of small and 9 pairs of large endcap disk layers.
Charged particles of pseudorapidity $\abs{\eta}<3$ are reconstructed with the
combined system. For particles with transverse momentum of $1<\pt<10\GeV$, the
track resolutions are typically 1.5% in and 20–75$\mum$ in the transverse
impact parameter [20]. The barrel and endcap detectors are extended to the
forward region with two calorimeters which use steel as the absorber and
quartz fibers as the sensitive material. These hadron forward (HF)
calorimeters are located 11.2m from the interaction region, one on each side,
and provide coverage in the range $3.0<\abs{\eta}<5.2$. These detectors are
segmented into multiple $0.175{\times}0.175$ ($\Delta\eta{\times}\Delta\phi$)
“towers”, where $\phi$ is azimuthal angle in radians. Muons are measured in
gas-ionization detectors embedded in the steel flux-return yoke outside the
solenoid. Events of interest are selected using a two-tiered trigger system.
The first level, composed of custom hardware processors, uses information from
the calorimeters and muon detectors to select events at a rate of around
100kHz within a fixed latency of about 4 [21]. The second level, known as the
high-level trigger, consists of a farm of processors running a version of the
full event reconstruction software optimized for fast processing, and reduces
the event rate to around 1kHz before data storage [22]. A more detailed
description of the CMS detector, together with a definition of the coordinate
system used and the relevant kinematic variables, can be found in Ref. [23].

With an integrated luminosity of 0.607 [24, 25], this analysis uses
$4.27\times 10^{9}$ minimum bias events that are triggered by requiring
signals above the readout threshold of 3in each of the HF calorimeters [22].
Background events due to beam-gas interactions and non-hadronic collisions are
filtered offline by applying the procedure described in Ref. [26]. The events
used in this analysis are required to have at least one primary interaction
vertex determined using two or more tracks [27] within a distance of 15cm from
the center of the nominal interaction point along the beam axis and to have at
least two calorimeter towers in each HF detector with energy deposits of more
than 4per tower. The shapes of the clusters in the pixel detector are required
to be compatible with those expected in PbPb collisions in order to suppress
the contamination from events with multiple collisions [28]. The combined
trigger and offline selection efficiency for inelastic events is greater than
95%. The event centrality is obtained from the transverse energy deposited in
both HF calorimeters, using the methodology described in Ref. [29]. The
analysis makes use of a minimum bias Monte Carlo PbPb sample, based on the 1.9
[30] event generator with a full detector simulation using [31].

##  0.3 Reconstruction of and candidates

The and candidates, denoted as $\textrm{V}^{0}$ candidates, used in this study
are reconstructed as in previous CMS analyses [32, 33, 34]. The
$\textrm{V}^{0}$ candidates are found by combining oppositely charged tracks
that pass criteria based on the “loose” selection discussed in Ref. [27]. The
charged tracks are assumed to be $\PGpp\PGpm$ in reconstruction and $\PGpm\Pp$
in reconstruction. For the latter, the higher momentum track is assumed to be
a proton since the proton carries nearly all of the momentum in the decay.
Each of the oppositely charged tracks must have hits in at least three layers
of the silicon tracker, and both tracks must have transverse and longitudinal
impact parameter significances (defined as the parameter value divided by its
uncertainty) with respect to the primary vertex greater than 1. The two tracks
are fitted to a common vertex and the $\chi^{2}$ per degree of freedom (dof)
from the fit must be less than 7. The distance of closest approach between the
two tracks is required to be less than 1cm. As a consequence of the long
lifetime of and particles, the significance of the $\textrm{V}^{0}$ decay
length, which is the three-dimensional distance between the primary and
$\textrm{V}^{0}$ vertices divided by its uncertainty, is required to be
greater than 2.5 to reduce combinatorial background contributions. To remove
candidates misidentified as particles and vice versa, the () candidates must
have a corresponding $\PGp\PGp\;(\Pp\PGp)$ mass more than 14 (7)(corresponding
to approximately 3 times the average resolution) away from the world-average
value [35] of the () mass. The angle $\theta$ between the $\textrm{V}^{0}$
momentum vector and the vector connecting the primary and $\textrm{V}^{0}$
vertices is required to satisfy $\cos\theta>0.999$. This reduces the
contribution from nuclear interactions, random combinations of tracks, and
particles originating from weak decays of and particles.

Further selection of $\textrm{V}^{0}$ candidates is performed with a boosted
decision tree (BDT) [36]. The selection is optimized separately for and
candidates. The discriminating variables include: the collision centrality,
the $\textrm{V}^{0}$ candidate and rapidity ($y$), the distance of closest
approach of the track pair, the three-dimensional decay length and
significance, $\cos\theta$, and the $\textrm{V}^{0}$ vertex fit $\chi^{2}$.
The included variables related to the $\textrm{V}^{0}$ daughters are ,
uncertainty in , $\eta$, the number of hits in the silicon tracker, the number
of pixel detector layers with hits, and the transverse and longitudinal impact
parameter significances with respect to the primary vertex. The BDT training
is performed using the simulated minimum bias sample separated into the signal
and background subsamples using the generator-level information. The mesons
are selected with $1<\pt<8.5\GeV$ and $\abs{y}<1$, while the baryons are
required to have $1.8<\pt<8.5\GeV$ and $\abs{y}<1$. The minimum and maximum
$y$ requirements are used to reduce background while the maximum requirement
is to reduce contributions from jets. The combined $\textrm{V}^{0}$
reconstruction and selection efficiencies are strongly dependent on the
centrality of the event and the of the $\textrm{V}^{0}$. Integrating over the
selected ranges, the combined efficiencies from the most central to peripheral
PbPb collisions are 1–3% for and 1–2% for . The $\textrm{V}^{0}$
reconstruction algorithm does not prevent a track from being used for more
than one $\textrm{V}^{0}$. While this is normally an infrequent occurrence,
selecting pairs of $\textrm{V}^{0}$ particles close together in phase space
makes it a significant contribution. To resolve this problem, for each
correlation measurement, a check of each pair of $\textrm{V}^{0}$ candidates
is performed and if two $\textrm{V}^{0}$ candidates are found to share one or
both daughter tracks, one of the $\textrm{V}^{0}$ candidates is randomly
selected to be removed from the event.

Fits to the invariant mass spectrum are performed using a sum of three
Gaussian functions with a common mean to describe the signal distribution and
a fourth-order polynomial to describe the background. These empirical
functions were chosen to provide a good description of the data. Peak and
sideband invariant mass regions are defined to select events dominated by
signal and background, respectively. Defining $\sigma$ as the average
resolution based on the Gaussian sum, the peak regions are selected to be
within ${\pm}2\sigma$ from the nominal $\textrm{V}^{0}$ mass and are given by
$486<\mathrm{M(\PGpp\PGpm)}<509\MeV$ and
$1111.5<\mathrm{M(\Pp\PGpm)}<1120.4\MeV$ for and candidates, respectively. The
sideband regions are selected to be more than $4\sigma$ from the nominal
$\textrm{V}^{0}$ mass and are given by
$23.5<\abs{\mathrm{M(\PGpp\PGpm)}-497.5}<62.5\MeV$ for candidates and
$1080<\mathrm{M(\Pp\PGpm)}<1107.5\MeV$ together with
$1124.2<\mathrm{M(\Pp\PGpm)}<1160\MeV$ for candidates. Examples of invariant
mass distributions for $\PGpp\PGpm$ and $\Pp\PGpm$ pairs, and their
corresponding fits in the 0–80% centrality range, are shown in Fig. 1.

Figure 1: The invariant mass of (left) and (right), and their corresponding
fits in the 0–80% centrality range. The circles are the data, and the fit is
shown with a solid (red) line for the total fit, and a dashed (green) line for
the background fit. The vertical dashed-dotted (pink) lines indicate the peak
region and the vertical dashed (blue) lines indicate the sideband regions.

##  0.4 Analysis method

The two-particle correlation is constructed as

$C^{\text{obs}}(q_{\text{inv}})=\mathcal{N}\frac{A^{\text{obs}}(q_{\text{inv}})}{B^{\text{obs}}(q_{\text{inv}})},$
(1)

where $C^{\text{obs}}(q_{\text{inv}})$ is the observed normalized pair yield,
corrected for detector effects, as a function of the invariant relative
momentum $q_{\text{inv}}$, defined as [1]

$\begin{split}&q_{\text{inv}}=\sqrt{{-Q^{\mu}Q_{\mu}}},\\\
&Q^{\mu}=(k_{1}-k_{2})^{\mu}-\frac{(k_{1}-k_{2})^{\mu}P_{\mu}}{P^{\mu}P_{\mu}}P^{\mu},\end{split}$
(2)

where $P=k_{1}+k_{2}$, and $k_{1}$ and $k_{2}$ are the four momenta of the
$\textrm{V}^{0}$ particles. Note that for two particles of the same mass, the
second term of $Q^{\mu}$ is zero.

The distribution $A^{\text{obs}}(q_{\text{inv}})$ is the signal distribution
that contains femtoscopic correlations formed by pairing the selected
$\textrm{V}^{0}$ particles from a given event. The reference distribution
$B^{\text{obs}}(q_{\text{inv}})$ is used to correct for phase space effects,
largely removing artifacts due to detector non-uniformities in the
$A^{\text{obs}}(q_{\text{inv}})$ distribution. The
$B^{\text{obs}}(q_{\text{inv}})$ distribution is constructed by mixing the
$\textrm{V}^{0}$ particles from different events [37]. In this procedure, the
$\textrm{V}^{0}$ particle from one event is paired with $\textrm{V}^{0}$
particles from 30 different events. To ensure that the 30 events used in the
mixing are similar to the signal events, the centrality and primary vertex of
each mixed event must be within 5% and 2cm, respectively, of those in the
corresponding signal event. The normalization factor $\mathcal{N}$ is the
ratio of the number of pairs in the reference distribution to that in the
signal distribution. Because of the background in the peak region of the
invariant mass distributions, the measured signal distribution
$(A^{\text{obs}}(q_{\text{inv}}))$ contains contributions from signal-signal
$(A^{\text{ss}}(q_{\text{inv}}))$, signal-background
$(A^{\text{sb}}(q_{\text{inv}}))$, and background-background
$(A^{\text{bb}}(q_{\text{inv}}))$ correlations. The measured
$A^{\text{obs}}(q_{\text{inv}})$ distribution can be written as

$A^{\text{obs}}(q_{\text{inv}})=f^{\text{ss}}A^{\text{ss}}(q_{\text{inv}})+f^{\text{sb}}A^{\text{sb}}(q_{\text{inv}})+f^{\text{bb}}A^{\text{bb}}(q_{\text{inv}}).$
(3)

The distributions $A^{\text{sb}}(q_{\text{inv}})$ and
$A^{\text{bb}}(q_{\text{inv}})$ are obtained from the peak-sideband and
sideband-sideband combinations, respectively. The small amount of background
(signal) contamination in the signal (sideband) region has a negligible effect
on the shape of $A^{\text{sb}}(q_{\text{inv}})$
$(A^{\text{bb}}(q_{\text{inv}}))$. All distributions,
$A^{\text{obs}}(q_{\text{inv}})$, $A^{\text{sb}}(q_{\text{inv}})$, and
$A^{\text{bb}}(q_{\text{inv}})$ are normalized to unity. The parameters,
$f^{\text{ss}}$, $f^{\text{sb}}$, and $f^{\text{bb}}$ are the signal-signal,
signal-background, and background-background fractions, extracted using an
invariant mass fit based on combinatorial analyses with

$\begin{split}&f^{\text{ss}}=\frac{\binom{\mathrm{s}}{2}}{\binom{\mathrm{s}+\mathrm{b}}{2}},\\\
&f^{\text{bb}}=\frac{\binom{\mathrm{b}}{2}}{\binom{\mathrm{s}+\mathrm{b}}{2}},\;\text{and}\\\
&f^{\text{sb}}=1-f^{\mathrm{ss}}-f^{\mathrm{bb}},\end{split}$ (4)

where $\binom{n}{2}=\frac{n(n-1)}{2}$ is the binomial coefficient, which
returns the number of ways that a pair can be chosen from $n$ objects. The
quantities $\mathrm{s}$ and $\mathrm{b}$ in the binomial coefficients are the
number of signal and background particles, respectively, obtained by
integrating the appropriate function from the fit to the invariant mass
distribution.

Once we have all the distributions $(A^{\text{obs}}(q_{\text{inv}})$,
$A^{\text{sb}}(q_{\text{inv}})$, and $A^{\text{bb}}(q_{\text{inv}}))$ and the
parameters $(f^{\text{ss}}$, $f^{\text{sb}}$, and $f^{\text{bb}})$, the
$A^{\text{ss}}(q_{\text{inv}})$ distribution can be extracted using Eq. (3),
with

$A^{\text{ss}}(q_{\text{inv}})=\Big{[}A^{\text{obs}}(q_{\text{inv}})-f^{\text{sb}}(A^{\text{sb}}(q_{\text{inv}}))-f^{\text{bb}}(A^{\text{bb}}(q_{\text{inv}}))\Big{]}/f^{\text{ss}}.$
(5)

The same procedure is followed for the reference distribution. After
extracting the $A^{\text{ss}}(q_{\text{inv}})$ and
$B^{\text{ss}}(q_{\text{inv}})$ distributions, the correlation distribution is
calculated as

$C^{\text{ss}}(q_{\text{inv}})=\mathcal{N}\frac{A^{\text{ss}}(q_{\text{inv}})}{B^{\text{ss}}(q_{\text{inv}})}.$
(6)

While the $C^{\text{ss}}(q_{\text{inv}})$ distribution is corrected for
detector effects and non-$\textrm{V}^{0}$ backgrounds, it still includes non-
femtoscopic background correlations, such as those associated with elliptic
flow [38], minijets [7], resonance decays [7], and energy-momentum
conservation [39]. The non-femtoscopic background contribution is modeled
using an empirically determined double Gaussian function

$\Omega(q_{\text{inv}})=N\left(1+\alpha_{1}\re^{-q_{\text{inv}}^{2}R^{2}_{1}}\right)\left(1-\alpha_{2}\re^{-q_{\text{inv}}^{2}R^{2}_{2}}\right),$
(7)

where $N$, $\alpha_{1}$, $\alpha_{2}$, $R_{1}$, and $R_{2}$ are fit
parameters. This function was selected for its reproduction of the
distributions in both real data at high $q_{\text{inv}}$ and simulated data
that do not include the correlations being measured.

Fits are performed to the $C^{\text{ss}}(q_{\text{inv}})$ distributions to
extract the source size and strong interaction scattering parameters. As the
$\textrm{V}^{0}$ particles are neutral, the Coulomb interaction is absent.
However, the correlations are sensitive to QS and FSI effects, with $s$-wave
interactions assumed to dominate for the small relative momenta of the
particle pairs analyzed. The correlation distribution for all pairs
($\PKzS\PKzS$, $\PGL\PKzS$, and $\PGL\PGL$) is interpreted in the LL model.
This model relates the two-particle correlation function to the source size
and also takes into account FSI effects [17]. The general correlation function
is

$C_{\text{total}}(q_{\text{inv}})=\Big{[}1+\Lam\Big{(}C_{\text{QS}}(q_{\text{inv}})+C_{\text{FSI}}(q_{\text{inv}})\Big{)}\Big{]}\Omega(q_{\text{inv}}),$
(8)

where $C_{\text{QS}}(q_{\text{inv}})$ is the QS function and
$C_{\text{FSI}}(q_{\text{inv}})$ is the FSI function. The parameter is
referred to as the incoherence parameter. In the absence of FSI effects,
equals unity for a perfectly incoherent Gaussian source. Effects such as
resonance decay violate the incoherent source assumption and can lead to
deviations of the parameter from the unity. Its value can also be affected by
non-Gaussian components of the correlations function and by the FSI between
particles.

Neglecting CP violation, the $\PKzS\PKzS$ system can be written as

$\mathinner{|{\PKzS\PKzS}\rangle}=\frac{1}{2}\Bigl{(}\mathinner{|{\PKz\PKz}\rangle}+\mathinner{|{\PKz\PAKz}\rangle}+\mathinner{|{\PAKz\PKz}\rangle}+\mathinner{|{\PAKz\PAKz}\rangle}\Bigr{)}.$
(9)

It can be shown [17, 40] that the resulting correlations follow Bose–Einstein
quantum statistics, with

$C_{\text{QS}}(q_{\text{inv}})=\re^{(-q_{\text{inv}}^{2}R_{\text{inv}}^{2})},$
(10)

where the source radius $R_{\text{inv}}$ reflects the size of the region over
which the particles are emitted.

The FSI for the $\PKzS\PKzS$ correlations is modeled by [17, 40]

$C_{\text{FSI}}(q_{\text{inv}})=\frac{1}{2}\Biggl{[}\left|\frac{f(k)}{R_{\text{inv}}}\right|^{2}+\frac{4\Re
f(k)}{\sqrt{\pi}R_{\text{inv}}}F_{1}(q_{\text{inv}}R_{\text{inv}})-\frac{2\Im
f(k)}{R_{\text{inv}}}F_{2}(q_{\text{inv}}R_{\text{inv}})\Biggr{]},$ (11)

where

$\begin{split}&k=q_{\text{inv}}/2,\\\
&F_{1}(z)=\frac{1}{2}\re^{-z^{2}}\int_{0}^{z}e^{x^{2}}\rd{x},\;\text{and}\\\
&F_{2}(z)=\frac{1-\re^{-z^{2}}}{z}.\end{split}$ (12)

The function $f(k)$ is the $\PKz\PAKz$ $s$-wave scattering amplitude, with
real and imaginary parts $\Re f(k)$ and $\Im f(k)$, respectively. This
amplitude is dominated by the near-threshold $s$-wave isoscalar resonance
$\PfDzP{980}$ and the $s$-wave isovector resonance , with the total scattering
amplitude given by an average of these contributions:
$f(k)=(f_{\PfDzP{980}}(k)+f_{\Paz}(k))/2$. The individual resonance amplitudes
depend on the resonance mass $m_{r}$, with $r=\PfDzP{980}$ or , the kaon mass
$m_{\PK}$, and the resonance couplings $\gamma_{r}$ ($\gamma^{\prime}_{r}$) to
the $\PKz\PAKz$ $(\PGp\PGp$ for $\PfDzP{980}$ and $\PGpz\PGh$ for $\Paz$
channels. Then,
$f_{r}(k)=\gamma_{r}/\\!\left[m^{2}_{r}-\zeta-i\gamma_{r}k-i\gamma^{\prime}_{r}k^{\prime}_{r}\right]$,
where $\zeta=4(m^{2}_{\PK}+k^{2})$ and $k^{\prime}_{r}$ denotes the momentum
in the second ($\PGp\PGp$ or $\PGpz\PGh$) decay channel with the corresponding
partial width $\Gamma^{\prime}=\gamma^{\prime}_{r}k^{\prime}_{r}/m_{r}$ (more
details can be found in Ref. [40]). The scattering amplitude is calculated
using the resonance mass and the coupling parameters from Refs. [41, 42, 43,
43, 44], taken from row C of Table 1 of Ref. [40].

For the correlations involving baryons, the $C_{\text{QS}}(q_{\text{inv}})$
and $C_{\text{FSI}}(q_{\text{inv}})$ functions are [17]

$\begin{split}&C_{\text{QS}}(q_{\text{inv}})=\alpha\re^{(-q_{\text{inv}}^{2}R_{\text{inv}}^{2})},\quad\text{and}\\\
&C_{\text{FSI}}(q_{\text{inv}})=(1+\alpha)\Bigg{[}\frac{1}{2}\frac{|f(k)|^{2}}{R_{\text{inv}}^{2}}\big{(}1-\frac{1}{2\sqrt{\pi}}\frac{d_{0}}{R_{\text{inv}}}\big{)}+\frac{2\Re{f(k)}}{\sqrt{\pi}R_{\text{inv}}}F_{1}(q_{\text{inv}}R_{\text{inv}})-\frac{\Im{f(k)}}{R_{\text{inv}}}F_{2}(q_{\text{inv}}R_{\text{inv}})\Bigg{]},\end{split}$
(13)

where $\alpha={-}1/2$ for $\PGL\PGL$ correlations for two identical fermions
and $\alpha=0$ for $\PGL\PKzS$ correlations as there are no QS effects for
non-identical particles [17]. The scattering amplitude $f(k)$ is parameterized
by a complex scattering length ($f_{0}$) and an effective range ($d_{0}$) with
$f(k)=[{1}/{f_{0}}+d_{0}k^{2}/2-ik]^{-1}$ [17]. The imaginary part of $f_{0}$
is responsible for inelastic processes (annihilation). For an attractive
interaction that is not strong enough to produce a bound state, the real part
of $f_{0}$ is positive, while a repulsive interaction corresponds to a
negative $\Re f_{0}$ of the order of the range of the repulsive potential. In
the presence of a bound state, $\Re f_{0}$ is also negative, but with a much
larger magnitude. The femtoscopic sign convention and notation for the
scattering length differ from those used in nuclear physics, where the
corresponding scattering length $a_{0}=-f_{0}$. As the $\PGL\PKzS$ and
$\PGL\PGL$ correlations each have only one spin state that contributes to the
$s$-wave scattering, Eq. (13) suffices to describe the FSI effects.

Fits to the correlation distribution of all the pairs were performed using Eq.
(8) with the non-femtoscopic background parameters ($N$, $\alpha_{1}$,
$\alpha_{2}$, $R_{1}$, and $R_{2}$) treated as free parameters. For
$\PKzS\PKzS$ correlations, the parameters of interest are $R_{\text{inv}}$ and
, with the scattering amplitude based on previous measurements [41, 42, 43,
43, 44]. The $\PGL\PKzS$ and $\PGL\PGL$ correlations include additional
parameters: $d_{0}$, $\Re f_{0}$, and $\Im f_{0}$. The $\Im f_{0}$ term for
$\PGL\PGL$ correlations is set to zero since there are no baryon-baryon
annihilation processes.

Histograms of the correlation distributions are generated in the range
$0<q_{\text{inv}}<6\GeV$ with 0.02wide bins for the $\PKzS\PKzS$ and
$\PGL\PKzS$ correlations and 0.04wide bins for the $\PGL\PGL$ correlations.
The fits exclude the first $q_{\text{inv}}$ bin to avoid a potential bias from
the method used to address cases where $\textrm{V}^{0}$ candidates share
daughter tracks. Studies using simulated events indicate that only this first
bin is affected by this remediation. Least-square fits are performed to the
experimental data with the uncertainties in the fit parameters calculated
using the minos technique [45]. Examples of correlation measurements and their
fits and corresponding $\chi^{2}/\mathrm{dof}$ values, are presented in Figs.
2 and 3. The $\PKzS\PKzS$ correlations, shown in Fig. 2, are independently
fitted for each of the six centrality bins with $0<\kt<2.5\GeV$, where
$\kt\equiv\abs{\vec{p}_{\mathrm{T},1}+\vec{p}_{\mathrm{T},2}}/2$ is the
average transverse momentum of the particle pair. While the LL model assumes a
Gaussian source function, results from charged-particle correlations have
demonstrated that this assumption breaks down for peripheral collisions [46].
This is likely the cause of the poor fit at low $q_{\text{inv}}$ for
centralities above 40%. The $\PGL\PKzS$ (left) and $\PGL\PGL$ (right)
correlations, shown in Fig. 3, involve fewer events and, therefore, only a
single fit is performed for each, with the data integrated over the centrality
range 0–80% and with no restriction on .

Figure 2: The correlation distributions and fits for $\PKzS\PKzS$ pairs in
different centrality ranges, starting from 0–10% centrality to 50–60%
centrality, with $0<\kt<2.5\GeV$. In each plot, the red circles are the data,
the blue solid line is the fit using Eq. (8), and the green dotted line is the
non-femtoscopic background from Eq. (7). The $\chi^{2}$ and $\mathrm{dof}$
values are for the full $q_{\text{inv}}$ range. The insert plots show the data
and the fit for the $q_{\text{inv}}<0.4\GeV$ region, with the $\chi^{2}$ and
number of bins evaluated in that region.

Figure 3: The correlation distributions and fits for $\PGL\PKzS$ (left) and
$\PGL\PGL$ (right) pairs with 0–80% centrality and no restriction on . In each
plot, the red circles are the data, the blue solid line is the fit using Eq.
(8), and the green dotted line is the non-femtoscopic background from Eq. (7).
The $\chi^{2}$ and $\mathrm{dof}$ values are for the full $q_{\text{inv}}$
range. The insert plots show the data and the fit for the
$q_{\text{inv}}<0.4\GeV$ region, with the $\chi^{2}$ and number of bins
evaluated in that region.

##  0.5 Systematic uncertainties

The systematic uncertainties for the fit parameters are based on the changes
found in the parameter values after individually varying each of the analysis
criteria, as discussed below. In cases with more than one variation for a
single source, the maximum deviation from the nominal value is used. The total
systematic uncertainty is obtained by adding the uncertainties from each
source in quadrature. The BDT discriminant is varied so as to adjust the
signal-to-background ratio, with the signal yield changing by $\pm$15% in the
process. The nominal method to account for $\textrm{V}^{0}$ candidates sharing
daughter tracks is to remove one of the $\textrm{V}^{0}$ candidates at random,
which is then not used by any pair. Two alternative approaches are used, one
in which both $\textrm{V}^{0}$ candidates are removed and another in which,
for events with multiple $\textrm{V}^{0}$ candidate pair combinations, only
the pairs in which the two particles share a daughter are removed. The
systematic uncertainties related to $\textrm{V}^{0}$ signal and background
modeling are investigated by varying the background shape from a fourth- to a
third-order polynomial and the signal shape from a sum of three Gaussian
functions to a sum of two or four Gaussian functions. An alternative non-
femtoscopic background function
$\Omega(q_{\text{inv}})=N(1+B\re^{-|q_{\text{inv}}/\sigma|^{2}})(1+\epsilon
q_{\text{inv}})$ is used to assess the uncertainty associated with the choice
of the non-femtoscopic background function. The selection requirements used to
construct the mixed event sample are varied to require centrality matching of
3 and 7% instead of the nominal 5% and the primary vertex position matching
with 1 and 3cm instead of the nominal 2cm. The effect of the centrality
resolution has been checked and found to be negligible. The peak region
requirement is changed from ${<}2.0\sigma$ to ${<}1.5\sigma$ and
${<}2.5\sigma$ and the sideband region selection from ${>}4.0\sigma$ to
${>}3.5\sigma$ and ${>}4.5\sigma$. The upper limit of the $q_{\text{inv}}$ fit
ranges is changed by $\pm$1and the lower limit is changed to include the first
bin. At low , the tracking efficiency is strongly dependent on . Therefore,
the simulated sample is used to explore possible effects of the tracking
efficiency. Based on these studies, it is found that the $\textrm{V}^{0}$
reconstruction efficiency for the detection of two $\textrm{V}^{0}$ particles
is well described by taking the product of the efficiency for each
$\textrm{V}^{0}$. It is also found that the fit results are only weakly
affected by the $\textrm{V}^{0}$ reconstruction efficiency. This is understood
as a consequence of the signal and reference samples being similarly affected
by the $\textrm{V}^{0}$ efficiency. Differences in the Monte Carlo and
experimental spectra could influence the cancellation of efficiency-dependent
effects in the signal and background correlations. Therefore, a systematic
uncertainty for the efficiency is assessed by comparing the results for the
default simulated sample to one in which the $\textrm{V}^{0}$ distribution is
reweighted to match the data. For the $\PKzS\PKzS$ correlations, an additional
systematic uncertainty is found from varying the mass and coupling parameters
for the $\PfDzP{980}$ and resonances by using rows A, B, and D of Table 1 of
Ref. [40]. The systematic uncertainties are summarized in Table 0.5.

Summary of absolute systematic uncertainties in $\PKzS\PKzS$, $\PGL\PKzS$ and
$\PGL\PGL$ correlation measurements. The values for $R_{\text{inv}}$, $d_{0}$,
$\Re f_{0}$, and $\Im f_{0}$ are infm.  Uncertainty source $\PKzS\PKzS$
$\PGL\PKzS$ $\PGL\PGL$ $R_{\text{inv}}$ $R_{\text{inv}}$ $d_{0}$ $\Re f_{0}$
$\Im f_{0}$ $R_{\text{inv}}$ $d_{0}$ $\Re f_{0}$ BDT cut 0.04–0.18 0.01–0.04
0.19 0.75 0.10 0.07 0.03 0.06 0.43 0.05 0.31 Duplicate $\textrm{V}^{0}$
removal 0.06–0.40 0.01–0.08 0.35 0.92 0.10 0.19 0.11 0.01 1.14 0.05 0.14 Mass
fit function 0.00–0.01 0.00–0.01 0.09 0.05 0.01 0.03 0.03 0.02 0.04 0.01 0.02
Non-femtoscopic func. 0.02–0.16 0.01–0.12 0.02 0.17 0.05 0.07 0.03 0.02 1.02
0.14 0.93 Reference sample 0.03–0.08 0.01–0.05 0.22 0.48 0.12 0.12 0.03 0.10
1.12 0.20 0.76 Peak region 0.00–0.07 0.01–0.02 0.43 0.10 0.05 0.17 0.08 0.22
1.21 0.08 0.35 Sideband region 0.00–0.03 0.00–0.01 0.02 0.02 0.01 0.01 0.00
0.01 0.03 0.01 0.04 Fitting range 0.01–0.11 0.01–0.04 0.20 0.18 0.03 0.08 0.04
0.04 1.79 0.20 0.60 Efficiency 0.03–0.03 0.01–0.01 0.08 0.29 0.06 0.09 0.03
0.02 0.05 0.03 0.04 $\PfDzP{980}$/param. 0.07–0.39 0.03–0.05 Total uncertainty
0.29–0.47 0.08–0.16 0.69 1.34 0.21 0.32 0.16 0.25 2.91 0.34 1.43

##  0.6 Results

The size of the particle emitting source $R_{\text{inv}}$ and the parameter
extracted from the $\PKzS\PKzS$ correlations for $0<\kt<2.5\GeV$ are shown as
a function of centrality in Fig. 4. It is observed that the $R_{\text{inv}}$
value decreases from central to peripheral events, as expected from a simple
geometric picture of the collisions. Over the full centrality range of 0–80%,
$R_{\text{inv}}=3.30\pm 0.10\stat\pm 0.37\syst$fm. The transverse mass can be
calculated as $\mT=\sqrt{\smash[b]{(m_{\text{inv}}/2)^{2}+\kt^{2}}}$, where
$m_{\text{inv}}$ is the invariant mass of the two-particle system [15]. The
average $\left<\mT\right>$ is evaluated from the transverse mass distribution
using two-particle pairs with $q_{\text{inv}}<0.4\GeV$, accounting for
background using the binomial analysis as done for the $q_{\text{inv}}$
distributions. Our results for $R_{\text{inv}}$ agree with the ALICE
$\PKzS\PKzS$ results from PbPb collisions at $\sqrtsNN=2.76\TeV$ at a similar
value [47]. The parameter is seen to decrease from about 0.45 to 0.25 as the
collisions become more peripheral. This decrease could arise from a relative
increase in the contribution from resonance decays or a source function that
becomes increasingly non-Gaussian as the collisions become more peripheral.
The assumption of a Gaussian source function in the LL model may also be
responsible for the relatively poor fits at low $q_{\text{inv}}$ for the most
peripheral collisions, as seen in Fig. 2.

Figure 4: The $R_{\text{inv}}$ (left) and parameter (right) as a function of
centrality. For each data point, the line and shaded area indicate the
statistical and systematic uncertainty, respectively.

Table 0.6 includes the extracted $R_{\text{inv}}$ and parameters as well as
$\left<\mT\right>$ for $\PKzS\PKzS$, $\PGL\PKzS$, and $\PGL\PGL$ combinations
in the 0–80% centrality range. A significant decrease is seen in
$R_{\text{inv}}$ as the $\left<\mT\right>$ increases. Qualitatively similar
results have been found, both for a given pair type in bins of and when
comparing multiple pair types [1]. Because of the different minimum
requirements for and particles, the variation in $\left<\mT\right>$ includes
both $\left<\pt\right>$ and particle mass differences. The anticorrelation of
$R_{\text{inv}}$ and $\left<\mT\right>$ has been interpreted as indicating the
presence of an expanding source [1].

Table 0.6 also includes the strong interaction scattering parameters $d_{0}$,
$\Re f_{0}$, and $\Im f_{0}$ obtained from the $\PGL\PKzS$ and $\PGL\PGL$
correlations. Figure 5 shows $d_{0}$ and $\Im f_{0}$ versus $\Re f_{0}$ in the
left and right panels, respectively, with the current results shown as red
stars and squares for $\PGL\PKzS$ and $\PGL\PGL$, respectively. The displayed
uncertainties are one-dimensional and are not based on a two-dimensional
contour.

Extracted values of the $R_{\text{inv}}$, $\Re f_{0}$, $\Im f_{0}$, $d_{0}$, ,
and $\left<\mT\right>$ parameters from the $\PKzS\PKzS$, $\PGL\PKzS$, and
$\PGL\PGL$ combinations in the 0–80% centrality range. The first and second
uncertainties are statistical and systematic, respectively.  Parameter
$\PKzS\PKzS$ $\PGL\PKzS$ $\PGL\PGL$ $R_{\text{inv}}$ (fm) $3.30\pm 0.10\pm
0.37$ $2.1^{+1.4}_{-0.5}\pm 0.7$ $1.3^{+0.4}_{-0.2}\pm 0.3$ $\Re f_{0}$ (fm)
$-0.76^{+0.29}_{-0.19}\pm 0.21$ $0.74^{+0.59}_{-0.16}\pm 0.34$ $\Im f_{0}$
(fm) $-0.07^{+0.48}_{-0.11}\pm 0.32$ $d_{0}$ (fm) $2.3^{+0.7}_{-0.5}\pm 1.3$
$4.2^{+5.7}_{-2.1}\pm 2.9$ $0.38\pm 0.02\pm 0.08$ $0.34^{+0.41}_{-0.12}\pm
0.16$ $1.5^{+1.2}_{-1.1}\pm 1.4$ $\left<\mT\right>$ () 1.53 2.09 2.60

Figure 5: The measured values of $d_{0}$ versus $\Re f_{0}$ (left) and $\Im
f_{0}$ versus $\Re f_{0}$ (right) from this analysis along with other
measurements and predictions as described in the text. For each data point,
the lines and the boxes indicate the (one-dimensional) statistical and
systematic uncertainties, respectively.

The negative value of $\Re f_{0}$ observed for the $\PGL\PKzS$ correlations,
combined with its relatively small magnitude, suggests a repulsive $\PGL\PKzS$
interaction. The uncertainty associated with the $\Im f_{0}$ value for the
$\PGL\PKzS$ correlations prevents any claim concerning possible inelastic
processes. The value of $\Re f_{0}$ found for $\PGL\PKzS$ correlations differs
from that reported by the ALICE Collaboration (teal diamonds) [15], which is
also for PbPb collisions but at $\sqrtsNN=2.76\TeV$. The uncertainties are too
large to determine if $d_{0}$ and $\Im f_{0}$ also differ between the two
results.

The positive $\Re f_{0}$ value obtained for the $\PGL\PGL$ correlations
suggests an attractive interaction that is not strong enough to produce a
bound state such as the H-dibaryon [9, 48]. This result disagrees with the
finding from the STAR Collaboration in gold-gold (AuAu) collisions at
$\sqrtsNN=200\GeV$ (blue circle). The negative $\Re f_{0}$ value of
${-}1.10\pm 0.37\stat^{+0.68}_{-0.08}\syst$fm found by STAR, combined with its
magnitude, imply a repulsive interaction. It is noted, however, that a
theoretical study of the STAR data which considers collective flow and feed-
down effects (shown as a shaded region at $d_{0}\approx 5$fm, $\Re
f_{0}\approx 0.9$fm) suggests that these data are consistent with the
$\PGL\PGL$ interaction being attractive [9]. An exclusion plot by the ALICE
Collaboration for the $\PGL\PGL$ scattering parameters obtained using the
$\PGL\PGL$ correlations from $\Pp\Pp$ collisions at $\sqrt{s}=7$ and 13, as
well as $\Pp\text{Pb}$ collisions at $\sqrtsNN=5.02\TeV$, also suggests an
attractive interaction [10]. In addition, our results are consistent with two
theoretical calculations (black triangles) that reproduce the $\PGL\PGL$
binding energy of ${}^{\;\;\;\;6}_{\PGL\PGL}{\mathrm{He}}$, as extracted from
the NAGARA event [49, 50].

##  0.7 Summary

The $\PKzS\PKzS$, $\PGL\PKzS$, and $\PGL\PGL$ femtoscopic correlations are
studied using lead-lead (PbPb) collision data at a center-of-mass energy per
nucleon pair of $\sqrtsNN=5.02\TeV$, collected by the CMS Collaboration. This
is the first report on $\PGL\PGL$ correlations in PbPb collisions at the CERN
LHC. The source size $R_{\text{inv}}$ and the incoherence parameter were
extracted for $\PKzS\PKzS$ correlations in six centrality bins covering the
0–60% range. The value of $R_{\text{inv}}$ decreases from 4 to 1fm going from
central to peripheral collisions and agrees with results from the ALICE
Collaboration at a similar transverse mass. Along with the $R_{\text{inv}}$
and parameters, the strong interaction scattering parameters, , the complex
scattering length and effective range, were extracted from $\PGL\PKzS$ and
$\PGL\PGL$ correlations in the 0–80% centrality range. These scattering
parameters indicate that the $\PGL\PKzS$ interaction is repulsive and that the
$\PGL\PGL$ interaction is attractive. The scattering parameters obtained from
$\PGL\PKzS$ correlations differ from those reported by the ALICE
Collaboration. The positive real scattering length obtained from the
$\PGL\PGL$ correlation disfavors the existence of a bound H-dibaryon state.
The $\PGL\PGL$ scattering parameters help to constrain baryon-baryon and, more
specifically, hyperon-hyperon interaction models. These measurements provide
an additional input to understand the nature of the strong interaction between
pairs of strange hadrons.

## References

  * [1] M. A. Lisa, S. Pratt, R. Soltz, and U. Wiedemann, “Femtoscopy in relativistic heavy ion collisions: Two decades of progress”,  Ann. Rev. Nucl. Part. Sci. 55 (2005) 357, 10.1146/annurev.nucl.55.090704.151533, arXiv:nucl-ex/0505014.
  * [2] V. G. J. Stoks, R. A. M. Klomp, M. C. M. Rentmeester, and J. J. de Swart, “Partial-wave analysis of all nucleon-nucleon scattering data below 350 MeV”,  Phys. Rev. C 48 (1993) 792, 10.1103/PhysRevC.48.792.
  * [3] J. J. de Swart and C. Dullemond, “Effective range theory and the low energy hyperon-nucleon interactions”,  Anna. Phys. 19 (1962) 458, 10.1016/0003-4916(62)90185-9.
  * [4] R. Engelmann, H. Filthuth, V. Hepp, and E. Kluge, “Inelastic -interactions at low momenta”,  Phys. Lett. 21 (1966) 587, 10.1016/0031-9163(66)91310-2.
  * [5] F. Eisele et al., “Elastic scattering at low energies”,  Phys. Lett. B 37 (1971) 204, 10.1016/0370-2693(71)90053-0.
  * [6] B. Sechi-Zorn, B. Kehoe, J. Twitty, and R. A. Burnstein, “Low-energy -proton elastic scattering”,  Phys. Rev. 175 (1968) 1735, 10.1103/PhysRev.175.1735.
  * [7] CMS Collaboration, “Bose–Einstein correlations in pp, pPb, and PbPb collisions at $\sqrtsNN=0.9\text{--}7\TeV$”,  Phys. Rev. C 97 (2018) 064912, 10.1103/PhysRevC.97.064912, arXiv:1712.07198.
  * [8] J. Schaffner-Bielich, M. Hanauske, H. Stöcker, and W. Greiner, “Phase transition to hyperon matter in neutron stars”,  Phys. Rev. Lett. 89 (2002) 171101, 10.1103/PhysRevLett.89.171101, arXiv:astro-ph/0005490.
  * [9] K. Morita, T. Furumoto, and A. Ohnishi, “$\PGL\PGL$ interaction from relativistic heavy-ion collisions”,  Phys. Rev. C 91 (2015) 024916, 10.1103/PhysRevC.91.024916, arXiv:1408.6682.
  * [10] ALICE Collaboration, “Study of the - interaction with femtoscopy correlations in pp and pPb collisions at the LHC”,  Phys. Lett. B 797 (2019) 134822, 10.1016/j.physletb.2019.134822, arXiv:1905.07209.
  * [11] R. L. Jaffe, “Perhaps a stable dihyperon”,  Phys. Rev. Lett. 38 (1977) 195, 10.1103/PhysRevLett.38.195.
  * [12] H. Takahashi et al., “Observation of a ${}^{\;\;\;\;6}_{\PGL\PGL}{\mathrm{He}}$ double hypernucleus”,  Phys. Rev. Lett. 87 (2001) 212502, 10.1103/PhysRevLett.87.212502.
  * [13] K. Nakazawa and H. Takahashi, “Experimental study of double- hypernuclei with nuclear emulsion”,  Prog. Theor. Phys. Supplement 185 (2010) 335, 10.1143/PTPS.185.335.
  * [14] Belle Collaboration, “Search for an $H$-dibaryon with a mass near $2{m}_{\PGL}$ in and decays”,  Phys. Rev. Lett. 110 (2013) 222002, 10.1103/PhysRevLett.110.222002, arXiv:1302.4028.
  * [15] ALICE Collaboration, “$\PGL\PK$ femtoscopy in Pb-Pb collisions at $\sqrtsNN=2.76\TeV$”,  Phys. Rev. C 103 (2021) 055201, 10.1103/PhysRevC.103.055201, arXiv:2005.11124.
  * [16] C. Loizides, J. Kamin, and D. d’Enterria, “Improved Monte Carlo Glauber predictions at present and future nuclear colliders”,  Phys. Rev. C 97 (2018) 054910, 10.1103/PhysRevC.97.054910, arXiv:1710.07098.
  * [17] R. Lednický and V. L. Lyuboshitz, “Final state interaction effect on pairing correlations between particles with small relative momenta”,  Sov. J. Nucl. Phys. 35 (1982) 770.
  * [18] HEPData record for this analysis, 2022.  10.17182/hepdata.133573.
  * [19] Tracker Group of the CMS Collaboration, “The CMS phase-1 pixel detector upgrade”,  JINST 16 (2021) P02027, 10.1088/1748-0221/16/02/P02027, arXiv:2012.14304.
  * [20] CMS Collaboration, “Track impact parameter resolution for the full pseudorapidity coverage in the 2017 dataset with the CMS phase-1 pixel detector”, CMS Detector Performance Note CMS-DP-2020-049, 2020.
  * [21] CMS Collaboration, “Performance of the CMS level-1 trigger in proton-proton collisions at $\sqrt{s}=13\TeV$”,  JINST 15 (2020) P10017, 10.1088/1748-0221/15/10/P10017, arXiv:2006.10165.
  * [22] CMS Collaboration, “The CMS trigger system”,  JINST 12 (2017) 01020, 10.1088/1748-0221/12/01/P01020, arXiv:1609.02366.
  * [23] CMS Collaboration, “The CMS experiment at the CERN LHC”,  JINST 3 (2008) S08004, 10.1088/1748-0221/3/08/S08004.
  * [24] CMS Collaboration, “Precision luminosity measurement in proton-proton collisions at $\sqrt{s}=13\TeV$ in 2015 and 2016 at CMS”,  Eur. Phys. J. C 81 (2021) 800, 10.1140/epjc/s10052-021-09538-2, arXiv:2104.01927.
  * [25] CMS Collaboration, “CMS luminosity measurement using nucleus-nucleus collisions at $\sqrtsNN=5.02\TeV$ in 2018”, CMS Physics Analysis Summary CMS-PAS-LUM-18-001, 2022.
  * [26] CMS Collaboration, “Charged-particle nuclear modification factors in PbPb and pPb collisions at $\sqrtsNN=5.02\TeV$”,  JHEP 04 (2017) 039, 10.1007/JHEP04(2017)039, arXiv:1611.01664.
  * [27] CMS Collaboration, “Description and performance of track and primary-vertex reconstruction with the CMS tracker”,  JINST 9 (2014) P10009, 10.1088/1748-0221/9/10/p10009, arXiv:1405.6569.
  * [28] CMS Collaboration, “Transverse-momentum and pseudorapidity distributions of charged hadrons in pp collisions at $\sqrt{s}=0.9$ and 2.36 TeV”,  JHEP 02 (2010) 041, 10.1007/JHEP02(2010)041, arXiv:1002.0621.
  * [29] CMS Collaboration, “Observation and studies of jet quenching in PbPb collisions at nucleon-nucleon $\sqrtsNN=2.76\TeV$”,  Phys. Rev. C 84 (2011) 024906, 10.1103/PhysRevC.84.024906, arXiv:1102.1957.
  * [30] C. Gale, S. Jeon, and B. Schenke, “Hydrodynamic modeling of heavy ion collisions”,  Int. J. Mod. Phys. A 28 (2013) 1340011, 10.1142/S0217751X13400113, arXiv:1301.5893.
  * [31] GEANT4 Collaboration, “—a simulation toolkit”,  Nucl. Instrum. Meth. A 506 (2003) 250, 10.1016/S0168-9002(03)01368-8.
  * [32] CMS Collaboration, “Strange hadron collectivity in pPb and PbPb collisions”, 2022. arXiv:2205.00080.
  * [33] CMS Collaboration, “Strange particle production in pp collisions at $\sqrt{s}=0.9$ and 7 TeV”,  JHEP 05 (2011) 064, 10.1007/JHEP05(2011)064, arXiv:1102.4282.
  * [34] CMS Collaboration, “Long-range two-particle correlations of strange hadrons with charged particles in pPb and PbPb collisions at LHC energies”,  Phys. Lett. B 742 (2015) 200, 10.1016/j.physletb.2015.01.034, arXiv:1409.3392.
  * [35] Particle Data Group Collaboration, “Review of particle physics”,  Prog. Theor. Exp. Phys. 2022 (2022) 083C01, 10.1093/ptep/ptac097.
  * [36] H. Voss, A. Höcker, J. Stelzer, and F. Tegenfeldt, “TMVA, the toolkit for multivariate data analysis with ROOT”, in  XIth International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT), p. 40.  2007\.  arXiv:physics/0703039.  10.22323/1.050.0040.
  * [37] G. I. Kopylov, “Like particle correlations as a tool to study the multiple production mechanism”,  Phys. Lett. B 50 (1974) 472, 10.1016/0370-2693(74)90263-9.
  * [38] A. Kisiel, “Non-identical particle correlation analysis in the presence of non-femtoscopic correlations”,  Acta Phys. Polon. B 48 (2017) 717, 10.5506/APhysPolB.48.717.
  * [39] ALICE Collaboration, “pp, p-, and - correlations studied via femtoscopy in pp reactions at $\sqrt{s}=7\TeV$”,  Phys. Rev. C 99 (2019) 024001, 10.1103/PhysRevC.99.024001, arXiv:1805.12455.
  * [40] STAR Collaboration, “Neutral kaon interferometry in Au+Au collisions at $\sqrtsNN=200\GeV$”,  Phys. Rev. C 74 (2006) 054902, 10.1103/PhysRevC.74.054902, arXiv:nucl-ex/0608012.
  * [41] A. D. Martin and E. N. Ozmutlu, “Analyses of $\PK\PK$ production and scalar mesons”,  Nucl. Phys. B 158 (1979) 520, 10.1016/0550-3213(79)90180-9.
  * [42] A. Antonelli, “Radiative decays”,  eConfC 020620 (2002) THAT06, 10.48550/ARXIV.HEP-EX/0209069, arXiv:hep-ex/0209069.
  * [43] N. N. Achasov and V. V. Gubin, “Analysis of the nature of the $\PGf\to\PGg\PGp\Pgh$ and $\PGf\to\PGg\PGpz\PGpz$ decays”,  Phys. Rev. D 63 (2001) 094007, 10.1103/PhysRevD.63.094007, arXiv:hep-ph/0101024.
  * [44] N. N. Achasov and A. V. Kiselev, “New analysis of the KLOE data on the $\PGf\to\PGh\PGpz\PGg$ decay”,  Phys. Rev. D 68 (2003) 014006, 10.1103/PhysRevD.68.014006, arXiv:hep-ph/0212153.
  * [45] F. James and M. Roos, “Minuit: A system for function minimization and analysis of the parameter errors and correlations”,  Comput. Phys. Commun. 10 (1975) 343, 10.1016/0010-4655(75)90039-9.
  * [46] PHENIX Collaboration, “Lévy-stable two-pion Bose–Einstein correlations in $\sqrtsNN=200\GeV$ Au+Au collisions”,  Phys. Rev. C 97 (2018) 064911, 10.1103/PhysRevC.97.064911, arXiv:1709.05649.
  * [47] ALICE Collaboration, “One-dimensional pion, kaon, and proton femtoscopy in Pb-Pb collisions at $\sqrtsNN=2.76\TeV$”,  Phys. Rev. C 92 (2015) 054908, 10.1103/PhysRevC.92.054908, arXiv:1506.07884.
  * [48] STAR Collaboration, “$\PGL\PGL$ correlation function in Au+Au collisions at $\sqrtsNN=200\GeV$”,  Phys. Rev. Lett. 114 (2015) 022301, 10.1103/PhysRevLett.114.022301, arXiv:1408.4360.
  * [49] E. Hiyama et al., “Four-body cluster structure of $\mathrm{A}=7\text{--}10$ double- hypernuclei”,  Phys. Rev. C 66 (2002) 024007, 10.1103/physrevc.66.024007, arXiv:nucl-th/0204059.
  * [50] I. N. Filikhin and A. Gal, “Faddeev–Yakubovsky calculations for light $\PGL\PGL$ hypernuclei”,  Nucl. Phys. A 707 (2002) 491, 10.1016/s0375-9474(02)01008-4, arXiv:nucl-th/0203036.