A femtoscopic correlation analysis tool using the Schrödinger equation (CATS)

We present a new analysis framework called “Correlation Analysis Tool using the Schrödinger equation” (CATS) which computes the two-particle femtoscopy correlation function C(k), with k being the relative momentum for the particle pair. Any local interaction potential and emission source function can be used as an input and the wave function is evaluated exactly. In this paper we present a study on the sensitivity of C(k) to the interaction potential for different particle pairs: p–p, p–Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\Lambda }$$\end{document}, K-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {K^-}$$\end{document}–p, K+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {K^+}$$\end{document}–p, p–Ξ-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\Xi }^-$$\end{document} and Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\Lambda }$$\end{document}–Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\Lambda }$$\end{document}. For the p–p Argonne v18\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{18}$$\end{document} and Reid Soft-Core potentials have been tested. For the other pair systems we present results based on strong potentials obtained from effective Lagrangians such as χ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi $$\end{document}EFT for p–Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\Lambda }$$\end{document}, Jülich models for K(K¯)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{K}(\bar{\mathrm{K}})$$\end{document}–N and Nijmegen models for Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\Lambda }$$\end{document}–Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\Lambda }$$\end{document}. For the p–Ξ-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\Xi }^-$$\end{document} pairs we employ the latest lattice results from the HAL QCD collaboration. Our detailed study of different interacting particle pairs as a function of the source size and different potentials shows that femtoscopic measurements can be exploited in order to constrain the final state interactions among hadrons. In particular, small collision systems of the order of 1 fm, as produced in pp collisions at the LHC, seem to provide a suitable environment for quantitative studies of this kind.


