S-shell $\Lambda\Lambda$ hypernuclei based on chiral interactions

We generalize the Jacobi no-core shell model (J-NCSM) to study double-strangeness hypernuclei. All particle conversions in the strangeness $S=-1,-2$ sectors are explicitly taken into account. In two-body space, such transitions may lead to the coupling between states of identical particles and of non-identical ones. Therefore, a careful consideration is required when determining the combinatorial factors that connect the many-body potential matrix elements and the free-space two-body potentials. Using second quantization, we systematically derive the combinatorial factors in question for $S=0,-1,-2$ sectors. As a first application, we use the J-NCSM to investigate $\Lambda \Lambda$ s-shell hypernuclei based on hyperon-hyperon (YY) potentials derived within chiral effective field theory at leading order (LO) and up to next-to-leading order (NLO). We find that the LO potential overbinds $^{\text{ }\text{ }\text{ } \text{}6}_{\Lambda \Lambda}\text{He}$ while the prediction of the NLO interaction is close to experiment. Both interactions also yield a bound state for $^{\text{ }\text{ }\text{ } \text{}5}_{\Lambda \Lambda}\text{He}$. The $^{\text{}\text{ }\text{ }\text{}4}_{\Lambda \Lambda}\text{H}$ system is predicted to be unbound.