Introduction
Femtoscopy has been mainly employed so far to study the properties of the particle emitting source in heavy-ion collisions by analyzing particle pairs with low relative momentum undergoing a known interaction [1,2]. In heavy-ion collisions pion femtoscopy has been the most common tool to get insight into the time-space evolution of the produced medium [3][4][5][6][7][8][9][10][11][12][13]. a e-mail: dimitar.mihaylov@mytum.de b e-mail: valentina.mantovani-sarti@tum.de Unlike heavy-ion collisions, nucleon-nucleon (NN) collisions are not affected by the formation of a medium such as the Quark-Gluon-Plasma. This means that the time evolution of the source can be neglected and the hadron interaction is not influenced by collective processes.
More generally, femtoscopy allows to describe any hadronhadron correlation involving both (anti)mesons and (anti) baryons. Recent femtoscopy studies in pp at √ s = 7 TeV and p-Pb at √ s NN = 5.02 TeV showed that in these systems a smaller source size radius is extracted compared to heavy-ion collisions [14][15][16]. When mesons are considered the correlation function is affected by a mini-jet background over the whole k range [14,15]. This effect is also visible in baryon-antibaryon correlations (B-B), where partons from the colliding protons hadronize in particle cones containing a B-B pair. The production of such pairs is enhanced due to the intrinsic baryon number conservation and therefore shows a strong kinematic correlation which overlaps to the true femtoscopic signal in the final correlation function. Since these mini-jets are not present in the correlation function of B-B/B-B [17], those systems are well suited to study the final state interaction.
The femtoscopic formalism includes the computation of the correlation function for a given source function and interaction potential [2] where k = (| p 1 − p 2 |)/2 is the reduced relative momentum in the center of mass of the pair (p 1 +p 2 = 0), r is the relative distance between the two particles, S(r) is the source function and ψ k (r) is the two-particle wave function.
The above formula is based on the assumptions that space and momentum emission coordinates are uncorrelated and the emission process is time independent [2]. The dynamics of the processes is explicitly embedded on one hand in the Fig. 1 The probability density function of the relative distance r for p-p pairs produced in pp collisions at 7 TeV. The blue solid line is the prediction from the EPOS transport model, the dashed orange line, the green dot-dashed line and the dotted red line are Gaussian sources with different radii. For better visibility the source with r 0 = 0.55 fm is scaled down by 0. 6. The red shaded area shows the typical range of the strong interaction where an approximate solution to the Schrödinger equation may fail emitting source S(r), which depends on the colliding system, and on the other hand in the interaction potential between the two particles expressed by the relative wave function ψ k (r).
The source could be either parametrized using a specific distribution function, e.g. a Gaussian or a Cauchy function, or extracted directly from transport models such as EPOS [18]. The pair wave function can be obtained by solving the Schrödinger equation for a specific interaction potential and a popular tool to do that is the Correlation Afterburner (CRAB) [19]. This afterburner only delivers the asymptotic solution at large relative distances which eventually leads to a wrong correlation function at small distances. For heavyion collisions the source radius is typically between 2 and 6 fm [3][4][5][6][7]. However in pp collisions at LHC energies the extracted source can be even smaller than 1.5 fm [14,15], and the EPOS transport model predicts a non-Gaussian source which peaks at distances of around 1 fm (see Fig. 1).
This calls for a tool capable of providing an exact solution of the Schrödinger equation valid also at short distances.
Another common approach based on solving the Schrödinger equation and used to study p-p correlations is the Koonin model [20]. Here a Gaussian source distribution is assumed and the relative wave function is obtained by employing a Coulomb potential and a Reid Soft-Core potential for the strong interaction [21].
Analytical solutions to compute C(k) do exist, but they are limited by different approximations. One of the most popular analytical models was developed by Lednický and Lyuboshits [22]. The emission source is assumed to have a Gaussian shape and the wave function is modeled within the effective range approximation, using the scattering length (a 0 ) and the effective range of the potential (r e f f ). Such a parametrization is only sensitive to the asymptotic region of the interaction. Due to its long range nature, the Coulomb potential is normally not considered in the Lednický model and can at most be treated in an approximate way [13].
All the above-mentioned issues motivated the development of a new femtoscopy tool called "Correlation Analysis Tool using the Schrödinger equation" (CATS).
CATS is capable of evaluating numerically the full wave function without using any approximations and computing C(k) for different source functions and interaction potentials. CATS is faster than most of the conventional numerical tools and also more flexible since it can be used together with an external fitter.
The main purpose of this work is to test the sensitivity of the correlation function to different potentials and provide qualitative examples.
In particular for p-p and p-interactions we only investigate the most common potentials used in femtoscopy so far and compare them to more up-to-date potentials.
The p-p correlation function has been studied by comparing the Reid Soft-Core [21] potential adopted in most femtoscopic analyses to the most recent Argonne v 18 [23] potential which provides a better description of scattering data.
Similarly, for the p-combination we compare the commonly used Usmani potential [24] to the most recent NLO potential obtained in χ EFT calculations [25] which significantly improves the agreement with available data.
For the K ± -p interaction we employ the Jülich model [26,27] and we investigate the interplay between the Coulomb and the strong interaction in the correlation function.
At the moment, the limited amount of experimental data for the -interaction do not allow to discriminate among different potentials. In this paper we show that future femtoscopic analyses for the -correlation function performed with CATS open the possibility to differentiate between different types of potentials.
For this pair, following the results presented in [28], we adopt four potentials: two versions of hard-core Nijmegen models (ND56 [29], NF50 [30]), a quark-model potential with meson exchange effect included (fss2 [31]) and another hard-core Nijmegen potential NF42 [30] which allows for a bound state.
The first three ones have been chosen since they deliver the best agreement with the recent -STAR measurements [12] while the latter has been selected since it is connected to a possible bound state and it shows a different behavior as a function of the relative distance r between the two particles.
For p-pairs we adopt the most recent potential obtained by the HAL QCD lattice collaboration [32,33].
The most modern potentials can also be implemented in CATS, however this would require their transformation to a local form, nevertheless this is beyond the scope of the current study.
This paper is organized as follows. In Sect. 2.1 we introduce the formalism of the Schrödinger equation as embedded in the algorithm, while in Sect. 2.2 we introduce the interaction potentials used for the later investigations of the p-p, p-, K(K)-p, p-− and -systems.
In Sect. 3 we show our results for the correlation functions of the above-mentioned pairs using different interaction potentials and source sizes. A comparison with the emitting source obtained from the EPOS transport model is presented for p-p and p-pairs. In addition we show results for p-− and -correlation functions including experimental effects such as momentum resolution and feed-down decays. Moreover we present a comparison between the extracted source radius for p-p pairs obtained from the HADES collaboration in p-Nb reactions at √ s NN = 3.18 GeV and a CATS fit of the same data.
Finally in Sect. 4 we present conclusions and future outlooks. In addition in Appendix A we provide technical information related to the installation of CATS and the numerical methods used in the code.

Overview
CATS is designed to compute the two-particle correlation function for any emission distribution and interaction potential. This is achieved by numerically solving the Schrödinger equation (SE) and by evaluating the convolution of the resulting wave function with a source distribution (see Eq. (1)).
The two-particle stationary SE reads where μ = m 1 m 2 m 1 +m 2 is the reduced mass of the system. By assuming a central interaction potential V , the total wave function is separated into a radial term and an angular term given by the spherical harmonics In particular the radial equation to be solved is where u(r ) = r R(r ).
The overall interaction potential V (r ) is the sum of a short range strong contribution and a long range Coulomb potential. Since we are looking for scattering states as a function of the relative momentum k and not for bound states, the energy E of the system needs to be fixed by using the relation E =h 2 k 2 2μ . The scattered state will then approach the asymptotic solution outside the range of the strong potential leading to a free or a pure Coulomb wave function depending on the charge of the particles involved.
In order to properly match the exact solution to the asymptotic form we have to evaluate the phase shifts of the corresponding wave function. For this purpose the SE is solved by expanding the total wave function in partial waves where P l (cos θ) are the Legendre polynomials. The sum over the partial waves runs until the convergence of the solution is reached, and from there we obtain the total wave function ψ k (r) to be used in Eq. (1).
Taking into account the Fermi statistics of the particle pairs and fixing the isospin (I ) configuration, only specific states 2S+1 L J are allowed depending on the particle species, where we use the standard spectroscopic notation with S denoting the total spin, L the orbital angular momentum and J the total angular momentum of the pair. Moreover, since NN, nucleon-hyperon (NY), kaon-nucleon (KN) and hyperonhyperon (YY) potentials might involve spin-orbit couplinḡ L ·S, the total angular momentumJ =L +S has to be taken as a good quantum number to characterize the eigenstates and has to be accounted for in the total degeneracy of the states.
The total correlation function that CATS provides as an output is then given by where each C (I,S,L ,J ) is evaluated by means of Eq. (1) and weighted by w (I,S,L ,J ) . In experimental conditions with unpolarized particles we should take into account the degeneracy in S, as well as in I . Moreover for L > 0 states also the degeneracy in J has to be considered. How to compute the corresponding weights w (I,S,L ,J ) is explained in Appendix A.
In CATS the source function is characterized in two ways. One possibility is to define an analytical function which models the emission source. It is typically based on a Gaussian distribution of the x,y and z single particle coordinates. The emission source S(r, θ, φ) is defined as the probability density function (pdf) to emit a particle pair at a certain relative distance r and relative polar and azimuthal angles θ and φ. For an uncorrelated particle emission the Gaussian source reads where r 0 is the size of the source. Since no angular dependence is involved, to obtain the probability of emitting two particles at a distance r a trivial integration over the solid angle is necessary Another possibility is to sample the relative distance between the particle pairs directly from the output of an event generator such as EPOS.
In Fig. 1 we compare the shapes of differently sized Gaussian sources to the EPOS source for p-p pairs obtained in a simulation of pp collisions at 7 TeV. Notably the EPOS source predicts that most of the particle pairs will be emitted at distances below 2 fm, which is the typical range of the strong interaction. Hence the asymptotic solution of the wave function will no longer be valid in that region, highlighting the necessity of an exact treatment of the problem. Moreover we see that EPOS predicts a non-Gaussian source.

NN and NY interaction potentials
The NN interaction has been widely constrained thanks to a large amount of scattering data along with data on bound states such as the deuteron [34]. In the common Yukawa picture the NN potential is described by means of model-independent components as the long range One-Pion-Exchange term (OPE) and the ππ exchange term for the intermediate attraction [35]. The vector meson (typically ω) exchange accounts for the repulsive core.
This repulsive contribution can be phenomenologically modeled with a Woods-Saxon function.
All these features are included in the most common NN potentials such as Bonn [36,37], Urbana [38], Paris [39], Nijmegen [40], Reid [21] and Argonne [23] potentials. In addition, the results obtained within a chiral effective field theory approach (χ EFT) [41][42][43][44][45] and within the latest version of the Nijmegen model ESC08 [46] have to be mentioned, since they provide similar results for the scattering parameters and both models show a nice description of the available data. The phase shifts are taken from [34] For the study of the p-p correlation function in this paper we employ the Argonne v 18 [23] and the Reid Soft-Core (RSC) [21] potentials. The latter is the default choice in femtoscopic analyses [19,20]. In the Koonin model only the RSC singlet state 1 S 0 is considered, while in the CRAB calculations also the triplet 3 P 2 state is included but without spinorbit or tensor terms. Nevertheless, these terms are crucial in order to reproduce the 3 P 2 phase shift, as can be seen in Fig. 2. Here the Argonne v 18 , which includes all s-and pwave states and both the spin-orbit and tensor terms, results in a much better agreement with the partial wave data in the momentum range we are interested in.
It should be noted that a more recent version of the RSC potential [47] is available. Nevertheless it is beyond the purpose of the current work to study in details the subtle differences between a large variety of NN potentials, since we do not expect a large discrepancy in the p-p correlation function.
The range of relative momenta k relevant for femtoscopic studies stops at k max ≈ 200 MeV/c and the dominant contribution in the phase shifts in this regime is mainly coming from s-and p-wave states [34]. For this reason at the moment we are neglecting the coupled channel to the 3 F 2 wave.
Unfortunately in the nucleon-hyperon sector only 36 scattering data points have been measured [50][51][52][53][54][55] and the available phenomenological potentials are fitted in order to reproduce these data along with the accessible hypernuclei binding energies. In particular, for the p-interaction we adopt the potential introduced by Bodmer, Usmani, and Carlson [24], which reproduces the available cross section data as shown in Fig. 3. This potential is similar to a NN Argonne interaction since it includes a repulsive Woods-Saxon core and a two-pion intermediate exchange term. At the moment we just consider s-waves, which account for both the singlet 1 S 0 and the triplet 3 S 1 state.  Table 1 The p-scattering parameters of the Usmani and the chiral potentials. We denote with (s) and (t) the singlet and triplet state respectively. The scattering parameters for the Usmani potential are evaluated here [24]. The parameters for the the chiral potential are taken from [25] for a cutoff scale of 600 MeV To test the sensitivity of the correlation function to the pinteraction, we have also used the results obtained by nextto-leading order (NLO) χ EFT calculations [25]. In this work the authors show that the agreement with the available data in the S = − 1 sector is considerably improved with respect to previous leading order results. A local potential for this model is not yet available, nevertheless we were provided with the radial wave function and used CATS to compute the resulting correlation functions for different sources. The scattering parameters of the two potentials used to model the p-interaction are listed in Table 1.
Concerning the S = − 2 sector the interaction of hyperons with nuclei is not well constrained experimentally or phenomenologically [56]. Recently a -hypernucleus candidate has been detected [57] and ongoing measurements suggest that the N-interaction is weakly attractive [58].
Recent lattice QCD calculations from the HAL QCD collaboration showed preliminary results on the p-− correlation function based on local potentials in the s-wave singlet and triplet states, both in I = 0 and I = 1 channels [32,33]. The interaction potentials have been obtained in (2 + 1)flavor lattice QCD simulations close to the physical masses of pions and kaons. In particular in the isoscalar 1 S 0 channel it has been shown that the p-− interaction is deeply attrac-tive at intermediate distances and relatively weakly repulsive at shorter distances. The isoscalar spin triplet state 3 S 1 has a shallow attractive well and a very repulsive core.
The I = 1 channel does not present any attraction in the spin singlet state while a mildly attractive region enclosing a repulsive inner core has been found in the spin triplet state. These contributions might not be relevant for large source sizes [32] but could impact the overall correlation function for smaller emitting sources as obtained in elementary collisions.