Introduction
The scarcity of hyperon-nucleon (YN; Y=Λ, Σ) data and the almost complete lack of direct empirical information on the hyperon-hyperon (YY) and ΞN systems a e-mail: h.le@fz-juelich.de b e-mail: j.haidenbauer@fz-juelich.de c e-mail: meissner@hiskp.uni-bonn.de d e-mail: a.nogga@fz-juelich.de poses an enormous challenge for theorists in the attempt to derive baryon-baryon (BB) interactions in the strangeness sector on a microscopic level. By exploiting SU(3) flavor symmetry, several sophisticated YN potentials have been derived [1][2][3][4][5], which all describe the available YN data on an adequate quantitative level. The situation however remains largely unsatisfactory for the strangeness S = −2 sector, at least for the foreseeable future, because it is practically impossible to perform direct YY scattering experiments and so far there have been no two-body S = −2 bound states observed. Data on ΛΛ-and Ξ hypernuclei are therefore an indispensable source of information that can provide valuable additional constraints for constructing YY interactions. The latter requires solving the exact S = −2 many-body Hamiltonian with microscopic twoand higher-body BB interactions as input.
In the present work, we utilize the Jacobi no-core shell model (J-NCSM) [6,7] to study double-strangeness hypernuclei. Historically, since the first observations of ΛΛ hypernuclei, 10 ΛΛ Be [8], 6 ΛΛ He [9] and especially after publication of the so-called Nagara event [10,11], various approaches have been employed to study doublystrange hypernuclei [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. For example, Nemura et al. used the so-called stochastic variational method in combination with phenomenological effective central ΛN and ΛΛ potentials to investigate ΛΛ s-shell hypernuclei [12]. Thereby, it was assumed that the effects of tensor forces, three-body forces and Λ − Σ conversion are effectively included in such central potentials by fitting the binding energies of the s-shell (core) Λ-hypernuclei. In their later work [13], the channel coupling, e.g., ΛN − ΣN or ΛΛ−ΞN was fully taken into account, but again the YY interaction consisted of only central potentials. Hiyama and co-workers have successfully applied the Jacobian-coordinate Gaussian expansion method to ΛΛ arXiv:2103.08395v1 [nucl-th] 15 Mar 2021 hypernuclei with A = 6 − 11, which treats hypernuclear systems as three-, four-, or five-cluster structures [14,15,23]. The authors have advanced the approach in order to allow for all possible rearrangement channels so that any changes in the dynamic structure due to Λ interactions can be taken into account. The interactions between a Λ and a cluster are approximated using the simulated G-matrix YN potentials that are derived from a set of one-boson-exchange potentials. Here, the ΛN − ΣN and ΛΛ−ΞN couplings were not treated explicitly but the tuning parameters of the simulated potentials are chosen to reproduce some experimental separation energies such as those of 5 Λ He and 6 ΛΛ He. It is therefore rather difficult to relate the properties of the employed potentials to the free-space BB interactions. Filikhin and Gal have solved the Faddeev-Yakubovsky equations formulated for three-(ΛΛα), four-cluster (ΛΛαα) and ΛΛnp components [16][17][18]. Their calculations are also based on simulated potentials similar to those used in the works of Hiyama et al. but the ΛΛ interactions were mainly restricted to the s-wave. Lately, Contessi and co-workers [25] have combined the stochastic variational method with pionless effective field theory (EFT) interactions at LO to investigate the consistency of A = 4 − 6 ΛΛ hypernuclei.
Recently, the Jacobi NCSM has been successfully employed by us in studies of single-Λ hypernuclei up to A = 7 [27,28]. In these investigations, the full complexity of the underlying nucleon-nucleon (NN) and YN interactions (tensor forces, channel coupling) could be incorporated. Now we extend the Jacobi NCSM to S = −2 systems. Also here all channel couplings, i.e., ΛN − ΣN , and ΛΛ − ΛΣ − ΣΣ − ΞN are explicitly considered. As first application, we use the approach to obtain predictions of the chiral leading order (LO) [29] and next-to-leading order (NLO) [30,31] YY interactions for ΛΛ s-shell hypernuclei. Chiral EFT [32] is a very successful tool for describing the NN interaction (see [33] and references therein) and allows for accurate calculations of nuclear observables [34][35][36][37]. The YN interaction derived within the chiral EFT approach up to NLO likewise leads to realistic results for (s-and p-shell) hypernuclei [28,[38][39][40][41] and for nuclear matter [31,42]. It is therefore of great interest to study the predictions of the chiral YY potentials for ΛΛ hypernuclei.
The paper is structured in the following way. We present in Section 2 details of the technical realization. Some relevant formulas are also provided in the appendix. First results for the 6 ΛΛ He, 5 ΛΛ He and 4 ΛΛ H hypernuclei are discussed in Section 3. Finally, conclusions and an outlook are given in Section 4. In this section, we generalize our Jacobi no-core shell model (J-NCSM) formalism [28] to S = −2 hypernuclei. Adding a second Λ hyperon to single-strangeness systems complicates the numerical realization in many ways. All particle conversions that involve a Ξ hyperon, for instance ΛΛ ↔ ΞN, ΣΣ ↔ ΞN or ΛΣ ↔ ΞN change the total number of nucleons in the system by one. The latter must be explicitly taken into account for the many-body Hamiltonian and for the basis states. Furthermore, particle conversions in both S = −1 and S = −2 sectors can also lead to couplings between states of identical and non-identical hyperons. Because of that, special attention is required when evaluating the Hamiltonian matrix elements. These issues will thoroughly be addressed in this section and in Appendix A. We start with the construction of many-body basis states first. Since the total number of nucleons in the system can change depending on the strange particles, we split the basis functions into two orthogonal sets: one set that involves two singly strange hyperons referred to as |α * (Y1Y2) , and the other that contains a doubly strange Ξ hyperon denoted as |α * (Ξ) . The former is constructed by coupling the antisymmetrized states of A − 2 nucleons, |α (A−2)N , to the states describing a system of two hyperons, Here the inequality Y 1 ≤ Y 2 indicates the fact that we distinguish among the three two-hyperon states |ΛΛ , |ΛΣ and |ΣΣ but do not consider the |ΣΛ state explicitly. The notations in Eq. (1) are the same as introduced in Ref. [7,28]. For example, the symbol α (A−2)N stands for all quantum numbers characterizing the antisymmetrized states of A − 2 nucleons: the total number of oscillator quanta N A−2 , total angular momentum J A−2 , isospin T A−2 and state index ζ A−2 as well. Similarly, α Y1Y2 stands for a complete set of quantum numbers describing the subcluster of two hyperons Y 1 and Y 2 . The principal quantum number n λ of the harmonic oscillator (HO) together with the orbital angular λ describe the relative motion of the (A − 2)N core with respect to the centerof-mass (C.M.) of the Y 1 Y 2 subcluster. The orders, in which these quantum numbers are coupled, are shown after the semicolon. As for the transition coefficients of for standard nuclei and single Λ hypernuclei [7,28], the corresponding momenta or position vectors point to Y 1 and the A − 2 cluster, respectively. Analogously, in order to construct the basis |α * (Ξ) , one combines the antisymmetrized states of an (A−1)N system, |α (A−1)N , with the HO states, |Ξ , describing the relative motion of a Ξ hyperon with respect to the C.M. of the (A-1)N subcluster ( Here, also α (A−1)N denotes a set of quantum numbers representing an antisymmetrized state of A − 1 nucleons. The relative motion of a Ξ hyperon is labelled by the HO principal quantum number n Ξ , the orbital angular momentum l Ξ and spin s Ξ = 1 2 which combine together to form the total angular momentum I Ξ , and the isospin t Ξ = 1 2 . Again, following the definition of our coefficients of fractional parantage (CFPs) [7], the momentum or position vector points towards the spectator particel, i.e. the Ξ. Finally, the last lines in Eqs. (1,2) also show the graphical representations of the states.
With the basis states defined in Eq. (1,2), the S = −2 hypernuclear wave function |Ψ (πJT ) can be expanded as The expansion coefficients are obtained when diagonalizing the S = −2 Hamiltonian in the basis Eq. (1,2). For practical calculations, the model space is truncated by limiting the total HO energy quantum number N = Of course, by doing so, the computed binding energies will be N max -and HO ω-dependent. For extracting the converged results, we follow the two-step extrapolation procedure that has been extensively employed for nuclear and single-Λ hypernuclear calculations within the J-NCSM approach [7,28]. For energies, we first define E N for a given N m ax = N by minimizing the energy with respect to ω. Then we perform an exponential fit to E N to extrapolate to N → ∞.

S = −2 many-body Hamiltoninan
For the solution of the A-body Schrödinger equation, we use a standard, iterative Lanczos solver. In order to introduce the pertinent ingredients, we will in the following present the evaluation of an expectation value of H. The extension to the evalution of H |Ψ is then straightforward. Using the wave function in Eq. (3), one can write down the final expression for the energy expectation value as follows The last line in Eq. (5) is obtained by exploiting the hermiticity of the Hamiltonian. It should be clear from Eq. (5) that the part of the Hamiltonian that only involves the doubly-strange hyperon Ξ does not contribute to the matrix element α * (Y1Y2) | H | α * (Y1Y2) (in the first line). Likewise, α * (Ξ) |H | α * (Ξ) will not receive any contributions from the part of the Hamiltonian that contains two singly-strange hyperons Y 1 and Y 2 , whereas the last term is nonzero only for the transition potentials in the S = −2 channels. Therefore, in order to write down the explicit form of the S = −2 A-body Hamiltonian, we distinguish three parts of the Hamiltonian, namely H Y1Y2 , H Ξ and H S=−2 Y1Y2−ΞN , which contributes to the matrix elements in the first, second and third lines in Eq. (5), respectively. The first part of the Hamiltonian H Y1Y2 corresponds to a system consisting of A−2 nucleons and two singly-strange hyperons Y 1 and Y 2 , and has the following form, and m N are the Y 1 , Y 2 hyperon and nucleon rest masses, respectively. M (t Y1 , t Y2 ) denotes the total rest mass of while µ iY1 and µ Y1Y2 are the YN and YY reduced masses, respectively. The rest mass differences within the nucleonand hyperon-isospin multiplets are neglected.
are the nucleon-nucleon (NN), YN and YY potentials. Finally, the last term in Eq. (6) accounts for the difference in the rest masses of the hyperons arising due to particle conversions.
Likewise, the second Hamiltonian, H Ξ (involving a Ξ hyperon) corresponds to a system composed of a Ξ hyperon and A − 1 nucleons. Hence, where m Ξ is the Ξ hyperon rest mass and µ iΞ is the reduced mass of a Ξ and a nucleon. The total mass of the system is now given by is the ΞN potential. The ellipses in Eqs. (6,7) stand for those higher-body forces that are omitted here. The transition Hamiltonian H S=−2 Y1Y2,ΞN is simply given by the YY-ΞN transition potential The evaluation of the non-strange part, does not require any new transition coefficients, and can be performed analogously as done for the S = −1 systems [28]. Furthermore, the combinatorial factors that relate the A-body matrix elements α * (Y1Y2) | H S=0 Y1Y2 | α * (Y1Y2) and α * (Ξ) |H S=0 Ξ | α * (Ξ) to the two-nucleon matrix elements in the two-body sector are given by the binomial coefficients of A nucl 2 = A nucl (A nucl −1)/2 with A nucl = A−2 and A nucl = A−1, respectively, being the number of nucleons in the system (see Appendix A for the definition of the combinatorial factors).
The matrix elements of the double-strange part H S=−2 of the Hamiltonian, are evaluated analogously. Indeed, in order to calculate the last two terms in Eq. (11), one simply needs to expand the states |α * (Ξ) in the complete set of intermediate states |α * (ΞN ) that explicitly single out a ΞN pair, Here the transition coefficients α * (Ξ) |α * (ΞN ) can be computed using the expression Eq. (A.6) in Ref. [28]. It is easy to see that the last term in Eq. (11), , differs from the matrix element of the two-body ΞN Hamiltonian in the |ΞN basis by a combinatorial factor of A − 1. The factor that relates α * (Y1Y2) | H S=−2 Y1Y2−ΞN | α * (Ξ) to the two-body transition potential V Y1Y2−ΞN is, however, not obvious because of possible couplings between identical and non-identical two-body states, for instance, ΣΣ − ΞN or ΛΛ − ΞN . In Appendix A, we have shown that, in this case, the corresponding combinatorial factor is √ A − 1 (see Table 4).

Separation of a YN pair
Let us now discuss the evaluation of the second term in Eq. (9) that involves the singly-strange Hamiltonian H S=−1 of Eq. (6), in some details since it requires new sets of transition coefficients. Here, in order to compute the matrix elements α * (Y1Y2) | H S=−1 Y1Y2 | α * (Y1Y2) , one needs to employ other sets of intermediate states that explicitly separate out a YN pair . Obviously, each of the hyperons, Y 1 and Y 2 , can interact with a nucleon independently (as it is clearly seen from the expression for in Eq. (6)). It is then instructive to exploit two separate intermediate sets, namely | α * (Y1N ) * (Y2) and | α * (Y2N ) * (Y1) . The first set, | α * (Y1N ) * (Y2) , is needed when computing the matrix elements of the first two terms of H S=−1

Y1Y2
where Y 1 is the active hyperon while Y 2 plays the role of a spectator. Similarly, the second set, | α * (Y2N ) * (Y1) , is useful for evaluating the two remaining terms in Eq. (6) where the roles of Y 1 and Y 2 hyperons have been interchanged (i.e.,Y 2 is now the active particle). The construction of these bases is straightforward. For example, the first set can be formed by combining the hyperon states |Y 2 , depending on the Jacobi coordinate of the Y 2 hyperon relative to the C.M. of the ((A − 3)N + Y 1 N ) subcluster, with the |α * (Y1N ) states constructed in Eq. (9) in [28]. Thus, and, similarly In both of these basis states, we have one momentum/position of the spectator pointing towards the spectator, the one of the pair pointing towards the hyperon and the third momentum/position pointing towards the A−3 cluster. Clearly, each of the above two auxiliary sets is complete with respect to the basis states |α * (Y1Y2) in Eq. (1). This in turn allows for the following expansions or, Obviously, when Y 1 and Y 2 are identical, the two auxiliary sets Eqs. (14,15) are the same, and there is no need to distinguish between the two expansions. In any case, the expansion coefficients in Eqs. (16,17) are very similar to each other and can be computed analogously. In the following, we focus on the transition coefficients of the first expansion. For computing the overlap, (α * (Y1N ) ) * (Y2) |α * (Y1Y2) , we make use of another set of auxiliary states, |(α * (Y1) ) * (Y2) that explicitly single out the Y 1 and Y 2 hyperons. These states are obtained by coupling the hyperon states |Y 2 to the basis states of the ((A − 2)N + Y 1 ) system, |α * (Y1) A−1 , defined in Eq. (4) in [28], The third line in Eq. (18) is to illustrate how the quantum number of the three subclusters: (A-2) nucleons, Y 1 and Y 2 hyperons, are combined to form the intermediate states with the definite quantum numbers N , J and T . Exploiting the completeness of the auxiliary states | α * (Y1) * (Y2) , the transition coefficient in Eq. (16) then becomes where a summation over the states | α * (Y1) * (Y2) is im- (19) is essentially given by the transition coefficients of a system consisting of (A − 2) nucleons and the Y 1 hyperon (see Eq. (11) in [28]), whereas the second term α * (Y1) * (Y2) |α * (Y1Y2) can quickly be deduced from Eq. (11) in [7], with, .
The transition coefficients for the second expansion in Eq. (17) are computed analogously. Taking into account the expansions Eqs. (16,17), the matrix element The subscript in each term on the right-hand side of Eq. (21) specifies the hyperon spectator. The first contribution is further given by The expression for the second term in Eq. (21) is obtained from Eq. (22) by interchanging the roles of the Y 1 and Y 2 hyperons in the intermediate states, Although Eqs. (22,23) are very similar to the expression for computing the Hamiltonian matrix elements in S = −1 systems, the presence of a hyperon spectator Y 2 (Y 1 ) makes it rather difficult to determine the proper combinatorial factors that relate the many-body matrix ele- to the YN Hamiltonian matrix elements in the two-body sector. These factors are also provided in Table 3 in Appendix A. From Table  3, one can clearly see that the corresponding factors depend not only on the total number of nucleons but also on the two hyperons Y 1 and Y 2 in the intermediate states.

Results
In this section, as a first application, we report results for the ΛΛ s-shell hypernuclei 2 ), and 6 ΛΛ He(0 + , 0). To zeroth approximation, these systems can be regarded as a ΛΛ pair in the 1 S 0 state being attached to the corresponding corenuclei predominantly in their ground states. While the quantum numbers of 5 ΛΛ He, (J + , T ) = ( 1 2 + , 1 2 ), are obvious, those for the 4 ΛΛ H hypernucleus are chosen according to our observations that the state with (J + , T ) = (1 + , 0) is the lowest-lying level and in many calculations the one closest to binding of all A = 4 S = −2 hypernuclei. Therefore, we will report our results for this state below.
For all calculations presented here, we employ BB interactions that are derived within chiral EFT [32]. The high-order semilocal momentum-space regularized potential with a regulator of Λ N = 450 MeV (N 4 LO+(450)) [33], SRG-evolved to λ N N = 1.6 fm -1 , is adopted for describing the NN interaction. The next-toleading order potential NLO19 [4] with a chiral cutoff of Λ Y = 650 MeV and an SRG parameter of λ Y N = 0.868 fm -1 is used for the YN interaction. We remark that the chosen NN and YN potentials successfully predict the empirical Λ-separation energies for 3 Λ H, 4 Λ He(1 + ) and 5 Λ He but slightly underbind 4 Λ He(0 + ) [28]. To describe the two-body interactions in the S = −2 sector, we utilize the chiral YY interactions at LO [29] and up to NLO [30,31], with a chiral cutoff of Λ Y Y = 600 MeV.
One of our primary aims here is to establish the predictions of these chiral YY potentials for double-Λ s-shell hypernuclei. Ultimately, it is expected that results from such a study may provide useful additional constraints for constructing realistic S = −2 BB interaction potentials, given the scarcity of direct empirical information on the underlying two-body systems (ΛΛ, ΞN , ...). Due to the latter circumstance, in the chiral approach (as well as in meson-exchange and/or constituent quark models) the assumption of SU(3) f symmetry is an essential prerequisite for deriving pertinent potentials. For example, in chiral EFT the shortdistance dynamics is represented by contact terms which involve low-energy constants (LECs) that need to be determined from a fit to data [32]. SU(3) symmetry strongly limits the number of independent LECs [3]. However, at NLO, there are two LECs which are only present in the S = −2 sector, and which contribute to the interaction in the spin-and isospin zero channel, specifically to the 1 S 0 partial wave of ΛΛ. They correspond to the SU(3) singlet irreducible representation, see Ref. [30], and are denoted byC 1 and C 1 , respectively, in that work. These have been fixed by considering the extremely sparse and uncertain YY data (i.e., a total cross section for Ξ − p − ΛΛ [43] and the upper limits of elastic and inelastic Ξ − p cross sections [44]). Clearly, such poor empirical data do not allow for a reliable quantitative determination of the unknown strength of the two contact terms in question. Nevertheless, it turned out that reasonable choices for the C 1 's can be made [30,31] and the YY cross sections predicted by the two NLO potentials are fairly consistent with the experiments. Furthermore, the ΛΛ 1 S 0 scattering lengths predicted by these interactions are compatible with values inferred from empirical information [45,46]. The LO interaction yields a somewhat large scattering length in comparison to those values and it also exhibits a rather strong regulator dependence [29].
It should be pointed out that our initial NLO interaction for S = −2 [30] and the updated version [31] differ only in the antisymmetric SU(3) f component which means essentially only in the strength of the ΞN interaction in the 3 S 1 partial wave. This has an impact on the corresponding in-medium properties of the Ξ. Specifically, the updated version from 2019 [31] yields a moderately attractive Ξ single-particle potential that is roughly in line [47] with recent experimental evidence that the existence of bound Ξ-hypernuclei is very likely [48]. With regard to ΛΛ systems, we observe that the two realizations yield very similar binding energies for the double-Λ s-shell hypernuclei. This indicates that, in general, the actual strength of the spin-triplet ΞN interaction has little influence on few-body observables related to ΛΛ. In the following, we therefore present results for the LO and the updated NLO interactions for a chiral cutoff of Λ Y Y = 600 MeV. In order to speed up the convergence, both YY potentials are also SRGevolved. We use a wide range of the SRG flow parameters, namely 1.4 ≤ λ Y Y ≤ 3.0 fm -1 , to quantify the contribution of possible SRG-induced YYN three-body forces.
3.1 6 ΛΛ He(0 + , 0) The 6 ΛΛ He hypernucleus is so far the lightest double-Λ system being unambiguously established. Since the observation of the Nagara event [10], its ΛΛ separation energy, defined as B ΛΛ ( 6 ΛΛ He) = E( 4 He) − E( 6 ΛΛ He), has been exploited as a crucial constraint for constructing effective potentials that are then employed in manybody calculations like the Gaussian expansion method [13,49] or the cluster Faddeev-Yakubovsky approach [16,17]. The re-analysis of the Nagara event using the updated Ξ mass yielded a slightly smaller ΛΛ separation energy, B ΛΛ ( 6 ΛΛ He) = 6.91 ± 0.16 MeV [11,50], as compared to the initially estimated value of B ΛΛ ( 6 ΛΛ He) = 7.25 ± 0.19 [10]. This, in turn, may have direct consequences for theoretical predictions for potentially observable bound states of other s-shell ΛΛ hypernuclei, particularly the A = 4 double-Λ system [25,51], see also the discussion in 3.3. We note that the information about B ΛΛ ( 6 ΛΛ He) has not been directly utilized in order to constrain the LECs appearing in the chiral LO and NLO potentials. It is therefore of enormous interest to explore this double-Λ system using the two chiral interactions to scrutinize their consistency with the measured ΛΛ separation energy.
As mentioned earlier, in order to eliminate the effect of the finite-basis truncation on the binding energies, we follow the two-step extrapolation procedure as explained in [28]. The ω-and N -space extrapolations for E( 6 ΛΛ He) are illustrated in panels (a) and (b) of Fig. 1, respectively. Here, for illustration purposes, we present results for the NLO potential with λ Y Y = 1.8 fm -1 but stress that the convergence trend is similar for all other values of λ Y Y , and for the LO interaction. Also, the behavior of E( 6 ΛΛ He) with respect to ω and N resembles that of the binding energy of the parent hypernucleus 5 Λ He [28]. Furthermore, panel (b) clearly demonstrates a nice convergence pattern of the binding energy E( 6 ΛΛ He) computed for model spaces up to N max = 14. Likewise, the ΛΛ-separation energy B ΛΛ ( 6 ΛΛ He), displayed in panel (c), is also wellconverged for N max = 14 (practically with the same speed as that of E( 6 ΛΛ He)). Note that, for single-Λ hypernuclei, the separation energy B Λ converges somewhat faster than the individual binding energies. For S = −2 systems, we are also interested in the so-called ΛΛ excess binding energy which provides information about the strength of the ΛΛ interaction.B andĒ are spin averaged Λ-separation and binding energies of the hypernuclear core if the core supports several spin states. Cleary, this difference is also affected by the spin-dependent part of the Λcore interaction, dynamical changes in the core-nucleus structure as well as the mass-polarization effect [8,15]. For 6 ΛΛ He, the spin-dependent part of the Λ-core interaction vanishes because of the spin zero the parent nucleus 4 He, hence the difference will reflect the net contributions of the ΛΛ interactions and the 4 He core-distortion 1 (polarization) effects. In panel (d), we exemplify the model-space extrapolation for ∆B ΛΛ ( 6 ΛΛ He). Interestingly, ∆B ΛΛ converges with respect to N visibly faster than both the ΛΛ-separation and the binding energies.
Being able to accurately extract B ΛΛ ( 6 ΛΛ He) and ∆B ΛΛ ( 6 ΛΛ He), we are in a position to study the impact of the two chiral interactions on these quantities. The converged results for B ΛΛ and ∆B ΛΛ , calculated for a wide range of the SRG flow parameter λ Y Y , are presented in the left and right plots of Fig. 2, respectively. Evidently, the LO YY potential (blue triangles) produces too much attraction (more than 2 MeV as can be seen in the right panel), which, as a consequence, leads to overbinding by about 1.5 MeV in 6 ΛΛ He as can be seen in the left panel. On the other hand, the moderately attractive NLO interaction predicts a ΛΛ excess energy of ∆B ΛΛ ≈ 1.1 MeV, that is only slightly larger than the empirical value of ∆B exp ΛΛ = 0.67 ± 0.17 MeV [11,50]. For completeness, let us mention that the 1 Our preliminary results for the RMS distances of an NN pair and point-nucleon radii in 6 ΛΛ He, 5 Λ He and 4 He are very similar to each other which implies that the distortions of the 4 He core are rather small. However, we also note that Hiyama et al. in their study for A = 7 − 10 double-strangeness systems using the Gaussian-basis coupled cluster method found that the dynamical changes in the nuclear core structures are quite visible [15]. Further studies are necessary in order to clarify the discrepancy. 12 16 20 24 28 [MeV]    (d) ∆B ΛΛ ( 6 ΛΛ He) as a function of N . Fig. 1 Binding energy E, ΛΛ-separation energy B ΛΛ and separation-energy difference ∆B ΛΛ for 6 ΛΛ He computed using the YY NLO(600) interaction [31] that is SRG evolved to a flow parameter of λ Y Y = 1.8 fm -1 . The SMS N 4 LO+(450) potential [33] with λ N N = 1.6 fm -1 and the NLO19(650) potential [4] with λ Y N = 0.868 fm -1 are employed for the NN and YN interactions, respectively.   [11]. Same NN and YN interactions as in Fig. 1 [29]) and a = −0.66 fm (NLO [30]), respectively. It is rather remarkable that both, B ΛΛ ( 6 ΛΛ He) and ∆B ΛΛ ( 6 ΛΛ He), exhibit a rather weak dependence on the SRG YY parameter λ Y Y . With an order of 100 keV, it is at least one order of magnitude smaller than the variation of, say, B Λ ( 5 Λ He) with respect to the SRG YN flow parameter λ Y N [28]. The insensitivity of the ΛΛseparation energy to the SRG evolution indicates that the SRG-induced YYN forces are negligibly small. This is probably the result of a rather weak ΛΛ interaction.
Finally, we benchmark the probabilities of finding one Σ (P ΛΣ ) or two Σ (P ΣΣ ), or the Ξ hyperon (P Ξ ) in the ground-state wave function of 6 ΛΛ He obtained for the two chiral potentials. Such probabilities are summarized in Table 1 for several values of λ Y Y . Overall, the P ΛΣ and P ΣΣ probabilities are fairly small, but almost stable with respect to the SRG evolution of the YY interaction. Also, their dependence on the two considered potentials is practically negligible. We remark that the probability of finding a Σ in 5 Λ He for the employed NN and YN interactions is also very small, P Σ ( 5 Λ He) = 0.07%. In contrast, P Ξ is more sensitive to the evolution and also strongly influenced by the interactions. Surprisingly, the updated NLO potential, that yields a more attractive Ξ-nuclear interaction [31], predicts a considerably smaller Ξ probability (less than 0.2 % for λ Y Y = 3.0 fm -1 ) as compared to the value of P Ξ = 1.1% obtained for the LO at the same λ Y Y . This reflects our observation in the S = −1 sector that there is no simple one-to-one connection between the probabilities of finding a hyperon particle (Σ, Ξ) and the interaction strength.

5
ΛΛ He( 1 2 + , 1 2 ) The next system that we investigate is 5 ΛΛ He. Although the existence of 5 ΛΛ He has not been experimentally confirmed yet, most of the many-body calculations employ-ing effective potentials that reproduce the separation energy B ΛΛ ( 6 ΛΛ He) predict a particle-stable bound state of 5 ΛΛ He [13,16,51]. However, there are visible discrepancies among the values of B ΛΛ ( 5 ΛΛ He) predicted by different numerical approaches or different interaction models. Additionally, it has been observed in Faddeev cluster calculations that there is an almost linear correlation between the calculated values of B ΛΛ for the 5 ΛΛ He ( 5 ΛΛ H) and 6 ΛΛ He hypernuclei [16]. Such a behavior was also seen in the study based on pionless EFT [25]. It will be of interest to see whether one observes a similar correlation using other realizations of the chiral interactions. However, at this stage, we postpone that question to a future investigation and focus on the different effects of the LO and NLO potentials on B ΛΛ ( 5 ΛΛ He) instead. The ω-and N -extrapolation of the binding energy E, ΛΛ-separation energy B ΛΛ and the separation-energy difference ∆B ΛΛ of 5 ΛΛ He are illustrated in Fig. 3. Here, the results are shown for the NLO potential with a flow parameter of λ Y Y = 1.8 fm -1 and for model spaces up to N max = 16. Note that in the case of 4 Λ He, the energy calculations were performed for model spaces up to N max ( 4 Λ He) = 22 in order to achieve a good convergence. Calculations with such large model spaces are currently not feasible for 5 ΛΛ He because of computermemory constraints. Nonetheless, the illustrative results in Fig. 3 clearly indicate that well-converged results are achieved for this double-Λ hypernucleus already for model spaces up to N max = 16. Moreover, the employed two-step extrapolation procedure also allows for a reliable estimate of the truncation uncertainty. Let us further remark that, when calculating the excess energy we do not simply assign the ground-state Λ-separation energy B Λ ( 4 Λ He, 0 + ) to B Λ ( 4 Λ He) but rather use a spin-averaged value B Λ ( 4 Λ He) of the ground-state doublet [15] By doing so, the computed quantity ∆B ΛΛ ( 5 ΛΛ He) will be less dependent on the spin-dependence effect of the Λ-core interactions, and, therefore, can be used as a measure of the ΛΛ interaction strength, provided that the nuclear contraction effects are small. The results for B ΛΛ ( 5 ΛΛ He) and ∆B ΛΛ ( 5 ΛΛ He) calculated for the two interactions and a wide range of flow parameter, 1.4 ≤ λ Y Y ≤ 3.0 fm -1 are shown in Fig. 4. Overall, we observe a very weak dependence of these two quantities on the SRG flow parameter, like for 6 ΛΛ He, reinforcing the insignificance of SRG-induced YYN forces.     Binding energy E, ΛΛ-separation energy B ΛΛ and separation-energy difference ∆B ΛΛ for 5 ΛΛ He computed using the YY NLO(600) interaction that is SRG evolved to a flow parameter of λ Y Y = 1.8 fm -1 . Same NN and YN interactions as in Fig. 1     Again, the LO interaction predicts a much larger ΛΛseparation energy and a more significant ΛΛ interaction strength than the one at NLO. In either case, the ΛΛ excess energy ∆B ΛΛ computed for 5 ΛΛ He, slightly exceeds the corresponding one for 6 ΛΛ He, by about 0.23 and 0.5 MeV for the LO and NLO interactions, respectively. The main deviations should come from the nuclear-core distortion and the suppression of the ΛΛ − ΞN coupling in 6 ΛΛ He as discussed in [18,52,53]. However, it is necessary to carefully study the impact of the employed interactions on the results before a final conclusion can be drawn. We further note that Filikhin and Gal [16] in their Faddeev cluster calculations, based on potentials that simulate the low-energy s-wave scattering parameters of some Nijmegen interaction models, obtained an opposite relation, namely ∆B ΛΛ ( 5 ΛΛ He) < ∆B ΛΛ ( 6 ΛΛ He). As a consequence, our results do also not fit into the correlation of ∆B ΛΛ ( 5 ΛΛ He) and ∆B ΛΛ ( 6 ΛΛ He) shown in the same work. We will need to study more interactions in future to understand whether such a correlation can also be established using chiral interactions.
It is also very interesting to point out that the ΛΛseparation energies B ΛΛ for both 5 ΛΛ He and 6 ΛΛ He predicted by the NLO potential are surprisingly close to the results obtained by Nemura et al., B ΛΛ ( 5 ΛΛ He) = 3.66 MeV, B ΛΛ ( 6 ΛΛ He) = 7.54 MeV, using the modified Nijmegen YY potential (mND s ) [13]. Finally, we provide in Table 2 the probabilities of finding a Σ (P ΛΣ ), double Σ (P ΣΣ ), or a Ξ (P Ξ ) in the 5 ΛΛ He ground-state wave function, computed with the two potentials and several SRG values, λ Y Y = 1.4, 2.0 and 3.0 fm -1 . Apparently, all the probabilities including also P Ξ exhibit a rather weak sensitivity to the flow parameter λ Y Y . The two interactions seem to have little impact on the Σ-probabilities (P ΛΣ and P ΣΣ ) but strongly influence P Ξ . Like in the 6 ΛΛ He system, here, the LO potential yields considerably larger Ξ-probabilities as compared to the values predicted by the NLO interaction. It also clearly sticks out from Tables 1 and 2 that the probabilities of finding a Σ or Ξ hyperon in 5 ΛΛ He are visibly larger than the corresponding ones in 6 ΛΛ He. This is indeed consistent with the Σ-probabilities in the groundstate wave functions of their parent hypernuclei (e.g., P Σ ( 4 Λ He) = 0.43 % and P Σ ( 5 Λ He) = 0.07 %), and more importantly, is consistent with the suppression of particle conversions such as ΛΛ − ΞN in p-shell hypernuclei [52].

4
ΛΛ H(1 + , 0) Our final exploratory s-shell hypernucleus is 4 ΛΛ H. This system has been the subject of many theoretical and experimental studies. It turned out that theoretical predictions of the stability of 4 ΛΛ H against the 3 Λ H + Λ breakup are very sensitive to the interpretations of double-strangeness hypernuclear data, in particular, the 6 ΛΛ He hypernucleus [51]. Indeed, Nemura et al. [13] observed a particle-stable but loosely bound state of 4 ΛΛ H (just only about 2 keV below the 3 Λ H + Λ threshold for the mND s potential) using the fully coupled-channel stochastic variational method in combination with effective YY potentials that are fitted to reproduce the initially extracted value of B ΛΛ ( 6 ΛΛ He) = 7.25 ± 0.19 MeV [10]. The study by Filikhin and Gal [17] indicated, however, that there is a sizable model dependence. The authors found no bound state within an exact four-body (Faddeev-Yakubovsky) calculation for the ΛΛpn system, but a particle-stable 4 ΛΛ H hypernucleus when solving the (three-body) Faddeev equation for the ΛΛd cluster system. A more recent calculation by Contessi et al. [25], based on the pionless EFT interaction at LO, showed that the existence of a bound state in 4 ΛΛ H is not compatible with the corrected value of B ΛΛ ( 6 ΛΛ He) = 6.91 ± 0.16 MeV. Although the observation of 4 ΛΛ H was reported in an experiment at BNL [54], it has been recently invalidated by a thorough re-examination of the recorded events [55]. Nevertheless, the existence of a stable 4 ΛΛ H hypernucleus cannot be completely ruled out and the search for its experimental confirmation or exclusion is still ongoing.
In view of the previous calculations, it is interesting to see whether the chiral YY potential at NLO, that predicts similar results for A = 5 − 6 ΛΛ hypernuclei as the mND s interaction [13], also results in a loosely bound state for 4 ΛΛ H. It is well-known that NCSM calculations for very loosely bound systems like the hypertriton converge very slowly. Hence, in order to unambiguously answer that question, converged results for the binding energy of the parent 3 Λ H and the groundstate energy of 4 ΛΛ H are crucial. In panels (a) and (b) of Fig. 5 [41]. We conclude that a model space truncation of N max = 32 for the energy calculations in 4 ΛΛ H should be sufficient in order to draw conclusions about the stability of the system against Λ emission.
The extrapolated ground-state energies E( 4 ΛΛ H) for the NLO (red circles) and LO (blue triangles) potentials evolved to a wide range of flow parameters are displayed in panel (c) of Fig. 5. Here, the dashed black line together with the grey band represent the computed E( 3 Λ H) and the estimated uncertainty. Calcula-tions with the NLO potential seem to converge more slowly than the ones for the LO interaction. The NLO potential clearly leads to an unbound 4 ΛΛ H hypernucleus. Although our results for A = 5 and 6 are similar to the ones of Ref. [13], our results for A = 4 do not support the existence of a bound 4 ΛΛ H state. The LO results for 4 ΛΛ H likely hint at a particle-unstable system with respect to the hypertriton 3 Λ H. Admittedly, in order to draw a definite conclusion on the actual situation, the uncertainties of the calculation would have to be reduced. However, since the LO interaction considerably overbinds 6 ΛΛ He, very likely it overpredicts the actual attraction in the A = 4 system, too. Interestingly, in pionless EFT [25] a ΛΛ scattering length practically identical to that of our LO interaction was found as limit for which the 4 ΛΛ H system becomes bound.

Conclusions and outlook
In this work, we have generalized the J-NCSM formalism in order to include strangeness S = −2 hyperons. Using the second quantization approach, we systematically derived the necessary combinatorial factors that relate the Hamiltonian matrix elements in a many-body basis to the corresponding ones in a two-body basis for the S = 0, −1 and −2 sectors. A generalization to higher-strangeness sectors will be straightforward. We then applied the J-NCSM approach to compute predictions of the chiral YY interactions at LO and NLO for ΛΛ s-shell hypernuclei. To speed up the convergence, the two interactions are also evolved via SRG. Unlike for the S = −1 systems, here, we observed a very small effect of the SRG YY evolution on the ΛΛ-separation energies, implying negligible contributions of SRG-induced YYN forces. Furthermore, we found that the binding energy for 6 ΛΛ He predicted by the YY NLO potential is close to the empirical value while the LO interaction overbinds the system. Both interactions also yield a particle-stable 5 ΛΛ He hypernucleus, whereas 4 ΛΛ H is found to be unstable against a breakup to 3 Λ H + Λ. However, for a final conclusion, a more elaborate study that involves a more careful estimate of uncertainties stemming from various NN, YN and YY interactions is definitely necessary. It will be also very interesting to study the predictions of the chiral YY interactions for other s-shell ΛΛ systems such as 4 ΛΛ n or 4 ΛΛ He, as well as for p-shell hypernuclei. Finally, investigating possible Tjon-line like correlations for B ΛΛ of different systems is also of importance.

Appendix A: Many-body Schrödinger equation in second quantization
Generally, baryon-baryon (BB) interactions in the S = −2 sector can lead to couplings between states with identical particles and with non-identical particles, for example ΣΣ → N Ξ. Such transitions make it not straightforward to properly determine the combinatorial factors of free-space two-body potentials that are embedded in the A-body Hamiltonian matrix elements. In this appendix, we demonstrate that these factors can systematically be deduced by comparing the Schrödinger equation for A-body systems with the free-space twobody Schrödinger equation, provided that these equations are derived in a consistent way. We show explicit examples for systems of two and three particles, and then generalize to the A-baryon problems. We note that Glöckle and Miyagawa [56] have also derived a system of coupled-Faddeev equations for three-baryon systems taking into account full particle conversions. However it is not clear to us how to read off the involved combinatorial factors based on their equations. The authors of Ref. [57] have formulated the problem taking all permutations of particles explictly into account. This is however not consistent with the approach of BB interactions used in [29,30,58]. For directly taking these interactions into account, we therefore require to derive the combinatorial factors consistent with these interactions.
To derive the general Schrödinger equation, we will work with second quantization. The many-body Hamiltonian then has the form, where k i stands for a set of quantum numbers characterizing the particle state, i.e., momentum, spin, isospin as well as particle species λ i (N, Λ, Σ or Ξ). When it is necessary to separate the particle species λ i from other quantum numbers, we use k i = λ iki . Let us further assume that the potential matrix elements V k 1 k 2 ,k1k2 in Eq. (A.1) are antisymmetric under exchanges of two indices, i.e., V k 1 k 2 ,k1k2 = −V k 1 k 2 ,k2k1 = −V k 2 k 1 ,k1k2 = V k 2 k 1 ,k2k1 . Note that, there is no ordering imposed for quantum numbers of the incoming particles k 1 and k 2 or of the outgoing pair k 1 and k 2 in Eq. (A.1).
Appendix A.1: Two-body Schrödinger equation We start with the derivation of the Schrödinger equation in a two-particle basis. For that, we define the ordered two-body antisymmetrized basis states as with the right-hand side being the states in first quantization. Here, p 1 and p 2 also stand for the sets of quantum numbers (momentum, spin, isospin and particle species) describing particles 1 and 2, respectively. The completeness relation of the basis Eq. (A.2) for bases with particle species where the inequality p 1 < p 2 accounts for the ordering of the states in Eq. (A.2) where the leading sorting key is assumed to be particle species. Note that by exploiting the antisymmetry of the basis functions, the left hand side of Eq. (A.3) is equivalent to Hence, the summation over the ordered particle species on the left hand side of Eq. (A.3) can be replaced by a normal summation over all particle species but with a factor of 1 2 . For the case of two identical particles, i.e., λ 1 = λ 2 , the completeness relation becomes following similar lines. The factor 1 2 can also be absorbed into the definition of the states when one rewrites Eq. (A.5) as follows Now, exploiting the anticommutator relation for the creation and annihilation operators, the kinetic and potential matrix elements in the basis Eq. (A.2) are easily obtained In . (A.9) Here, it will be sufficient to consider only those components of Ψ (p 1 p 2 ) with p 1 < p 2 . Since the basis states are antisymmetric, the other components of Ψ (p 1 p 2 ) with p 1 > p 2 will differ from the ones with p 1 < p 2 by a simple phase factor. Plugging Eq. (A.7) into Eq. (A.9) and using p 1 < p 2 , one arrives at a general two-body Schrödinger equation We note that there is a factor of 2 in front of the potential matrix elements, which drops out for the case of the two-identical particle basis, i.e., λ 1 = λ 2 . In that case, we use p1<p2 → 1/2 p1,p2 and equation Eq. (A.10) becomes To better understand the prefactors of the potential matrix elements present in Eqs. (A.10,A.11), let us consider some explicit bases. In the first example, the basis consists of two two-particle states, one with identical particles and one with distinguishable particles, e.g., |{ΛΛ} and |{N Ξ} . Then, the completeness relation is obtained by combining Eqs. (A.3) and (A.6) leading to the following expression for the norm of the wave function (A.13) Therefore, we absorb the 1 √ 2 -factor into the amplitude of states by introducing a new set of the wave-function components, so that the Schrödinger equation Eqs. (A.10,A.11) for the two newly defined components possesses a symmetric form  where, for readability, we have omitted the dependence onp andp . Similarly, for the case where the basis consists of four states {|ΛΛ} , |{ΣΣ} , |{ΛΣ} and |{ΞN } , one analogously defines a new set of wavefunction components for which the Schrödinger equation again possesses a symmetric form One sees that there is a √ 2-factor for the transition between states of identical and of distinguishable particles, and a factor of 2 for the transition between states of nonidentical particles. It is important to mention that these factors are already included in the definition of the two-body potentials derived from chiral EFT [29,30] or phenomenological models [58] (see, e.g., Eq. (2) of [30]). We therefore denote these initial two-body potentials V λ1λ2,λ 1 λ 2 with an appropriate factor of √ 2 or 2 or 1 to be our new potentialṼ λ1λ2,λ 1 λ 2 . Expressing in terms of the new potentialsṼ , the Hamiltonian Eq. (A.18) now has a more intuitive form In the next step, we are going to derive a similar Schrödinger equation in a three-body basis. Then, by comparing the obtained equation with the one for twobody basis, we will be able to determine the corresponding combinatorial factors for the potentials in each strangeness sector.
Appendix A.2: Three-body Schrödinger equation We define the ordered three-body basis states in second quantization and its completeness relations as (A.20) The kinetic and potential matrix elements in the basis Eq. (A.20) read and The contributions from the potential operator are a little bit more cumbersome, but can be reduced to a compact form by exploiting the antisymmetry properties under the exchange of two indices of the potential as well as of the wave function. For example, the first three terms in Eq. (A. 22) give Analogously, the next three terms in Eq. (A.22) yield and, the three remaining terms result in can be written as which, as one expects, differs from the Schrödinger equation in the two-body basis Eq. (A.10) by the kinetic energy of the third particle and the two-body interactions between particles 1-3 and 2-3. Again, the factor of 2 in front of the potential vanishes when the incoming particles are identical and the summations include all statesp 1 ,p 2 etc. For illustration purposes, let us consider Eq. (A.28) in an explicit basis consisting of four states, |{N ΛΛ} , |{N ΣΣ} , |{N ΛΣ} and |{N N Ξ} . The norm of the wave function in this four-particlestate basis can be calculated as follows Based on Eq. (A.29), we define a new set of wavefunction components The Schrödinger equation Eq. (A.28), applying the wave function components in Eq. (A.30), now has a symmetric form, with T being a diagonal matrix In the last step, we have expressed the potential matrix elements in terms ofṼ as given in Eq. (A.19). Eqs. (A.31-A.33) define the combinatorial factors of the two-body potentials present in the three-body Hamiltonian. In following, we want to generalize this result to an A-body system.

Appendix A.3: A-body Schrödinger equation
With the preparation of the A = 3 system, we are now able to generalize the combinatorial factors to arbitrary A. For the kinetic energy, the generalization is trivial and leads to the sum of the single particle kinetic energies since no particle conversion can take place for this operator. Interactions are more involved. To the general A-body matrix element {p 1 . . . p A }|V |{p 1 . . . p A } of the n-particle interaction   different permutations of V k 1 ...k n ,k1...kn contribute. Therein, the first 1 n! is just from the definition of V . Following the same steps that lead to Eq. (A.25), these terms can be rearranged such that the application to an arbitrary state Ψ can be written as For this form, we assume that Ψ is represented using the ordered states p 1 < . . . < p A . Then only one of the (A − n)! different spectator permutations contributes. One of the A n terms is needed to make the sorting on the spectator particles and on the interacting particles independent from each other as done in Eq. (A.25). The other one is explicitly taken care of by the sum over i 1 < i 2 . . . < i n .
If the interacting particles are (partly) identical, we will again replace the sum over p i1 < . . . < p in by (partly) full sums and add the appropriate combinatorial factor (e.g. 1/2! in the case of two identical particles). Note that this factors depend on the kind of particles in the incoming Ψ state.
We again introduce rescaled wave functions by studying the norm of the states similar to Eq. (A.29). The appropriate factors for states with p particles species and n 1 , . . . , n p particles of each species are n 1 ! . . . n p !. The potential matrix element needs to be multiplied with (divided by) this factor for incoming states (outgoing states) to reexpress Eq. (A.36) in terms of Φ states. We note that the potential matrix in terms of these states is symmetrical. In summary, the potential matrix elements then reads n! n 1 ! . . . n p ! n 1 ! . . . n p ! V p i 1 ,...,p in ,pi 1 ,...,pi n . (A.37) Note that here the factor does not include the additional factor required when identical particles are involved in the sum of Eq. (A.36).
We then simplify the expressions by identifying nparticles that contribute identically to Eq. (A.36). The sum over i 1 < i 2 . . . < i n can then be reduced and tuples of outgoing states involving the same kind of particles can be combined by the appropriate factor. Finally, we build the ratio of the factors for the A-body and n-body systems to find the correct combinatorial factors that enter our J-NCSM calculations.
The result for the two matrix elements is