KN andKN interaction
As for the hyperon-nucleon case, the (anti)kaon-nucleon system is a powerful tool to study the short range strong interaction with the exchange of strangeness as a new degree of freedom.
In particular kaons are suited to study the inner part of nuclei since at low energies they interact rather weakly with the nucleons [56]. On the contrary antikaons scattering on nuclei is characterized by absorption processes that can lead to the production of hyperons such as and [56,59].
From a theoretical point of view, based on the successful description of the NN interaction, several meson-baryon exchange models, such as the Jülich model [26,27], have been applied so far in the description of the available KN data. In this approach the short range repulsion in K + -N is not well reproduced by the simple presence of the ω meson. In order to reproduce the s-wave scattering parameters an additional repulsive contribution is needed. The inclusion of spin-dependent terms, such as one gluon exchange (OGE) or exchange of the scalar-isovector meson a(980), improves the agreement in the s-channel [67].
It is also worth mentioning that different theoretical approaches based on quark models [68][69][70] and χ EFT interactions [71] have also been used to investigate the kaonnucleon interaction.
In this work we adopted the above mentioned Jülich model, which besides single particle exchange terms also contains box diagrams allowing for intermediate states such as N, , K(K), K * (K * ). A local potential was not available for this model, nevertheless we were provided with the exact solution of the radial wave function, which was included in CATS for the computation of the correlation function.
Besides the wave function corresponding to the strong K(K)-N interaction potential, we also take into account the Coulomb interaction in an approximate way by multiplying the (strong) wave function with the Gamow factor A c (η) [13] Table 2 The potential parameters used in the two-range Gaussian parametrization in Eq. (10) and the corresponding scattering parameters a 0 (fm) and r e f f (fm). The values are taken from [28] Model where η = 1 2 a s k KN,KN −1 and a s = (m red z 1 z 2 e 2 ) −1 , with z 1,2 being the charge numbers of the two particles. A more elaborate approximation based on the renormalization of the asymptotic behaviour of the strong wave function has been developed in the past [72] and will be used in future analyses.

YY interaction
In the YY sector the available experimental data on the binding energy of hypernuclei, despite recent data on 6 He [73,74], do not allow to set tight constraints on the nature of the underlying interaction [56].
Recently the STAR collaboration employed the femtoscopy technique to study -correlations in Au-Au collisions at √ s NN = 200 GeV [12]. The reported shallow repulsive interaction in this work is still under debate, since in an alternative analysis, that considered the contribution of the residuals in a more sophisticated way, a shallow attractive interaction was confirmed [28]. One of the goals of this work is to show that a detailed study of the -correlation function might be a sensitive tool to extract information on the interaction of the two particles.
For this purpose, we select four different potentials: three meson-exchange Nijmegen models with a short range hardcore (ND56 [29], NF50 and NF42 [30]) and one quark-model potential which includes meson exchange effects as well (fss2 [31]). The models ND56, NF50 and fss2 deliver the best agreement with the -correlation data measured by STAR. The NF42 has been chosen to show the sensitivity of the correlation function to the presence of a possible bound state, which is allowed by this potential.
The above mentioned potentials have been translated in a local form by using a two-range Gaussian parametrization, following the approach used in [28] The parameters are fixed to values that reproduce the scattering parameters obtained within each model. They are sum-  10) is taken from [28] marized in Table 2. In Fig. 4 the potentials are plotted as a function of the relative distance r . It is worthwhile mentioning that only the NF42 model provides a positive scattering length a 0 ≈ 3.7 fm due to the presence of a bound state, while all the remaining potentials lead to negative (attractive) values of a 0 and effective ranges between 4 and 5 fm.

Results
In Fig. 5 a comparison between the correlation functions obtained by CRAB and CATS for p-p pairs is shown. In particular two different transport models are used to model the source: UrQMD [75] for simulating p-Nb reactions at √ s NN = 3.18 GeV and EPOS [18] for simulating pp reactions at √ s = 7 TeV. The p-Nb system is expected to have a Gaussian source of around 2 fm [76], while the pp system produces a much smaller source, comparable to a Gaussian source of below 1 fm (see Fig. 1). As already pointed out the CRAB tool provides only an approximate solution at small distances. Figure 5 shows that in the case of the larger source (the p-Nb system) the two models are in agreement within few percent. However for pp collisions at 7 TeV the source is much smaller and consequently a much larger discrepancy, of up to 20%, is observed between CRAB and the exact CATS solution. In Fig. 6 we show the correlation function C(k) obtained by CATS for p-p pairs (upper panel) and p-pairs (lower panel) with different interaction potentials and Gaussian source sizes.
In particular for p-p pairs we perform the calculations for the RSC potential as used in the Koonin and CRAB models and compare them to the Av 18 potential. The profile of the resulting correlation functions reproduce the main features coming from the interplay between the attractive part of the NN potential and the short-range repulsive core along with contributions from Coulomb and Fermi-Dirac statistics. The repulsive interaction is visible between 60 < k < 200 MeV since C(k) is below unity. This effect is enhanced for smaller sources.
The strength of the p-p correlation signal increases by almost a factor 2.5 going from r 0 = 2 fm to r 0 = 0.85 fm. The differences among the potentials are negligible for most of the source sizes. Nevertheless it is evident that for sources below 1 fm deviations up to 10% are present at k > 150 MeV.
The p-pair is neither affected by Coulomb interaction nor by quantum statistics. As a result the shape of the corresponding correlation function C(k) is only affected by the strong interaction potential. Here results are shown for the Usmani potential and the χ EFT NLO potential (see Sect. 2.2). Both of these potentials possess a phenomenological repulsive core which is dominated by an intermediate attraction as k decreases. This results in a rise of C(k) at lower momenta. The two potentials produce very similar correlation functions and the differences are mostly associated with the slightly more attractive 3 S 1 state in the Usmani potential. This can be easily understood by comparing the scattering parameters (see Table 1) of the Usmani and the chiral potential. While the scattering lengths in the singlet state are almost identical, they do differ in the triplet state, where the Usmani potential is slightly more attractive, resulting in an enhanced correlation.
In addition for the NLO correlation function we compare the exact solution to the approximate Lednický model. As expected we observe a significant deviation between the exact and approximate solution at small source radii. The Lednický model produces an enhanced correlation signal because it contains the NLO scattering parameters that account only for the average attractive interaction. On the contrary the correlation function obtained feeding the exact NLO wave function in CATS shows the repulsive component expected within this model.
It is clear that for sources ≥ 2 fm the correlation function is not sensitive to the repulsive core both for p-p and ppairs. However, in the case of smaller sources, as shown in the last panel of Fig. 6, the correlation functions split up and remain different even for relative momenta k > 150 MeV. This region is commonly used in femtoscopy analyses as a correlation free baseline to which the overall normalization of C(k) is applied. Due to to the possible non-flat correlation in that region particular care has to be considered in including the correct potential in the correlation function.
In Fig. 7 the emitting source has been directly taken from the EPOS model simulating pp collisions at √ s = 7 TeV. In this transport code the evaluated source peaks at very low distances and has a long tail, as already shown in Fig. 1. Moreover the EPOS source presents a non-trivial angular depen-dence of the emission source. The sources for both the p-p and p-systems in EPOS are the same. However, if a Gaussian parametrization is employed to reproduce the correlation functions obtained with EPOS, different source sizes r 0 for both systems are needed. This suggests that when investigating experimental baryon-baryon correlations in small systems, particular care has to be taken in the assumption of the profile of the emitting source.
Results for K − -p and K + -p correlations for different source sizes are shown in the upper and lower panel of Fig. 8. The correlation function including only the Coulomb interaction is plotted in order to see the effect of the strong potential.
The K − -p correlation function is deeply affected by the strong interaction already at k ≈ 200 MeV and it dominates the low k region. As expected both these effects are enhanced when the source gets smaller.
The opposite scenario is depicted for the K + -p pair, where the Coulomb and the strong repulsion lead to an anticorrelated signal at low momenta.
In Fig. 9 we show the expected theoretical p-− correlation function. We compare C(k) for a pure Coulomb interaction with the correlation function including the preliminary strong potential from the HAL QCD collaboration. The included channels are the (I = 0, S = 0), (I = 0, S = 1), (I = 1, S = 0) and (I = 1, S = 1), where only s-waves are considered. There is a clear enhancement of the correlation signal at small relative momenta k, which is a result of the attractive I = 0 channel. This effect is further increased for small emitting sources.
In Fig. 10 we present results for -analysis with several interaction potentials (see Sect. 2.4). As a comparison the correlation function obtained only with the quantum statistics included is presented. The latter is responsible for the limit value C(k) = 0.5 as the relative momentum k tends to zero.
The correlation function shows a strong sensitivity to the interaction potential and moreover the shape of C(k) is deeply affected by the source size. This behaviour can  Fig. 4 and their corresponding scattering parameters listed in Table 2. The most binding potentials, NF42 and ND56, lead to a correlation function that is significantly enhanced for smaller sources as it becomes sensitive to the deep attractive contribution. The weakly attractive potentials, NF50 and fss2, also show an enhancement compared to the quantum statistics baseline, nevertheless the correlation signal is much weaker due to the presence of a large repulsive core.
When making predictions about the correlation function it is important to take into account the expected experimental limitations. In particular the momentum resolution and the residual correlations contributing to the femtoscopic signal can distort the correlation function [77]. To estimate how the p-− and -correlation functions would be affected by experimental corrections, we use the linear decomposition of the correlation function as described in [77]. In short one assumes that only a fraction λ of the particle pairs represent the actual correlation of interest and the rest of the pairs are either misidentified or stemming from the decays of resonances (feed-down). For the purpose of a qualitative prediction we will assume that all feed-down contributions have a flat correlation distribution. In that case the experimental correlation function C ex p (k) is expressed by In Fig. 11 we plot the estimated experimental correlation functions for p-− and -pairs. For both systems we apply the momentum smearing matrix used for the p-correlation in [77]. The λ coefficient for the -correlation is taken from the same reference and has the value of 32%. To obtain a reasonable λ coefficient for p-− we use the proton feed-down and impurities as evaluated in [77], while for the − particle we assume 90% purity and feed-down fractions of 32% from (1530) [78] and ≈ 1% from − . The corresponding λ coefficient for the genuine p-− correlation has a value of 52%. One can see in Fig. 11 that for the -case the effect of the different potentials on the correlation function is significantly damped by the inclusion of expected experimental effects. This implies that very large statistics are needed to reach a high precision in the determination of the scattering parameters or in testing different models. For the p-− case the effect of the strong potential is clearly distinguish-  In order to apply the CATS framework to experimental data we present, in Fig. 12, the result for the p-p correlation function obtained in p-Nb collisions at √ s NN = 3.18 GeV [76]. We evaluated C(k) in CATS using the Av 18 potential with s-and p-waves included and performed a fit to the data by multiplying C(k) with a normalization constant N . The source is assumed to be Gaussian and the source size is a free fit parameter. The HADES analysis extracts a radius of r 0 = 2.02 ± 0.01(stat) fm and our fit value is 2.01 ± 0.01(stat) fm, which is in perfect agreement with the published data.

Conclusions and outlooks
In this paper we presented the new femtoscopy analysis tool CATS, which allows for an exact computation of the correlation function for different sources and potentials. To obtain a solution valid for any source size CATS solves the Schrödinger equation exactly, thus providing the possibility to investigate not only heavy-ion reactions but small collision systems as well.
In the present work we show the theoretical correlation functions for the p-p, p-, K(K)-p, p-− and -particle pairs, by using different models for the interaction potential.
A comparison between CATS and CRAB reveals that for an emitting source of around 2 fm the correlation functions are quite similar. Nevertheless, for small colliding systems, such as pp at TeV energies, the exact treatment of the Schrödinger equation strongly deviates from the approximate solution of CRAB. Thus CATS provides a unique opportunity to investigate the short-range interaction between particles in small systems.
We show that both the p-p and p-correlation functions are affected by the repulsive core of their interaction potentials, which is mostly evident for small emitting sources. This effect is particularly relevant for momenta above 100 MeV, where the exact solution deviates significantly compared to approximate models like Lednický.
Moreover for the p-p and p-correlations we test different emission sources by using the EPOS transport model. In Fig. 1 we show that the emitting source provided by EPOS is completely different from the traditionally Gaussian form used so far in femtoscopic analyses. The results presented here suggest the need for a more detailed analysis of the profile of the emitting source beyond the traditional Gaussian parametrization, whenever dealing with small systems.
The investigation of the K(K)-p correlation function shows a relevant modification due to the strong interaction. This effect increases with a decreasing source size.
The p-− correlation function is calculated by employing a preliminary local potential obtained from recent lattice HAL QCD results. We observe that for small emitting sources, as expected in pp collision systems, the strong interaction modifies the correlation function significantly. This implies that future experimental femtoscopic data on the p-− correlation should be sensitive to this effect. This statement is still valid if one includes momentum resolution corrections and residual effects in the correlation function.
For the -correlation function we show results for different interaction potentials, based on Nijmegen bosonexchange models and quark models. Starting from source sizes of r 0 = 2 fm, the corresponding C(k) shows a relevant sensitivity to the strong potentials. As the source radius decreases the correlation function is modified by the attractive contribution common to all of the above-mentioned potentials. We also show that the inclusion of the expected experimental effects makes the separation of the different potentials extremely challenging and will require a lot of statistics.
CATS was designed such that it can easily be used with an external fitter. To demonstrate this we perform a refit of the pp correlation function measured by the HADES experiment in p-Nb reactions at √ s NN = 3.18 GeV. The obtained results are in a perfect agreement with the published data.
The experimental database of femtoscopy is constantly growing, thus in the near future it will be essential to perform global analyses over all available data in order to get the best possible constraints on the interaction potentials. The flexibility of CATS allows it to be used as the base of more complex frameworks and is therefore ready to fulfill the requirements of the upcoming experimental analyses. like to focus on the two most important inputs, namely the emission source and the interaction potential.
Appendix A.1: The source function Commonly the source function is considered to have only r dependence, r being the relative distance between the two particles. However CATS uses a more generic definition which allows an additional dependence on the relative momentum k and the scattering angle θ = (r, k). The specific implementation of the source and its parameters in CATS allows to link it to an external framework to perform data-fitting (see Fig. 12).
Alternately the source can be inserted into CATS from an OSCAR 1997 A file [79]. This format is commonly used by transport codes in relativistic heavy ion collisions and it contains event-by-event information about the position and the momentum four-vectors of all particles produced in the simulation. Note that obtaining the source from simulations leads to statistical uncertainties which result in a limited precision when determining the source function. Thus CATS also computes the uncertainty associated with the source function and propagates it to the correlation function.

Appendix A.2: Computing the wave function
The heart of CATS is the numerical solver for computing the wave function. The two body problem can be reduced to a one-body problem simply by using the reduced mass μ = (m 1 m 2 )/(m 1 +m 2 ) and the reduced relative momentum k = 1 2 |p 1 − p 2 |, where p 1 and p 2 are the momenta of the two particles evaluated in their center of mass (CM) frame of reference.
To obtain the SE solution we follow the standard procedure of separation of variables and partial waves decomposition (Eq. (A.1)) information on which, including the standard notation, can be found in any quantum mechanics textbook (e.g. [80]) For low-energy problems, as the one at hand, this series converges quite fast. The Legendre polynomials P l are easily computed using the GNU scientific library (GSL) [81], hence the only thing left to evaluate is the radial wave function R k,l (r ), which satisfies the radial Schrödinger equation (RSE) d 2 u k,l (r ) dr 2 = 2μV I,s,l, j (r ) where u k,l (r ) = r R k,l (r ) and V I,s,l, j (r ) is the interaction potential between the particle species of interest. As marked by the subindices the potential can depend on all quantum numbers describing the pair and consequently the solutions for the radial wave functions will differ between the different states. How to combine the different solutions to obtain the final correlation function will be explained at the end of this subsection, but first let us turn the focus to solving the Schrödinger equation for fixed values of k, I , s, l and j.
To avoid the overuse of indices, the radial wave functions u k,l (r ), R k,l (r ) and the potential V I,s,l, j (r ) will be written as u(r ), R(r ) and V (r ) respectively. In order to solve Eq. (A.2) numerically we employed a simple Euler-method equipped with an adaptive grid which automatically adjusts its step size when computing u(r ). The basic idea is to make the step smaller when computing regions of the wave function that have a strong non-linear behavior, i.e. a large second derivative. The discretized version of Eq. (A.2) using the Euler method and a dynamic step size reads where u i is the value of u(r ) at the i-th point on the discrete r -grid, Δ i = r i+1 − r i is the distance on the grid to the next evaluation point. The relation u i = F i u r is given by combining Eqs. (A.2) and (A.4), thus in order to adjust the step size one could simply keep the last term in Eq. (A.3) constant. This ensures that Δ 2 i ∝ 1/u i and the non-linear increase of u i+1 is kept constant throughout the computation. We have explored the possibility to select a more sophisticated numerical method, e.g. a Numerov method. However the inclusion of an adaptive grid tends to become quite complicated for higher-order numerical solvers. We have found out that the Euler method with an adaptive provides a faster solution, without compromising accuracy or precision, compared to the Numerov method without an adaptive grid. For this reason we have decided to employ the former.
One additional problem that arises when solving the Schrödinger equation is related to the boundary conditions. In particular Eq. (A.3) can be used to get the full wave function only if the first two points u 0,1 are known. Note that u(r ) = r R(r ), hence u(0) = u 0 = 0 can be used as the starting point of the solver. However the second point of the wave function depends on the potential and is thus a priori not known. Nevertheless if the second point u 1 is chosen randomly Eq. (A.3) will still yield the correct shape of the wave function at the expense of a wrong absolute normaliza- Fig. 13 An example for the CATS input when both the sand p-wave interactions are included. There are four distinct interaction channels: ψ 0 = ψ 1S0 + ψ 1P1 + ..., ψ 1 = ψ 3S1 + ψ 3P0 + ..., ψ 2 = ψ 3S1 + ψ 3P1 + ... and ψ 3 = ψ 3S1 + ψ 3P2 + .... The dots represent all partial waves beyond the p-wave, which are assumed to have no influence on C(k) and are thus considered irrelevant. The grayed out states are not allowed in the case of p-p correlation, but can be included if needed, e.g. for pcorrelations tion. To obtain the correct normalization a second boundary condition is needed. A straightforward approach to overcome this issue is to use the asymptotic solution. For a short range interaction potential V (r ), the evolution of the asymptotic wave function is governed by a phase shifted free or Coulomb wave. CATS is able to check when the solution has reached its asymptotic region and then simply match the numerical solution to the asymptotic one. This procedure allows to determine not only the normalization of u(r ) but the phase shift of the potential V (r ) as well (see Fig. 2). The total wave function ψ k (r) is given by the sum over all l partial waves Here the value of l max is determined by the condition for convergence, namely when l > l max the partial waves are practically zero in the asymptotic range. The sum is then split in two parts -first all numerically evaluated partial waves are summed up and than all remaining partial waves are added using their asymptotic solution obtained from GSL.
Another important thing to consider is the possibility to have different spin and isospin states. Moreover since the total angular momentum J is a good quantum number for potentials with spin-orbit terms, the degeneracy in J has to be taken into account as well. All of these contributions will result in different partial wave states 2S+1 L J and the user should provide CATS with all relevant information in order to get a correct total correlation function. The way this is achieved is by introducing the notation of interaction channels (see Fig. 13). An interaction channel is defined as the total wave function for a specific spin and isospin state and each partial wave that is included in the computation is considered to be in a fixed J state. The user should carefully examine all possible permutations to include each channel that contributes to the final correlation and should provide CATS with the probabilities that the particle pair can be found in particular channel.
We provide an example for the p-p interaction. For this system the isospin can only be 1, hence it is not considered in the calculation. However the total spin S of the system can be either 0 or 1. Since those states correspond to a singlet and a triplet configuration respectively, as long as there is no preferred polarization the probability to be in the singlet state is 1/4 and in the triplet 3/4. Those weights can be computed from the simple formula regarding spin degeneracy w (S) = (2S + 1) (2s 1 + 1)(2s 2 + 1) , where s 1 and s 2 are the spins of the individual particles. If one would consider only the s-waves in the computation the relevant partial waves are 1 S 0 and 3 S 1 . In fact for p-p, due to the quantum statistics of identical fermions, the 3 S 1 state is not allowed and therefore no potential will be considered in CATS. Nevertheless it is very important to note that the weights of 1/4 and 3/4 still need to be set, since even if the 3 S 1 partial wave is canceled out for p-p, the total wave function of the S = 1 is still included in C(k) with a flat contribution. The situation becomes a bit more complex if the p-waves are to be considered as well. In the S = 0 the 1 P 1 is the only p-partial wave and for p-p it is canceled out again due to the quantum statistics. However the S = 1 state has 3 possible p-waves, namely 3 P 0 , 3 P 1 and 3 P 2 . These states have the same S and L quantum numbers, but different J . Hence to compute their relative contributions one needs to consider the degeneracy in J which is given by w (J ) = (2J + 1) (2L + 1)(2S + 1) .
From the user's perspective, the inputs that CATS has to be provided with are the number of interaction channels, the number of partial waves to be included in each channel and finally the potentials for each partial wave state. In the example from Fig. 13 the user needs to define 4 channels each containing two partial waves. However, due to the Pauli blocking, the 1 P 1 and 3 S 1 states can be safely left undefined in the code.