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

Abstract

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 $\varLambda \varLambda$ 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 $^{\,\,\,{\,}6}_{\varLambda \varLambda }\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}_{\varLambda \varLambda }\text {He}$. The $^{\text {}\text { }\text { }\text {}4}_{\varLambda \varLambda }\text {H}$ system is predicted to be unbound.

1 Introduction

The scarcity of hyperon-nucleon (YN; Y=$\varLambda$, $\varSigma$) data and the almost complete lack of direct empirical information on the hyperon-hyperon (YY) and $\varXi N$ systems 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 $\varLambda \varLambda$- and $\varXi$ 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 two- and 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 $\varLambda \varLambda$ hypernuclei, $^{\text { }\text { } \text {}10}_{\varLambda \varLambda }\text {Be}$ [8], $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}$  [9] and especially after publication of the so-called Nagara event  [10, 11], various approaches have been employed to study doubly-strange 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 $\varLambda N$ and $\varLambda \varLambda$ potentials to investigate $\varLambda \varLambda$ s-shell hypernuclei [12]. Thereby, it was assumed that the effects of tensor forces, three-body forces and $\varLambda -\varSigma$ conversion are effectively included in such central potentials by fitting the binding energies of the s-shell (core) $\varLambda$-hypernuclei. In their later work [13], the channel coupling, e.g., $\varLambda N - \varSigma N$ or $\varLambda \varLambda - \varXi 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 $\varLambda \varLambda$ 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 $\varLambda$ interactions can be taken into account. The interactions between a $\varLambda$ 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 $\varLambda N -\varSigma N$ and $\varLambda \varLambda - \varXi 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}_{\varLambda }\text {He}$ and $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {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- ($\varLambda \varLambda \alpha$), four-cluster ($\varLambda \varLambda \alpha \alpha$) and $\varLambda \varLambda n p$ 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 $\varLambda \varLambda$ 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$ $\varLambda \varLambda$ hypernuclei.

Recently, the Jacobi NCSM has been successfully employed by us in studies of single-$\varLambda$ 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., $\varLambda N -\varSigma N$, and $\varLambda \varLambda - \varLambda \varSigma - \varSigma \varSigma - \varXi N$ are explicitly considered. As a 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 $\varLambda \varLambda$ 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  [27, 28, 38,39,40] and for nuclear matter [31, 41]. It is therefore of great interest to study the predictions of the chiral YY potentials for $\varLambda \varLambda$ hypernuclei.

The paper is structured in the following way. We present in Sect. 2 some details of the technical realization. Some relevant formulas are also provided in the appendix. First results for the $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}$, $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$ and $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ hypernuclei are discussed in Sect. 3. Finally, conclusions and an outlook are given in Sect. 4.

2 Numerical realization

2.1 Jacobi basis for $S=-2$ systems

In this section, we generalize our Jacobi no-core shell model (J-NCSM) formalism [28] to $S=-2$ hypernuclei. Adding a second $\varLambda$ hyperon to single-strangeness systems complicates the numerical realization in many ways. All particle conversions that involve a $\varXi$ hyperon, for instance $\varLambda \varLambda \leftrightarrow \varXi N, \, \varSigma \varSigma \leftrightarrow \varXi N \, \text {or} \,\, \varLambda \varSigma \leftrightarrow \varXi 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 $|\alpha ^{*{(Y_1 Y_2)}}\rangle$, and the other that contains a doubly strange $\varXi$ hyperon denoted as $|\alpha ^{*{(\varXi )}}\rangle$. The former is constructed by coupling the antisymmetrized states of $A-2$ nucleons, $|\alpha _{(A-2)N}\rangle$, to the states describing a system of two hyperons, $| Y_1 Y_2 \rangle$

with $Y_1, Y_2 =\varLambda , \varSigma$ and $Y_1 \le Y_2$. Here the inequality $Y_{1} \le Y_2$ indicates the fact that we distinguish among the three two-hyperon states $|\varLambda \varLambda \rangle$, $|\varLambda \varSigma \rangle$ and $|\varSigma \varSigma \rangle$ but do not consider the $|\varSigma \varLambda \rangle$ state explicitly. The notations in Eq. (1) are the same as introduced in Refs. [7, 28]. For example, the symbol $\alpha _{(A-2)N}$ stands for all quantum numbers characterizing the antisymmetrized states of $A-2$ nucleons: the total number of oscillator quanta ${\mathcal {N}}_{A-2}$, total angular momentum $J_{A-2}$, isospin $T_{A-2}$ and state index $\zeta _{A-2}$ as well. Similarly, $\alpha _{Y_1 Y_2}$ stands for a complete set of quantum numbers describing the subcluster of two hyperons $Y_1$ and $Y_2$. The principal quantum number $n_{\lambda }$ of the harmonic oscillator (HO) together with the orbital angular $\lambda$ describe the relative motion of the $(A-2)$N core with respect to the center-of-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 for standard nuclei and single $\varLambda$ 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 $|\alpha ^{*(\varXi )}\rangle$, one combines the antisymmetrized states of an $(A-1)$N system, $|\alpha _{(A-1)N}\rangle$, with the HO states, $|\varXi \rangle$, describing the relative motion of a $\varXi$ hyperon with respect to the C.M. of the (A-1)N subcluster

Here, also $\alpha _{(A-1)N}$ denotes a set of quantum numbers representing an antisymmetrized state of $A-1$ nucleons. The relative motion of a $\varXi$ hyperon is labelled by the HO principal quantum number $n_{\varXi }$, the orbital angular momentum $l_{\varXi }$ and spin $s_{\varXi }=\frac{1}{2}$ which combine together to form the total angular momentum $I_{\varXi }$, and the isospin $t_\varXi =\frac{1}{2}$. Again, following the definition of our coefficients of fractional parentage (CFPs) [7], the momentum or position vector points towards the spectator particle, i.e., the $\varXi$. Finally, the last lines in Eqs. (1,2) also show the graphical representations of the states.

With the basis states defined in Eqs. (1,2), the $S=-2$ hypernuclear wave function $|\varPsi (\pi J T)\rangle$ can be expanded as

$$\begin{aligned} |\varPsi (\pi J T) \big \rangle= & {} \sum _{\alpha ^{*{(Y_1 Y_2)}}} C_{\alpha ^{*(Y_1 Y_2)}} \, \big |\alpha ^{*(Y_1 Y_2)}({\mathcal {N}} J T)\big \rangle \nonumber \\&+ \sum _{\alpha ^{*(\varXi )}} C_{\alpha ^{*(\varXi )}} \,\big | \alpha ^{*(\varXi )}({\mathcal {N}}J T)\big \rangle . \end{aligned}$$

(3)

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 ${\mathcal {N}} = {\mathcal {N}}_{A-2} + 2n_{\lambda } + l_{\lambda } + {\mathcal {N}}_{Y_1 Y_2 } = {\mathcal {N}}_{A-1} + l_{\varXi } + 2n_{\varXi } \le {\mathcal {N}}_{max}$. Of course, by doing so, the computed binding energies will depend on $\mathcal {N}_{max}$ and on HO frequency $\omega$. There are several extrapolation methods that have been well tested for nuclear NCSM calculations [42,43,44,45]. Here we observe, like for the $S=-1$ systems, that the optimal HO frequencies for the $S=-2$ hypernuclear binding energies are generally not the same as the one for the parent nuclear binding energies. Therefore, in order to extract the converged results, we follow the two-step extrapolation procedure that has been extensively employed for the J-NCSM nuclear and single-$\varLambda$ hypernuclear calculations  [7, 28]. The energies $E({\mathcal {N}}, \omega )$ are computed for all accessible model spaces ${\mathcal {N}}_{max}$ and for a wide range of $\omega$. We then determine $E_{{\mathcal {N}}}$ for a given ${\mathcal {N}}_{max} = {\mathcal {N}}$ by minimizing the energies $E({\mathcal {N}}, \omega )$ with respect to $\omega$. In the second step, an exponential fit is applied to $E_{{\mathcal {N}}}$ in order to extrapolate to ${\mathcal {N}} \rightarrow \infty$. For excitation energies and separation energy differences, it is more appropriate to fit the results to a constant as will be discussed in Sect. 3.1. A similar approach has been followed in [28].

2.2 $S=-2$ many-body Hamiltoninan

For the solution of the A-body Schrödinger equation,

$$\begin{aligned} H \ | \varPsi \rangle = E \ | \varPsi \rangle \end{aligned}$$

(4)

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 \ | \varPsi \rangle$ is then straightforward. Using the wave function in Eq. (3), one can write down the final expression for the energy expectation value as follows

$$\begin{aligned}&\langle \varPsi (\pi J T) \,| \,H \, |\, \varPsi (\pi J T) \rangle \nonumber \\&\quad = \sum _{\begin{array}{c} \alpha ^{*(Y_1 Y_2)} \\ \alpha ^{\prime *(Y_1 Y_2)} \end{array}} C_{\alpha ^{*(Y_1 Y_2)}} C_{\alpha ^{\prime *(Y_1 Y_2)}} \, \langle \alpha ^{*(Y_1 Y_2)} |\, H \,|\, \alpha ^{\prime *(Y_1 Y_2)} \rangle \nonumber \\&\qquad +\sum _{\alpha ^{*(\varXi )}, \,\alpha ^{\prime *(\varXi )}} C_{\alpha ^{*(\varXi )}} C_{\alpha ^{\prime *(\varXi )}} \langle \alpha ^{*(\varXi )} \, | H \,|\, \alpha ^{\prime *(\varXi )} \rangle \nonumber \\&\qquad +2 \sum _{\begin{array}{c} \alpha ^{*(Y_1 Y_2)} \\ \alpha ^{\prime *(\varXi )} \end{array}} C_{\alpha ^{*(Y_1 Y_2)}} C_{\alpha ^{\prime *(\varXi )}} \, \langle \alpha ^{*(Y_1 Y_2)} |\, H \,|\, \alpha ^{\prime *(\varXi )} \rangle \,. \end{aligned}$$

(5)

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 $\varXi$ does not contribute to the matrix element $\langle \alpha ^{*(Y_1 Y_2)} |\, H \,|\, \alpha ^{\prime *(Y_1 Y_2)}\rangle$ (in the first line). Likewise, $\langle \alpha ^{*(\varXi )} \, | H \,|\, \alpha ^{\prime *(\varXi )} \rangle$ 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_{Y_1 Y_2}, H_{\varXi }$ and $H^{S=-2}_{Y_1 Y_2 -\varXi 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_{Y_1 Y_2}$ corresponds to a system consisting of $A-2$ nucleons and two singly-strange hyperons $Y_1$ and $Y_2$, and has the following form,

$$\begin{aligned} H_{Y_1 Y_2}= & {} H^{S=0}_{Y_1 Y_2} \,+ H^{S=-1}_{Y_1 Y_2} \,+ H^{S=-2}_{Y_1 Y_2} \nonumber \\= & {} \sum _{i < j=1}^{A-2} \Big ( \frac{2p^{2}_{ij}}{M(t_{Y_1}, t_{Y_2})} \,+ V^{S=0}_{ij} \Big )\nonumber \\&+ \sum _{i=1}^{A-2} \Big ( \frac{m_N + m(t_{Y_1})}{M(t_{Y_1}, t_{Y_2})} \,\frac{p^2_{iY_1}}{2\mu _{iY_1}} \,+ V^{S=-1}_{iY_1}\nonumber \\&+ \frac{m_N + m(t_{Y_2})}{M(t_{Y_1}, t_{Y_2})} \,\frac{p^2_{iY_2}}{2\mu _{iY_2}} \,+ V^{S=-1}_{iY_2} \Big )\nonumber \\&+ \frac{m(t_{Y_1}) + m(t_{Y_2})}{M(t_{Y_1}, {t_{Y2}})}\, \frac{p^{2}_{Y_1 Y_2}}{2\mu _{Y_1 Y_2}} \, +V^{S=-2}_{Y_1 Y_2}\nonumber \\&+ \big (m(t_{Y_1}) + m({t_{Y_2}}) - 2m_{\varLambda }\big ) +\cdots , \end{aligned}$$

(6)

with $Y_1, Y_2 = \varLambda , \varSigma$ and   $Y_1 \le Y_2$. Here, $m({t_{Y_1}}), m(t_{Y_2})$ and $m_N$ are the $Y_1$, $Y_2$ hyperon and nucleon rest masses, respectively. $M(t_{Y_1},t_{Y_2})$ denotes the total rest mass of the system $M(t_{Y_1},t_{Y_2}) = m(t_{Y_1}) + m(t_{Y_2}) + (A-2)m_N$, while $\mu _{i Y_{1}}$ and $\mu _{Y_{1} Y_{2}}$ are the YN and YY reduced masses, respectively. The rest mass differences within the nucleon- and hyperon-isospin multiplets are neglected. $V_{ij}^{S=0}$, $V_{iY}^{S=-1}$, and $V_{YY}^{S=-2}$ 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_{\varXi }$ (involving a $\varXi$ hyperon) corresponds to a system composed of a $\varXi$ hyperon and $A-1$ nucleons. Hence,

$$\begin{aligned} H_{\varXi }= & {} H^{S=0}_{\varXi } \, + \, H^{S=-2}_{\varXi } \nonumber \\= & {} \sum _{i < j=1}^{A-1} \Big ( \frac{2p^{2}_{ij}}{M({\varXi })} \,+ V^{S=0}_{ij} \Big )\nonumber \\&+ \sum _{i=1}^{A-1} \Big ( \frac{m_N + m_{\varXi }}{M({\varXi })} \,\frac{p^2_{\varXi i}}{2\mu _{\varXi i}} \,+ V^{S=-2}_{\varXi i }\Big )\nonumber \\&+ \big (m_{\varXi } + m_N - 2m_{\varLambda }\big )+ \cdots , \end{aligned}$$

(7)

where $m_{\varXi }$ is the $\varXi$ hyperon rest mass and $\mu _{i \varXi }$ is the reduced mass of a $\varXi$ and a nucleon. The total mass of the system is now given by $M(\varXi ) = m_{\varXi } + (A-1)m_N$. $V^{S=-2}_{\varXi i }$ is the $\varXi N$ potential. The ellipses in Eqs. (6,7) stand for those higher-body forces that are omitted here. The transition Hamiltonian $H^{S=-2}_{Y_1 Y_2 ,\varXi N}$ is simply given by the YY-$\varXi$N transition potential

$$\begin{aligned} \begin{aligned} H^{S=-2}_{Y_1 Y_2 ,\varXi N} = \sum _{i=1}^{A-1} V^{S=-2}_{Y_{1} Y_{2} , \varXi i}\,. \end{aligned} \end{aligned}$$

(8)

2.3 Evaluation of the $S=-2$ Hamiltonian matrix elements

Now, taking into account the explicit forms of the A-body Hamiltonian in Eqs. (6-8), all possible contributions to the matrix element $\langle \varPsi (\pi J T) \,| \,H \, |\, \varPsi (\pi J T) \rangle$ can then be split into three groups involving the non-strange $H^{S=0}$, single-strange $H^{S=-1}$ and double-strange $H^{S=-2}$ parts of the total Hamiltonian,

$$\begin{aligned} \langle \varPsi (\pi J T) \,|\, H \,| \, \varPsi (\pi J T) \rangle= & {} \langle \varPsi (\pi J T) \,|\, H^{S=0}\,| \, \varPsi (\pi J T) \rangle \nonumber \\&+ \langle \varPsi (\pi J T) \, |\, H^{S=-1} \, | \, \varPsi (\pi J T) \rangle \nonumber \\&+ \langle \varPsi (\pi J T) \,| \,H^{S=-2} \, |\, \varPsi (\pi J T) \rangle . \end{aligned}$$

(9)

The evaluation of the non-strange part,

$$\begin{aligned}&\langle \varPsi (\pi J T) \,|\, H^{S=0}\,| \, \varPsi (\pi J T) \rangle \nonumber \\&\quad = \sum _{\begin{array}{c} \alpha ^{*(Y_1 Y_2)} \\ \alpha ^{\prime *(Y_1 Y_2)} \end{array}} C_{\alpha ^{*(Y_1 Y_2)}} C_{\alpha ^{\prime *(Y_1 Y_2)}} \, \langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=0}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle \nonumber \\&\qquad +\sum _{\alpha ^{*(\varXi )}, \alpha ^{\prime *(\varXi )}} C_{\alpha ^{*(\varXi )}} C_{\alpha ^{\prime *(\varXi )}} \langle \alpha ^{*(\varXi )} \, | H^{S=0}_{\varXi } | \alpha ^{\prime *(\varXi )} \rangle \,,\nonumber \\ \end{aligned}$$

(10)

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 $\langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=0}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle$ and $\langle \alpha ^{*(\varXi )} \, | H^{S=0}_{\varXi } |\, \alpha ^{\prime *(\varXi )} \rangle$ to the two-nucleon matrix elements in the two-body sector are given by the binomial coefficients of $\left( {\begin{array}{c}A_{\text {nucl}}\\ 2\end{array}}\right) = A_{\text {nucl}}(A_{\text {nucl}}-1)/2$ with $A_{\text {nucl}} =A-2$ and $A_{\text {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,

$$\begin{aligned}&\langle \varPsi (\pi J T) \,| \,H^{S=-2} \, |\, \varPsi (\pi J T) \rangle \nonumber \\&\quad =\sum _{\begin{array}{c} \alpha ^{*(Y_1 Y_2)} \\ \alpha ^{\prime *(Y_1 Y_2)} \end{array}} \!\! C_{\alpha ^{*(Y_1 Y_2)}} C_{\alpha ^{\prime *(Y_1 Y_2)}} \, \langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-2}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle \nonumber \\&\qquad \!+ \sum _{\begin{array}{c} \alpha ^{*(Y_1 Y_2)} \\ \alpha ^{\prime *(\varXi )} \end{array}} \!2 C_{\alpha ^{*(Y_1 Y_2)}} C_{\alpha ^{\prime *(\varXi )}} \, \langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-2}_{Y_1Y_2 ,\varXi N} |\, \alpha ^{\prime *(\varXi )} \rangle \nonumber \\&\qquad \!+\sum _{\alpha ^{*(\varXi )}, \alpha ^{\prime *(\varXi )}} \! C_{\alpha ^{*(\varXi )}} C_{\alpha ^{\prime *(\varXi )}} \langle \alpha ^{*(\varXi )} \, | H^{S=-2}_{\varXi } |\, \alpha ^{\prime *(\varXi )} \rangle ,\nonumber \\ \end{aligned}$$

(11)

are evaluated analogously. Indeed, in order to calculate the last two terms in Eq. (11), one simply needs to expand the states $|\alpha ^{*(\varXi )}\rangle$ in the complete set of intermediate states $|\alpha ^{*(\varXi N)}\rangle$ that explicitly single out a $\varXi N$ pair,

$$\begin{aligned} \begin{aligned} |\alpha ^{*(\varXi )} \rangle = \sum _{\alpha ^{*(\varXi N)} } \langle \alpha ^{*(\varXi )} | \alpha ^{*(\varXi N)} \rangle \,|\alpha ^{*(\varXi N)} \rangle \,. \end{aligned} \end{aligned}$$

(12)

Here the transition coefficients $\langle \alpha ^{*(\varXi )} | \alpha ^{*(\varXi N)} \rangle$ can be computed using the expression Eq. (A.6) in Ref. [28]. It is easy to see that the last term in Eq. (11), $\langle \alpha ^{*(\varXi )} \, | H^{S=-2}_{\varXi } |\, \alpha ^{\prime *(\varXi )} \rangle ,$ differs from the matrix element of the two-body $\varXi N$ Hamiltonian in the $|\varXi N \rangle$ basis by a combinatorial factor of $A-1$. The factor that relates $\langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-2}_{Y_1Y_2 -\varXi N} |\, \alpha ^{\prime *(\varXi )}\rangle$ to the two-body transition potential $V_{Y_1 Y_2 - \varXi N}$ is, however, not obvious because of possible couplings between identical and non-identical two-body states, for instance, $\varSigma \varSigma - \varXi N$ or $\varLambda \varLambda -\varXi N$. In Appendix A, we have shown that, in this case, the corresponding combinatorial factor is $\sqrt{A-1}$ (see Table 4).

2.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),

$$\begin{aligned}&\langle \varPsi (\pi J T) \,|\, H^{S=-1}\,| \, \varPsi (\pi J T) \rangle \nonumber \\&\quad = \sum _{\begin{array}{c} \alpha ^{*(Y_1 Y_2)} \\ \alpha ^{\prime *(Y_1 Y_2)} \end{array}} \! \!\!C_{\alpha ^{*(Y_1 Y_2)}} C_{\alpha ^{\prime *(Y_1 Y_2)}} \, \langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-1}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle ,\nonumber \\ \end{aligned}$$

(13)

in some details since it requires new sets of transition coefficients. Here, in order to compute the matrix elements $\langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-1}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle$, 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 $H^{S=-1}_{Y_1 Y_2}$ in Eq. (6)). It is then instructive to exploit two separate intermediate sets, namely $|\big ( \alpha ^{*{(Y_1 N)}} \big )^{*(Y_2)}\rangle$ and $|\big ( \alpha ^{*{(Y_2 N)}} \big )^{*(Y_1)}\rangle$. The first set, $|\big ( \alpha ^{*{(Y_1 N)}} \big )^{*(Y_2)}\rangle$, is needed when computing the matrix elements of the first two terms of $H^{S=-1}_{Y_1 Y_2}$ where $Y_1$ is the active hyperon while $Y_2$ plays the role of a spectator. Similarly, the second set, $|\big ( \alpha ^{*{(Y_2 N)}} \big )^{*(Y_1)}\rangle$, is useful for evaluating the two remaining terms in Eq. (6) where the roles of the $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}\rangle$, 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 $|\alpha ^{*(Y_1 N)}\rangle$ 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 $|\alpha ^{*(Y_1 Y_2)}\rangle$ in Eq. (1). This in turn allows for the following expansions

$$\begin{aligned}&|\alpha ^{*(Y_1 Y_2)}\rangle \nonumber \\&\quad = \!\!\!\sum _{\begin{array}{c} { } \\ (\alpha ^{*{(Y_1 N)}} )^{*(Y_2)} \end{array}} \!\!\!\!\! \!\!\big \langle \big ( \alpha ^{*{(Y_1 N)}} \big )^{*(Y_2)} \big | \alpha ^{*(Y_1 Y_2)}\big \rangle \, \big |\big ( \alpha ^{*{(Y_1 N)}} \big )^{*(Y_2)}\big \rangle ,\nonumber \\ \end{aligned}$$

(16)

or,

$$\begin{aligned}&|\alpha ^{*(Y_1 Y_2)}\rangle \nonumber \\&\quad = \!\!\!\sum _{\begin{array}{c} { } \\ (\alpha ^{*{(Y_2 N)}} )^{*(Y_1)} \end{array}} \!\!\! \!\!\!\! \big \langle \big ( \alpha ^{*{(Y_2 N)}} \big )^{*(Y_1)} \big | \alpha ^{*(Y_1 Y_2)}\big \rangle \, \big |\big ( \alpha ^{*{(Y_2 N)}} \big )^{*(Y_1)}\big \rangle .\nonumber \\ \end{aligned}$$

(17)

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, $\langle ( \alpha ^{*{(Y_1 N)}})^{*(Y_2)} | \alpha ^{*(Y_1 Y_2)} \rangle$, we make use of another set of auxiliary states, $|(\alpha ^{*(Y_1)} )^{*(Y_2)}\rangle$, that explicitly single out the $Y_1$ and $Y_2$ hyperons. These states are obtained by coupling the hyperon states $|Y_{2}\rangle$ to the basis states of the $((A-2)N + Y_1)$ system, $|\alpha ^{*(Y_1)}\rangle _{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 ${\mathcal {N}}, J$ and T. Exploiting the completeness of the auxiliary states $|\big ( \alpha ^{*{(Y_1)}} \big )^{*(Y_2)}\big \rangle$, the transition coefficient in Eq. (16) then becomes

where a summation over the states $|\big ( \alpha ^{*{(Y_1)}} \big )^{*(Y_2)}\big \rangle$ is implied. The first overlap $\big \langle \big ( \alpha ^{*{(Y_1 N)}} \big )^{*(Y_2)} | \big ( \alpha ^{*{(Y_1)}} \big )^{*(Y_2)}\big \rangle$ in Eq. (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 $\big \langle \big ( \alpha ^{*{(Y_1)}} \big )^{*(Y_2)} | \alpha ^{*(Y_1 Y_2)}\big \rangle$ can quickly be deduced from Eq. (11) in [7],

$$\begin{aligned}&\! \big \langle \big (\alpha ^{*(Y_1)}\big )^{*(Y_2)}| \alpha ^{*(Y_1 Y_2)} \big \rangle \nonumber \\&\quad = \! \delta _{{T}^{\prime }_{A-2} T_{A-2}} \delta _{{J}^{\prime }_{A-2} J_{A-2}} \delta _{\mathcal {{N}^{\prime }}_{A-2} {\mathcal {N}}_{A-2}}\delta _{{\zeta }^{\prime }_{A-2} \zeta _{A-2}}\nonumber \\&\qquad \times {\hat{I}}_{Y_{1}} {\hat{I}}_{Y_{2}} \, {\hat{J}}_{Y_1 Y_2} {\hat{S}}_{Y_1 Y_2}\, {\hat{T}}_{Y_1 Y_2} \,{\hat{I}}_{\lambda } {\hat{J}}^{*(Y_1)}_{A-1}\, {\hat{T}}^{*(Y_1)}_{A-1}\nonumber \\&\qquad \times (-1)^{ 3J_{A-2} + 2T_{A-2} + T_{Y_1 Y_2} + S_{Y_1Y_2} +\lambda + t_{Y_1} + l_{Y_1} + t_{Y_2} + l_{Y_2} + I_{Y_1}}\nonumber \\&\qquad \times \sum _{S_{A-1}=J_{A-2} + s_{Y_1}} \!\! \!\!\! (-1)^{S_{A-1} +1} \, {\hat{S}}^{2}_{A-1} \left\{ \begin{array}{ccc} J_{A-2} &{} s_{Y_1} &{} S_{A-1}\\ l_{Y_1} &{} J^{*(Y_1)}_{A-1} &{} I_{Y_1} \end{array} \right\} \nonumber \\&\qquad \times \,\sum _{\begin{array}{c} {L = l_{Y_1} + l_{Y_2}}\\ {S= S_{A-1} +s_{Y_2}} \end{array}} \!\!\!\!\! \!\!\!\!{\hat{L}}^2 {\hat{S}}^2 \left\{ \begin{array}{ccc} l_{Y_1} &{} S_{A-1} &{} J^{*(Y_{1})}_{A-1}\\ l_{Y_2} &{} s_{Y_2} &{} I_{Y_2}\\ L &{} S &{}{J} \end{array}\right\} \! \left\{ \begin{array}{ccc} l_{Y_1 Y_2} &{} S_{Y_1 Y_2} &{} J_{Y_1 Y_2}\\ \lambda &{} J_{A-2} &{} I_{\lambda }\\ L &{} S &{} {J} \end{array}\!\right\} \nonumber \\&\qquad \! \times \left\{ \begin{array}{ccc} s_{Y_2} &{} s_{Y_1} &{} S_{Y_1 Y_2}\\ J_{A-2} &{} S &{} S_{A-1} \end{array} \right\} \left\{ \begin{array}{ccc} t_{Y_2} &{} t_{Y_1} &{} T_{Y_1Y_2}\\ T_{A-2} &{} {T} &{} T^{*(Y_1)}_{A-1} \end{array} \right\} \nonumber \\&\qquad \times \langle n_{Y_1} \,l_{Y_1}\, n_{Y_2}\, l_{Y_2} : L \,|\, n_{Y_1 Y_2} \,l_{Y_1Y_2} \, n_{\lambda } \, \lambda : L \rangle _{d}\,, \end{aligned}$$

(20)

with,

$$ d=\frac{(A-2) m_N \, m(t_{Y_2})}{ m(t_{Y_1}) \big ((A-2)\,m_N + m(t_{Y_1}) + m(t_{Y_2})\big )}\,.$$

Here, we use the notation ${{\hat{j}}} = \sqrt{2j+1}$ and abbreviate the summations running from $|J_1-J_2|$ to $J_1+J_2$ simply by $J_1+J_2$.

The transition coefficients for the second expansion in Eq. (17) are computed analogously. Taking into account the expansions Eqs. (16,17), the matrix element $\langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-1}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle$ in Eq. (13) is then decomposed into,

$$\begin{aligned}&\langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-1}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle \nonumber \\&\quad = \langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-1}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle _{Y_2} \nonumber \\&\qquad +\,\, \langle \alpha ^{*(Y_1 Y_2)} |\, H^{S=-1}_{Y_1Y_2} |\, \alpha ^{\prime *(Y_1 Y_2)} \rangle _{Y_1 }\,. \end{aligned}$$

(21)

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 elements and 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.

3 Results

In this section, as a first application, we report results for the $\varLambda \varLambda$ s-shell hypernuclei $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}(1^+, 0)$, $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}(\frac{1}{2}^+, \frac{1}{2})$, and $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}(0^+, 0)$. To zeroth approximation, these systems can be regarded as a $\varLambda \varLambda$ pair in the $^1S_0$ state being attached to the corresponding core-nuclei predominantly in their ground states. While the quantum numbers of $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$, $(J^+, T) = (\frac{1}{2}^+, \frac{1}{2}),$ are obvious, those for the $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {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 $\varLambda _{N}= 450$ MeV ($\hbox {N}^{4}$LO+(450))  [33], SRG-evolved to $\lambda _{NN} =1.6$ ${\mathrm{fm}}^{-1}$, is adopted for describing the NN interaction. The next-to-leading order potential NLO19  [4] with a chiral cutoff of $\varLambda _{Y} = 650$ MeV and an SRG parameter of $\lambda _{YN} =0.868$ ${\mathrm{fm}}^{-1}$ is used for the YN interaction. We remark that the chosen NN and YN potentials successfully predict the empirical $\varLambda$-separation energies for $^3_{\varLambda }\text {H}$, $^4_{\varLambda }\text {He}(1^+)$ and $^5_{\varLambda }\text {He}$, and underbind $^4_{\varLambda }\text {He}(0^+)$ only slightly [28]. Therefore, those potentials are an excellent starting point for the extension to $S=-2$. Eventually, in a future study, it will be interesting to also examine the dependence of $S=-2$ hypernuclei on the SRG evolution and the starting interactions in the $S=-1$ sector. Clearly, in such an investigation one has to ensure to maintain the favorable description of $S=-1$ hypernuclei. Therefore, we do not expect a significant impact of any variations subject to that pre-condition on separation energies for $S=-2$ hypernuclei. For 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 $\varLambda _{YY} =600$ MeV. Therefore, the predictions for $S=-2$ hypernuclei shown here are based on a set of interactions that are consistent with the available NN, YN and YY data and with the empirical separation energies of light $S=-1$ hypernuclei.

One of our primary aims here is to establish the predictions of these chiral YY potentials for double-$\varLambda$ 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 ($\varLambda \varLambda$, $\varXi N$, ...). Due to the latter circumstance, in the chiral approach (as well as in meson-exchange and/or constituent quark models) the assumption of $\hbox {SU(3)}_{\mathrm{f}}$ symmetry is an essential prerequisite for deriving pertinent potentials. For example, in chiral EFT the short-distance 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 $^1S_0$ partial wave of $\varLambda \varLambda$. They correspond to the SU(3) singlet irreducible representation, see Ref. [30], and are denoted by ${{\tilde{C}}}^{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 $\varXi ^{-} p - \varLambda \varLambda$  [46] and the upper limits of elastic and inelastic $\varXi ^{-} p$ cross sections [47]). 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 $\varLambda \varLambda$ $^1S_0$ scattering lengths predicted by these interactions are compatible with values inferred from empirical information  [48, 49]. 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 $\hbox {SU(3)}_{\mathrm{f}}$ component which means essentially only in the strength of the $\varXi N$ interaction in the $^3S_1$ partial wave. This has an impact on the corresponding in-medium properties of the $\varXi$. Specifically, the updated version from 2019  [31] yields a moderately attractive $\varXi$ single-particle potential that is roughly in line  [50] with recent experimental evidence that the existence of bound $\varXi$-hypernuclei is very likely  [51]. With regard to $\varLambda \varLambda$ systems, we observe that the two realizations yield very similar binding energies for the double-$\varLambda$ s-shell hypernuclei. This indicates that, in general, the actual strength of the spin-triplet $\varXi N$ interaction has little influence on few-body observables related to $\varLambda \varLambda$. In the following, we therefore present results for the LO and the updated NLO interactions for a chiral cutoff of $\varLambda _{YY} =600$ MeV. In order to speed up the convergence, both YY potentials are also SRG-evolved. We use a wide range of the SRG flow parameters, namely $1.4 \le \lambda _{YY} \le 3.0$ ${\mathrm{fm}}^{-1}$, to quantify the contribution of possible SRG-induced YYN three-body forces.

3.1 $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}(0^+, 0)$

The $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}$ hypernucleus is so far the lightest double-$\varLambda$ system being unambiguously established. Since the observation of the Nagara event [10], its $\varLambda \varLambda$ separation energy, defined as $B_{\varLambda \varLambda }( ^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}) = E(^4\text {He}) - E( ^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He})$, has been exploited as a crucial constraint for constructing effective potentials that are then employed in many-body calculations like the Gaussian expansion method  [13, 52] or the cluster Faddeev-Yakubovsky approach  [16, 17]. The re-analysis of the Nagara event using the updated $\varXi$ mass yielded a slightly smaller $\varLambda \varLambda$ separation energy, $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}) =6.91 {\pm } 0.16$ MeV  [11, 53], as compared to the initially estimated value of $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}) = 7.25 {\pm } 0.19$ [10]. This, in turn, may have direct consequences for theoretical predictions for potentially observable bound states of other s-shell $\varLambda \varLambda$ hypernuclei, particularly the $A=4$ double-$\varLambda$ system  [25, 54], see also the discussion in 3.3. We note that the information about $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {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-$\varLambda$ system using the two chiral interactions to scrutinize their consistency with the measured $\varLambda \varLambda$ separation energy.

Fig. 1

Binding energy E, $\varLambda \varLambda$-separation energy $B_{\varLambda \varLambda }$ and $\varLambda \varLambda$-excess energy $\varDelta B_{\varLambda \varLambda }$ for $^{\text { }\text { }6}_{\varLambda \varLambda }$He computed using the YY NLO(600) interaction [31] that is SRG evolved to a flow parameter of $\lambda _{YY} =1.8$ ${\mathrm{fm}}^{-1}$. The SMS $\hbox {N}^{4}$LO+(450) potential [33] with $\lambda _{NN}=1.6$ ${\mathrm{fm}}^{-1}$ and the NLO19(650) potential  [4] with $\lambda _{YN} =0.868$ ${\mathrm{fm}}^{-1}$ are employed for the NN and YN interactions, respectively. (a): Solid lines with different colors and symbols represent numerical results for different model spaces ${\mathcal {N}}$. Dashed lines are obtained using the ansatz Eq. (22) in [28]. (b-d): Red (horizontal) lines with shaded areas indicate the converged results and the corresponding uncertainties. Note that binding energies strictly converge from above whereas binding energy differences like separation and excess energies do not necessarily fulfill this constraint

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 $\omega$- and ${\mathcal {N}}$-space extrapolations for $E(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He})$ are illustrated in panels (a) and (b) of Fig. 1, respectively. The error bars shown in the figures of the ${\mathcal {N}}$-dependence of energies, panel (b), are given by the difference to the next model space. These error bars are not meant to provide a realistic uncertainty estimate but only to give weights for the following extrapolation to ${\mathcal {N}}\rightarrow \infty$. For illustration purposes, we present results for the NLO potential with $\lambda _{YY}=1.8$ ${\mathrm{fm}}^{-1}$ but stress that the convergence trend is similar for all other values of $\lambda _{YY}$, and for the LO interaction. Also, the behavior of $E(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He})$ with respect to $\omega$ and ${\mathcal {N}}$ resembles that of the binding energy of the parent hypernucleus $^5_{\varLambda }\text {He}$ [28]. Furthermore, panel (b) clearly demonstrates a nice convergence pattern of the binding energy $E(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He})$ computed for model spaces up to ${\mathcal {N}}_{max}=14$. We also perform an exponential fit to extrapolate the $\varLambda \varLambda$-separation energy, as done for $S=-1$ systems. Clearly, the result for $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He})$, displayed in panel (c), is also well-converged for ${\mathcal {N}}_{max}=14$ (practically with the same speed as that of $E(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He})$). Note that, for single-$\varLambda$ hypernuclei, the separation energy $B_{\varLambda }$ converges somewhat faster than the individual binding energies. For $S=-2$ systems, we are also interested in the so-called $\varLambda \varLambda$ excess binding energy

$$\begin{aligned} \begin{aligned} \varDelta B_{\varLambda \varLambda } (^{\text { }\text { } A}_{\varLambda \varLambda } \text {X} )&= B_{\varLambda \varLambda }(^{\text { }\text { } A}_{\varLambda \varLambda } \text {X} ) - 2 {{\bar{B}}}_{\varLambda }(^{A-1}_{\varLambda }\text {X})\\&\quad = 2 {{\bar{E}}}(^{A-1}_{\varLambda }\text {X}) - E(^{\text { }\text { }A}_{\varLambda \varLambda } \text {X}) - E(^{A-2}\text {X}) \end{aligned} \end{aligned}$$

(24)

which provides information about the strength of the $\varLambda \varLambda$ interaction. ${{\bar{B}}}$ and ${{\bar{E}}}$ are spin averaged $\varLambda$-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 $\varLambda$-core interaction, dynamical changes in the core-nucleus structure as well as the mass-polarization effect  [8, 15]. For $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }$He, the spin-dependent part of the $\varLambda$-core interaction vanishes because of the spin zero the parent nucleus $^4$He, hence the difference

$$\varDelta B_{\varLambda \varLambda } (^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He} ) = B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He} ) - 2B_{\varLambda }(^{5}_{\varLambda }\text {He}),$$

will reflect the net contributions of the $\varLambda \varLambda$ interactions and the $^4\text {He}$ core-distortion¹ (polarization) effects. In panel (d), we exemplify the model-space extrapolation for $\varDelta B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He} )$. We observe that $\varDelta B_{\varLambda \varLambda }$ converges with respect to ${\mathcal {N}}$ visibly faster than both, the $\varLambda \varLambda$-separation and the binding energies. For this quantity it is more appropriate to fit to a constant in order to determine the large ${\mathcal {N}}$ extrapolation.

Fig. 2

$B_{\varLambda \varLambda }(^{\text { }\text { }6}_{\varLambda \varLambda } \text {He} )$ (left) and $\varDelta B_{\varLambda \varLambda }{(^{\text { }\text { } 6}_{\varLambda \varLambda } \text {He} )}$ (right) as functions of the flow parameter $\lambda _{YY}$. Calculations are based on the YY LO(600) (blue triangles) and NLO(600) (red circles) potentials. Dash-dotted line with grey band represents the experimental value and the uncertainty of the Nagara event [11]. Same NN and YN interactions as in Fig. 1

Being able to accurately extract $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He} )$ and $\varDelta B_{\varLambda \varLambda }{(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He} )}$, we are in a position to study the impact of the two chiral interactions on these quantities. The converged results for $B_{\varLambda \varLambda }$ and $\varDelta B_{\varLambda \varLambda }$, calculated for a wide range of the SRG flow parameter $\lambda _{YY}$, 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 $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He}$ as can be seen in the left panel. On the other hand, the moderately attractive NLO interaction predicts a $\varLambda \varLambda$ excess energy of $\varDelta B_{\varLambda \varLambda } \approx 1.1$ MeV, that is only slightly larger than the empirical value of $\varDelta B^{exp}_{\varLambda \varLambda } = 0.67 {\pm } 0.17$ MeV  [11, 53]. For completeness, let us mention that the pertinent $\varLambda \varLambda$ $^1S_0$ scattering lengths are $a=-1.52$ fm (LO [29]) and $a=-0.66$ fm (NLO [30]), respectively.

It is rather remarkable that both, $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He} )$ and $\varDelta B_{\varLambda \varLambda }{(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He} )}$, exhibit a rather weak dependence on the SRG YY parameter $\lambda _{YY}$. With an order of 100 keV, it is at least one order of magnitude smaller than the variation of, say, $B_{\varLambda }(^5_{\varLambda }\text {He})$ with respect to the SRG YN flow parameter $\lambda _{YN}$  [28]. The insensitivity of the $\varLambda \varLambda$-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 $\varLambda \varLambda$ interaction.

Finally, we benchmark the probabilities of finding one $\varSigma$ $(P_{\varLambda \varSigma })$ or two $\varSigma$ $(P_{\varSigma \varSigma })$, or the $\varXi$ hyperon $(P_{\varXi })$ in the ground-state wave function of $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He}$ obtained for the two chiral potentials. Such probabilities are summarized in Table 1 for several values of $\lambda _{YY}$. Overall, the $P_{\varLambda \varSigma }$ and $P_{\varSigma \varSigma }$ 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 $\varSigma$ in $^5_{\varLambda }\text {He}$ for the employed NN and YN interactions is also very small, $P_{\varSigma }(^5_{\varLambda }\text {He}) = 0.07 \%$. In contrast, $P_{\varXi }$ is more sensitive to the evolution and also strongly influenced by the interactions. Surprisingly, the updated NLO potential, that yields a more attractive $\varXi$-nuclear interaction [31], predicts a considerably smaller $\varXi$ probability (less than 0.2 $\%$ for $\lambda _{YY} =3.0$ ${\mathrm{fm}}^{-1}$) as compared to the value of $P_{\varXi } =1.1 \%$ obtained for the LO at the same $\lambda _{YY}$. 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 $(\varSigma , \varXi )$ and the interaction strength.

Table 1 Probabilities $(\%)$ of finding a single and double $\varSigma$, and a $\varXi$ hyperons in the ground-state wavefunction of $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}$. Note that $P_{\varSigma }(^{5}_{\varLambda }\text {He})=0.07 \%$

Fig. 3

Binding energy E, $\varLambda \varLambda$-separation energy $B_{\varLambda \varLambda }$ and $\varLambda \varLambda$-excess $\varDelta B_{\varLambda \varLambda }$ for $^{\text { }\text { }5}_{\varLambda \varLambda }$He computed using the YY NLO(600) interaction that is SRG evolved to a flow parameter of $\lambda _{YY} =1.8$ ${\mathrm{fm}}^{-1}$. Same notation, NN and YN interactions as in Fig. 1

3.2 $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}(\frac{1}{2}^+, \frac{1}{2})$

The next system that we investigate is $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$. Although the existence of $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$ has not been experimentally confirmed yet, most of the many-body calculations employing effective potentials that reproduce the separation energy $B_{\varLambda \varLambda }(^{\text { }\text { }6}_{\varLambda \varLambda } \text {He})$ predict a particle-stable bound state of $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$  [13, 16, 54]. However, there are visible discrepancies among the values of $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {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_{\varLambda \varLambda }$ for the $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$ ($^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {H}$) and $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {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_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He})$ instead.

Fig. 4

$B_{\varLambda \varLambda }(^{\text { }\text { } 5}_{\varLambda \varLambda } \text {He} )$ (left) and $\varDelta B_{\varLambda \varLambda }{(^{\text { }\text { } 5}_{\varLambda \varLambda } \text {He} )}$ (right) as functions of the flow parameter $\lambda _{YY}$. Calculations are based on the YY LO(600) (blue triangles) and NLO(600) (red circles) potentials. Same NN and YN interactions as in Fig. 1

The $\omega$- and ${\mathcal {N}}$-extrapolation of the binding energy E, $\varLambda \varLambda$-separation energy $B_{\varLambda \varLambda }$ and the $\varLambda \varLambda$-excess energy $\varDelta B_{\varLambda \varLambda }$ of $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$ are illustrated in Fig. 3. Here, the results are shown for the NLO potential with a flow parameter of $\lambda _{YY}=1.8$ ${\mathrm{fm}}^{-1}$ and for model spaces up to ${\mathcal {N}}_{max} =16$. Note that in the case of $^4_{\varLambda }\text {He}$, the energy calculations were performed for model spaces up to ${\mathcal {N}}_{max}(^4_{\varLambda }\text {He})=22$ in order to achieve a good convergence. Calculations with such large model spaces are currently not feasible for $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$ because of computer-memory constraints. Nonetheless, the illustrative results in Fig. 3 clearly indicate that well-converged results are achieved for this double-$\varLambda$ hypernucleus already for model spaces up to ${\mathcal {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

$$\begin{aligned} \varDelta B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}) = B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}) -2 {B}_{\varLambda }(^4_{\varLambda }\text {He}) \ , \end{aligned}$$

(25)

we do not simply assign the ground-state $\varLambda$-separation energy $B_{\varLambda }(^4_{\varLambda }\text {He}, 0^+)$ to $B_{\varLambda }(^4_{\varLambda }\text {He})$ but rather use a spin-averaged value ${\overline{B}}_{\varLambda }(^4_{\varLambda }\text {He})$ of the ground-state doublet [15]

$$\begin{aligned} {\overline{B}}_{\varLambda }(^4_{\varLambda }\text {He}) = \frac{1}{4} B_{\varLambda }(^4_{\varLambda }\text {He}, 0^+) + \frac{3}{4} B_{\varLambda }(^4_{\varLambda }\text {He}, 1^+), \end{aligned}$$

(26)

with $B_{\varLambda }(^4_{\varLambda }\text {He}, 0^+ (1^+)) =1.708\, (0.904)$ MeV for the employed NN and YN potentials [28]. By doing so, the computed quantity $\varDelta B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He})$ will be less dependent on the spin-dependence effect of the $\varLambda$-core interactions, and, therefore, can be used as a measure of the $\varLambda \varLambda$ interaction strength, provided that the nuclear contraction effects are small. The results for $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He})$ and $\varDelta B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He})$ calculated for the two interactions and a wide range of flow parameter, $1.4 \le \lambda _{YY} \le 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 $^{\text { }\text { } 6}_{\varLambda \varLambda }\text {He}$, reinforcing the insignificance of SRG-induced YYN forces. Again, the LO interaction predicts a much larger $\varLambda \varLambda$-separation energy and a more significant $\varLambda \varLambda$ interaction strength than the one at NLO. In either case, the $\varLambda \varLambda$ excess energy $\varDelta B_{\varLambda \varLambda }$ computed for $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$, slightly exceeds the corresponding one for $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {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 $\varLambda \varLambda -\varXi N$ coupling in $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}$ as discussed in  [18, 55, 56]. 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 $\varDelta B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He} ) < \varDelta B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He})$. As a consequence, our results do also not fit into the correlation of $\varDelta B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He} )$ and $\varDelta B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He})$ shown in the same work. We will need to study more interactions in the future to understand whether such a correlation can also be established using chiral interactions.

Table 2 Probabilities (in percentage) of finding a $\varSigma$ $(P_{\varLambda \varSigma })$, double $\varSigma$ $(P_{\varSigma \varSigma })$and a $\varXi$ $(P_{\varXi })$ hyperons in $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$. $P_{\varSigma }(^4_{\varLambda }\text {He}) = 0.43\,\%$

It is also very interesting to point out that the $\varLambda \varLambda$-separation energies $B_{\varLambda \varLambda }$ for both $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}$ and $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}$ predicted by the NLO potential are surprisingly close to the results obtained by Nemura et al., $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda }\text {He}) = 3.66$ MeV, $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}) = 7.54$ MeV, using the modified Nijmegen YY potential ($\hbox {mND}_{\mathrm{s}}$)  [13]. Finally, we provide in Table 2 the probabilities of finding a $\varSigma$ ($P_{\varLambda \varSigma }$), double $\varSigma$ $(P_{\varSigma \varSigma })$, or a $\varXi$ $(P_{\varXi })$ in the $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda } \text {He}$ ground-state wave function, computed with the two potentials and several SRG values, $\lambda _{YY}=1.4, 2.0$ and 3.0 ${\mathrm{fm}}^{-1}$. Apparently, all the probabilities including also $P_{\varXi }$ exhibit a rather weak sensitivity to the flow parameter $\lambda _{YY}$. The two interactions seem to have little impact on the $\varSigma$-probabilities ($P_{ \varLambda \varSigma }$ and $P_{\varSigma \varSigma }$) but strongly influence $P_{ \varXi }$. Like in the $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda } \text {He}$ system, here, the LO potential yields considerably larger $\varXi$-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 $\varSigma$ or $\varXi$ hyperon in $^{\text { }\text { }\text { } \text {}5}_{\varLambda \varLambda } \text {He}$ are visibly larger than the corresponding ones in $^{\text { }\text { }6}_{\varLambda \varLambda } \text {He}$. This is indeed consistent with the $\varSigma$-probabilities in the ground-state wave functions of their parent hypernuclei (e.g., $P_{\varSigma }(^4_{\varLambda }\text {He} ) =0.43\,\%$ and $P_{\varSigma }(^5_{\varLambda }\text {He} ) =0.07\, \%$), and more importantly, is consistent with the suppression of particle conversions such as $\varLambda \varLambda - \varXi N$ in p-shell hypernuclei  [55].

Fig. 5

(a): Ground-state energies of $^{\text { }\text { } 4}_{\varLambda \varLambda }\text {He}$ as functions of $\omega$ for model space ${\mathcal {N}}=10 - 32$. Calculations are performed with the YY NLO(600) potential evolved to a flow parameter of $\lambda _{YY} =1.8$ ${\mathrm{fm}}^{-1}$. (b): model space extrapolation of $E(^{\text { }\text { }4}_{\varLambda \varLambda }\text {H})$ with the same YY interaction as in (a). (c): model space extrapolation of $E(^3_{\varLambda }\text {H})$. (d): Converged $E(^{\text { }\text { } 4}_{\varLambda \varLambda }\text {H})$ as functions of the flow parameter for the LO(600) (blue triangles) and NLO(600) (red circles) potentials. The dashed line with grey band represents the computed $E(^3_{\varLambda }\text {H})$ and the theoretical uncertainty, respectively. Same NN and YN interactions as in Fig. 1

3.3 $^{\text { }\text { }\text { }\text {}4}_{\varLambda \varLambda }\text {H}(1^+, 0)$

Our final exploratory s-shell hypernucleus is $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$. This system has been the subject of many theoretical and experimental studies. It turned out that theoretical predictions of the stability of $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ against the $^3_{\varLambda }\text {H} + \varLambda$ breakup are very sensitive to the interpretations of double-strangeness hypernuclear data, in particular, the $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}$ hypernucleus [54]. Indeed, Nemura et al.  [13] observed a particle-stable but loosely bound state of $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ (just only about 2 keV below the $^3_{\varLambda }\text {H} + \varLambda$ threshold for the $\hbox {mND}_{\mathrm{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_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}) = 7.25 {\pm } 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 $\varLambda \varLambda p n$ system, but a particle-stable $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ hypernucleus when solving the (three-body) Faddeev equation for the $\varLambda \varLambda 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 $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ is not compatible with the corrected value of $B_{\varLambda \varLambda }(^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}) = 6.91 {\pm } 0.16$ MeV. Although the observation of $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ was reported in an experiment at BNL [57], it has been recently invalidated by a thorough re-examination of the recorded events  [58]. Nevertheless, the existence of a stable $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {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$ $\varLambda \varLambda$ hypernuclei as the $\hbox {mND}_{\mathrm{s}}$ interaction [13], also results in a loosely bound state for $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {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_\varLambda$H and the ground-state energy of $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ are crucial. In panels (a) and (b) of Fig. 5, we examine the convergence of $E(^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H})$ in $\omega$- and ${\mathcal {N}}$-space, respectively, using model spaces to ${\mathcal {N}}_{max} =32$. The results are shown for the NLO(600) potential with a flow parameter of $\lambda _{YY}=2.4$ ${\mathrm{fm}}^{-1}$. For a better comparison, the ${\mathcal {N}}$-space extrapolation of $E(^3_{\varLambda }\text {H})$, computed with model spaces up to ${\mathcal {N}}=32$, is also presented in panel (c). As expected, due to the weak binding of the hypertriton, the binding energy calculations for both hypernuclei, $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ and $^3_{\varLambda }\text {H}$, converge very slowly when using HO bases. It also clearly sticks out that the optimal HO frequencies $\omega$ for large model space sizes are around $\omega _{opt} \approx 6$ MeV which is much smaller than the value of $\omega _{opt} \approx 16$ MeV for the $A=4,5$ systems. This again reflects the large spatial extension of the wave functions of $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ and $^3_{\varLambda }\text {H}$. Nevertheless, one can still observe a slightly faster convergence speed for $E(^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H})$ (especially with the LO potential) as for $E(^3_{\varLambda }\text {H})$. Moreover, our extrapolated value of $E(^3_{\varLambda }\text {H}) = -2.314 {\pm } 0.009$ MeV (for model space up to ${\mathcal {N}}=36$) agrees within 10 keV with the exact Faddeev result $E_{\mathrm{Fad}}(^3_{\varLambda }\text {H}) = -2.333 {\pm } 0.002$ MeV [27]. We conclude that a model space truncation of ${\mathcal {N}}_{max} =32$ for the energy calculations in $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ should be sufficient in order to draw conclusions about the stability of the system against $\varLambda$ emission.

The extrapolated ground-state energies $E(^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H})$ for the NLO (red circles) and LO (blue triangles) potentials evolved to a wide range of flow parameters are displayed in panel (d) of Fig. 5. Here, the dashed black line together with the grey band represent the computed $E(^3_{\varLambda }\text {H})$ and the estimated uncertainty. Calculations with the NLO potential seem to converge more slowly than the ones for the LO interaction. The NLO potential clearly leads to an unbound $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {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 $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ state. The LO results for $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ likely hint at a particle-unstable system with respect to the hypertriton $^3_{\varLambda }$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 $^{\text { }\text { }\text { } \text {}6}_{\varLambda \varLambda }\text {He}$, very likely it overpredicts the actual attraction in the $A=4$ system, too. Interestingly, in pionless EFT [25] a $\varLambda \varLambda$ scattering length practically identical to that of our LO interaction was found as limit for which the $^{\text { }\text { }\text { } \text {}4}_{\varLambda \varLambda }\text {H}$ system becomes bound.

4 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 $\varLambda \varLambda$ s-shell hypernuclei. In the actual calculation, the YY forces are combined with a set of BB interactions that is consistent with all available NN, $\varLambda N$ and $\varSigma N$ scattering data and with the empirical separation energies of light $S=-1$ hypernuclei. To speed up the convergence, the YY 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 $\varLambda \varLambda$-separation energies, implying negligible contributions of SRG-induced YYN forces. We found that the binding energy for $^{\text { } \text { }\text { }\text {} 6}_{\varLambda \varLambda }$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 $^{\text { } \text { }\text { }\text {} 5}_{\varLambda \varLambda }$He hypernucleus, whereas $^{\text { } \text { }\text { }\text {} 4}_{\varLambda \varLambda }$H is found to be unstable against a breakup to $^3_{\varLambda }\text {H} + \varLambda$. 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. Work in this direction is in progress. It will be also very interesting to study the predictions of the chiral YY interactions for other s-shell $\varLambda \varLambda$ systems such as $^{\text { } \text { }\text { }\text {} 4}_{\varLambda \varLambda }$n or $^{\text { } \text { }\text { }\text {} 4}_{\varLambda \varLambda }\text {He}$, as well as for p-shell hypernuclei. Finally, investigating possible Tjon-line like correlations for $B_{\varLambda \varLambda }$ 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 $\varSigma \varSigma \rightarrow N \varXi$. 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 two-body 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, Miyagawa and Kamada  [59, 60] 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. [61] have formulated the problem (for the ${{\bar{K}}} NN$-$\pi \varSigma N$ system) taking all permutations of particles explictly into account. This is however not consistent with the approach of BB interactions used in  [29, 30, 62]. 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,

$$\begin{aligned} {H} = \sum _{k_1 k_1^{\prime }} {T}_{ k_1^{\prime } k_1} a^{\dagger }_{k_1^{\prime }} a_{k_{1}} + \frac{1}{2} \sum _{\begin{array}{c} {k_1 k_2} \\ { k_1^{\prime } k_2^{\prime }} \end{array}} V_{k_1^{\prime } k_2^{\prime }, k_1 k_2} a^{\dagger }_{k_1^{\prime }} a^{\dagger }_{k_2^{\prime }} a_{k_2} a_{k_{1}},\nonumber \\ \end{aligned}$$

(A.1)

where $k_i$ stands for a set of quantum numbers characterizing the particle state, i.e., momentum, spin, isospin as well as particle species $\lambda _{i} \,(\text {N}, \varLambda , \varSigma \,\text {or}\, \varXi )$. When it is necessary to separate the particle species $\lambda _{i}$ from other quantum numbers, we use $k_i = \lambda _i \, {\tilde{k}}_i$. Let us further assume that the potential matrix elements $V_{k_1^{\prime } k_2^{\prime }, k_1 k_2}$ in Eq. (A.1) are antisymmetric under exchanges of two indices, i.e., $V_{k_1^{\prime } k_2^{\prime }, k_1 k_2 } = -V_{k_1^{\prime } k_2^{\prime }, k_2 k_1} = -V_{k_2^{\prime } k_1^{\prime }, k_1 k_2}= V_{k_2^{\prime } k_1^{\prime }, k_2 k_1}$. 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^{\prime }_1$ and $k^{\prime }_2$ in Eq. (A.1).

1.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

$$\begin{aligned} |\{p_1 p_2 \} \rangle \equiv a^{\dagger }_{p_{1}} a^{\dagger }_{p_2} | 0 \rangle = \frac{1}{\sqrt{2}} \big (|p_1\rangle | p_2\rangle - |p_2\rangle | p_1 \rangle \big ), \end{aligned}$$

(A.2)

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 $\lambda _1 \ne \lambda _2$ reads

$$\begin{aligned}&\sum _{p_1< p_2} |\{p_1 p_2 \} \rangle \langle \{p_1 p_2 \} | \nonumber \\&\quad \equiv \sum _{\lambda _1 < \lambda _2} \int d^3{{\tilde{p}}_1} d^3 {\tilde{p}}_2 |\{\lambda _1 {\tilde{p}}_1 \lambda _2 {\tilde{p}}_2 \} \rangle \langle \{\lambda _1 {\tilde{p}}_1 \lambda _2{\tilde{p}}_2 \} | = \mathbb {1},\nonumber \\ \end{aligned}$$

(A.3)

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

$$\begin{aligned}&\!\! \sum _{p_1< p_2} |\{p_1 p_2 \} \rangle \langle \{p_1 p_2 \} |\nonumber \\&\quad = \frac{1}{2} \Big \{ \sum _{p_1< p_2} |\{p_1 p_2 \} \rangle \langle \{p_1 p_2 \} | + \sum _{p_1< p_2} |\{p_1 p_2 \} \rangle \langle \{p_1 p_2 \} | \Big \}\nonumber \\&\quad = \frac{1}{2} \Big \{ \sum _{p_1 < p_2} |\{p_1 p_2 \} \rangle \langle \{p_1 p_2 \} | + \sum _{p_1 > p_2} |\{p_2 p_1 \} \rangle \langle \{p_2 p_1 \} | \Big \}\nonumber \\&\quad = \frac{1}{2} \sum _{p_1, p_2} |\{p_1 p_2 \} \rangle \langle \{p_1 p_2 \} |\,. \end{aligned}$$

(A.4)

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 $\frac{1}{2}$. For the case of two identical particles, i.e., $\lambda _{1} = \lambda _{2}$, the completeness relation becomes

$$\begin{aligned} \frac{1}{2} \int d^3{{\tilde{p}}_1} d^3 {\tilde{p}}_2 |\{\lambda _1 {\tilde{p}}_1 \lambda _1 {\tilde{p}}_2 \} \rangle \langle \{\lambda _1 {\tilde{p}}_1 \lambda _1 {\tilde{p}}_2 \} | = \mathbb {1} \end{aligned}$$

(A.5)

following similar lines. The factor $\frac{1}{2}$ can also be absorbed into the definition of the states when one rewrites Eq. (A.5) as follows

$$\begin{aligned} \int d^3{{\tilde{p}}_1} d^3 {\tilde{p}}_2 \frac{1}{\sqrt{2}}|\{\lambda _1 {\tilde{p}}_1 \lambda _1 {\tilde{p}}_2 \} \rangle \langle \{\lambda _1 {\tilde{p}}_1 \lambda _1 {\tilde{p}}_2 \} | \frac{1}{\sqrt{2}} = \mathbb {1}. \end{aligned}$$

(A.6)

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

$$\begin{aligned}&\!\!\langle \{p_1^{\prime } p_2^{\prime }\} | {T} | \{ p_1 p_2 \} \rangle \nonumber \\&\quad = \delta _{p_1^{\prime } p_1 } {T}_{p_2^{\prime } p_2} - \delta _{p^{\prime }_{1} p_2} {T}_{p_2^{\prime }p_1 } + \delta _{p_2^{\prime } p_2 }{T}_{p_1^{\prime } p_1} - \delta _{p^{\prime }_{2} p_1}{T}_{p_1^{\prime }p_2 } \nonumber \\&\quad = \delta _{p^{\prime }_1 p_1} \delta _{p^{\prime }_2 p_2} {t}_{p^{\prime }_2} - \delta _{p^{\prime }_1 p_2} \delta _{p^{\prime }_2 p_1} {t}_{p^{\prime }_2} + \delta _{p^{\prime }_2 p_2} \delta _{p^{\prime }_1 p_1} {t}_{p^{\prime }_1} \nonumber \\&\quad \qquad - \delta _{p^{\prime }_2 p_1} \delta _{p^{\prime }_1 p_2} {t}_{p^{\prime }_1} \nonumber \\&\quad \quad \!\! \langle \{p_1^{\prime } p_2^{\prime }\} | {V} | \{ p_1 p_2\} \rangle \nonumber \\&\quad = \frac{1}{2} \big ( {V}_{p^{\prime }_1 p_2^{\prime }, p_1 p_2} - {V}_{p^{\prime }_1 p_2^{\prime }, p_2 p_1} - {V}_{p^{\prime }_2 p_1^{\prime }, p_1 p_2} + {V}_{p^{\prime }_2 p_1^{\prime }, p_2 p_1} \big )\nonumber \\&\quad = 2{V}_{p_1^{\prime }p^{\prime }_2, p_1 p_2 }. \end{aligned}$$

(A.7)

In the second line of Eq. (A.7), we have exploited the fact that the kinetic operator is diagonal in the momentum basis. The Schrödinger equation,

$$\begin{aligned} \begin{aligned} {H} | \varPsi \rangle = {E} | \varPsi \rangle , \end{aligned} \end{aligned}$$

(A.8)

in the two-body basis Eq. (A.2) then reads

$$\begin{aligned} \begin{aligned} \sum _{p_1 < p_2} \langle \{ p^{\prime }_1 p^{\prime }_2 \} | {H} | \{ p_1 p_2\} \rangle \langle \{p_1 p_2\} | \varPsi \rangle = {E} \underbrace{\langle \{p_1^{\prime } p_2^{\prime }\} | \varPsi \rangle }_{\equiv \varPsi (p_1^{\prime }p_2^{\prime }) }. \end{aligned} \end{aligned}$$

(A.9)

Here, it will be sufficient to consider only those components of $\varPsi (p_1^{\prime } p_2^{\prime })$ with $p_1^{\prime } < p_2^{\prime }$. Since the basis states are antisymmetric, the other components of $\varPsi (p_1^{\prime } p_2^{\prime })$ with $p_1^{\prime } > p_2^{\prime }$ will differ from the ones with $p_1^{\prime } < p_2^{\prime }$ by a simple phase factor. Plugging Eq. (A.7) into Eq. (A.9) and using $p_1^{\prime } < p_2^{\prime }$, one arrives at a general two-body Schrödinger equation

$$\begin{aligned}&{t}_{p_1^{\prime }} \varPsi (p_1^{\prime } p_2^{\prime }) + {t}_{p_2^{\prime }} \varPsi (p_1^{\prime } p_2^{\prime }) + \sum _{p_1 < p_2 } 2 {V}_{p_1^{\prime } p_2^{\prime }, p_1 p_2} \varPsi (p_1 p_2) \nonumber \\&\quad ={E} \varPsi (p_1^{\prime } p_2^{\prime }). \end{aligned}$$

(A.10)

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., $\lambda _1 = \lambda _2$. In that case, we use $\sum _{p_1<p_2} \rightarrow 1/2 \ \sum _{p_1,p_2}$ and equation Eq. (A.10) becomes

$$\begin{aligned}&{t}_{p_1^{\prime }} \varPsi (p_1^{\prime } p_2^{\prime }) + {t}_{p_2^{\prime }} \varPsi (p_1^{\prime } p_2^{\prime }) + \sum _{p_1, p_2 } {V}_{p_1^{\prime } p_2^{\prime }, p_1 p_2} \varPsi (p_1 p_2) \nonumber \\&\quad ={E} \varPsi (p_1^{\prime } p_2^{\prime }). \end{aligned}$$

(A.11)

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., $|\{\varLambda \varLambda \} \rangle$ and $|\{N \varXi \}\rangle$. Then, the completeness relation is obtained by combining Eqs. (A.3) and (A.6)

$$\begin{aligned}&\int d^3{{\tilde{p}}_1} d^3 {\tilde{p}}_2 \,\Big \{ |\{\varXi {\tilde{p}}_1 N {\tilde{p}}_2 \} \rangle \langle \{\varXi {\tilde{p}}_1 N {\tilde{p}}_2 \} | \nonumber \\&\qquad + \frac{1}{\sqrt{2}}|\{\varLambda {\tilde{p}}_1 \varLambda {\tilde{p}}_2 \} \rangle \langle \{\varLambda {\tilde{p}}_1 \varLambda {\tilde{p}}_2 \} | \frac{1}{\sqrt{2}} \Big \} = \mathbb {1}, \end{aligned}$$

(A.12)

leading to the following expression for the norm of the wave function

$$\begin{aligned} \! \langle \varPsi | \varPsi \rangle = \!\int \!\!d^3{{\tilde{p}}_1} d^3 {\tilde{p}}_2 \Big \{ |\varPsi _{N \varXi }({\tilde{p}}{_1} {\tilde{p}}_2)|^2 + |\frac{1}{\sqrt{2}}\varPsi _{\varLambda \varLambda }({\tilde{p}}_1 {\tilde{p}}_2)|^2 \Big \}. \end{aligned}$$

(A.13)

Therefore, we absorb the $\frac{1}{\sqrt{2}}$-factor into the amplitude of states by introducing a new set of the wave-function components,

$$\begin{aligned} \varPhi _{\varLambda \varLambda }({\tilde{p}}_1 {\tilde{p}}_2)\, {=}\, \frac{1}{\sqrt{2}} \varPsi _{\varLambda \varLambda }({\tilde{p}}_1 {\tilde{p}}_2); \,\,\, \varPhi _{N \varXi }({\tilde{p}}_1 {\tilde{p}}_2) \,{=}\, \varPsi _{N \varXi }({\tilde{p}}_1 {\tilde{p}}_2), \end{aligned}$$

(A.14)

so that the Schrödinger equation Eqs. (A.10,A.11) for the two newly defined components possesses a symmetric form

$$\begin{aligned}&\left( \begin{array}{c c} 2 {t}_{\varLambda } + {V}_{\varLambda \varLambda , \varLambda \varLambda } &{} \sqrt{2} {V}_{\varLambda \varLambda , N \varXi } \\ \sqrt{2} {V}_{N \varXi , \varLambda \varLambda } &{} \quad {t}_{\varXi } + {t}_{N} + 2 {V}_{ N \varXi , N \varXi } \end{array} \right) \left( \begin{array}{c} \varPhi _{\varLambda \varLambda } \\ \varPhi _{N \varXi } \end{array} \right) \nonumber \\&\quad = E \left( \begin{array}{c} \varPhi _{\varLambda \varLambda } \\ \varPhi _{N \varXi } \end{array} \right) , \end{aligned}$$

(A.15)

where, for readability, we have omitted the dependence on ${\tilde{p}}$ and ${\tilde{p}}^{\prime }$. Similarly, for the case where the basis consists of four states $\{|\varLambda \varLambda \} \rangle$, $|\{\varSigma \varSigma \}\rangle$, $|\{ \varLambda \varSigma \}\rangle$ and $|\{ \varXi N \}\rangle$, one analogously defines a new set of wave-function components

$$\begin{aligned}&\varPhi _{\varLambda \varLambda } = \frac{1}{\sqrt{2}} \varPsi _{\varLambda \varLambda }\,; \quad \varPhi _{\varSigma \varSigma } = \frac{1}{\sqrt{2}} \varPsi _{\varSigma \varSigma }\,;\nonumber \\&\varPhi _{\varLambda \varSigma } = \varPsi _{\varLambda \varSigma }\,; \qquad \varPhi _{ N \varXi } = \varPsi _{N \varXi }\,, \end{aligned}$$

(A.16)

for which the Schrödinger equation again possesses a symmetric form

$$\begin{aligned} {H} \left( \begin{array}{c} \varPhi _{\varLambda \varLambda } \\ \varPhi _{\varSigma \varSigma } \\ \varPhi _{\varLambda \varSigma } \\ \varPhi _{N \varXi } \end{array} \right) = E \left( \begin{array}{c} \varPhi _{\varLambda \varLambda } \\ \varPhi _{\varSigma \varSigma } \\ \varPhi _{\varLambda \varSigma } \\ \varPhi _{N \varXi } \end{array} \right) , \end{aligned}$$

(A.17)

with

$$\begin{aligned} {H}= \left( \begin{array}{c c c c } 2 {t}_{\varLambda } + {V}_{\varLambda \varLambda , \varLambda \varLambda } &{} {V}_{\varLambda \varLambda , \varSigma \varSigma } &{} \sqrt{2} {V}_{\varLambda \varLambda , \varLambda \varSigma } &{} \sqrt{2} {V}_{\varLambda \varLambda ,N \varXi } \\ {V}_{\varSigma \varSigma , \varLambda \varLambda } &{} 2{t}_{\varSigma } + {V}_{\varSigma \varSigma , \varSigma \varSigma } &{} \sqrt{2}{V}_{\varSigma \varSigma , \varLambda \varSigma } &{} \sqrt{2} {V}_{\varSigma \varSigma , N \varXi } \\ \sqrt{2}{V}_{\varLambda \varSigma , \varLambda \varLambda } &{} \sqrt{2} {V}_{\varLambda \varSigma , \varSigma \varSigma } &{} {t}_{\varLambda } + {t}_{\varSigma } + 2{V}_{\varLambda \varSigma , \varLambda \varSigma } &{} 2 {V}_{\varLambda \varSigma ,N \varXi } \\ \sqrt{2} {V}_{N \varXi , \varLambda \varLambda } &{} \sqrt{2} {V}_{N \varXi , \varSigma \varSigma } &{} 2 {V}_{N \varXi , \varLambda \varSigma } &{} {t}_{\varXi } + {t}_{N} + 2 {V}_{N \varXi , N\varXi } \end{array} \right) \end{aligned}$$

(A.18)

One sees that there is a $\sqrt{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 [62] (see, e.g., Eq. (2) of [30]). We therefore denote these initial two-body potentials ${V}_{\lambda _1 \lambda _2,\lambda ^{\prime }_1 \lambda ^{\prime }_2}$ with an appropriate factor of $\sqrt{2}$ or 2 or 1 to be our new potential ${{\tilde{V}}}_{\lambda _1 \lambda _2,\lambda ^{\prime }_1 \lambda ^{\prime }_2}$. Expressing in terms of the new potentials ${\tilde{V}}$, the Hamiltonian Eq. (A.18) now has a more intuitive form

$$\begin{aligned} {H}= \left( \begin{array}{c c c c } 2 {t}_{\varLambda } + {{\tilde{V}}}_{\varLambda \varLambda , \varLambda \varLambda } &{} {{\tilde{V}}}_{\varLambda \varLambda , \varSigma \varSigma } &{} {{\tilde{V}}}_{\varLambda \varLambda , \varLambda \varSigma } &{} {{\tilde{V}}}_{\varLambda \varLambda ,N \varXi } \\ {{\tilde{V}}}_{\varSigma \varSigma , \varLambda \varLambda } &{} 2{t}_{\varSigma } + {{\tilde{V}}}_{\varSigma \varSigma , \varSigma \varSigma } &{} {{\tilde{V}}}_{\varSigma \varSigma , \varLambda \varSigma } &{} {{\tilde{V}}}_{\varSigma \varSigma , N \varXi } \\ {{\tilde{V}}}_{\varLambda \varSigma , \varLambda \varLambda } &{} {{\tilde{V}}}_{\varLambda \varSigma , \varSigma \varSigma } &{} {t}_{\varLambda } + {t}_{\varSigma } + {{\tilde{V}}}_{\varLambda \varSigma , \varLambda \varSigma } &{} {{\tilde{V}}}_{\varLambda \varSigma ,N \varXi } \\ {{\tilde{V}}}_{N \varXi , \varLambda \varLambda } &{} {{\tilde{V}}}_{N \varXi , \varSigma \varSigma } &{} {{\tilde{V}}}_{N \varXi , \varLambda \varSigma } &{} {t}_{\varXi } + {t}_{N} + {{\tilde{V}}}_{N \varXi , N\varXi } \end{array} \right) . \end{aligned}$$

(A.19)

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 two-body basis, we will be able to determine the corresponding combinatorial factors for the potentials in each strangeness sector.

1.2 Appendix A.2: Three-body Schrödinger equation

We define the ordered three-body basis states in second quantization and its completeness relations as

$$\begin{aligned} |\{p_1 p_2 p_3 \} \rangle\equiv & {} a^{\dagger }_{p_{1}} a^{\dagger }_{p_2} a^{\dagger }_{p_{3}} | 0 \rangle ;\nonumber \\&\sum _{p_1< p_2 < p_3 } |\{p_1 p_2 p_3 \} \rangle \langle \{p_1 p_2 p_3 \} | = \mathbb {1}.\nonumber \\ \end{aligned}$$

(A.20)

The kinetic and potential matrix elements in the basis Eq. (A.20) read

$$\begin{aligned}&\,\langle \{p_1^{\prime } p_2^{\prime } p_3^{\prime }\}| {T} | \{ p_1 p_2 p_3 \} \rangle \nonumber \\&\quad = \, {T}_{p_1^{\prime } p_1} \delta _{p_2^{\prime } p_2 } \delta _{p_3^{\prime } p_3 }- {T}_{p_1^{\prime } p_1} \delta _{p_2^{\prime } p_3 } \delta _{p_3^{\prime } p_2}+ {T}_{p_1^{\prime } p_2} \delta _{p_2^{\prime } p_3 } \delta _{p_3^{\prime } p_1 }\nonumber \\&\qquad - {T}_{p_1^{\prime } p_2} \delta _{p_2^{\prime } p_1 } \delta _{p_3^{\prime } p_3}+ {T}_{p_1^{\prime } p_3} \delta _{p_2^{\prime } p_1 } \delta _{p_3^{\prime } p_2 }- {T}_{p_1^{\prime } p_3} \delta _{p_2^{\prime } p_2 } \delta _{p_3^{\prime } p_1 }\nonumber \\&\qquad + {T}_{p_2^{\prime } p_1} \delta _{p_1^{\prime } p_3 } \delta _{p_3^{\prime } p_2 }- {T}_{p_2^{\prime } p_1} \delta _{p_1^{\prime } p_2 } \delta _{p_3^{\prime }p_3}+ {T}_{p_2^{\prime } p_2} \delta _{p_1^{\prime } p_1 } \delta _{p_3^{\prime } p_3 }\nonumber \\&\qquad - {T}_{p_2^{\prime } p_2} \delta _{p_1^{\prime } p_3 } \delta _{p_3^{\prime } p_1}+ {T}_{p_2^{\prime } p_3} \delta _{p_3^{\prime } p_1 } \delta _{p_1^{\prime } p_2 }- {T}_{p_2^{\prime } p_3} \delta _{p_1^{\prime } p_1 } \delta _{p_3^{\prime } p_2 }\nonumber \\&\qquad + {T}_{p_3^{\prime } p_1} \delta _{p_1^{\prime } p_2 } \delta _{p_2^{\prime } p_3 }- {T}_{p_3^{\prime } p_1} \delta _{p_1^{\prime } p_3 } \delta _{p_2^{\prime } p_2}+ {T}_{p_3^{\prime } p_2} \delta _{p_1^{\prime } p_3 } \delta _{p_2^{\prime } p_1 }\nonumber \\&\qquad - {T}_{p_3^{\prime } p_2} \delta _{p_1^{\prime } p_1 } \delta _{p_2^{\prime } p_3}+ {T}_{p_3^{\prime } p_3} \delta _{p_1^{\prime } p_1 } \delta _{p_2^{\prime } p_2 }- {T}_{p_3^{\prime } p_3} \delta _{p_1^{\prime } p_2 } \delta _{p_2^{\prime } p_1 },\nonumber \\ \end{aligned}$$

(A.21)

and

$$\begin{aligned}&\!\! \langle \{p_1^{\prime } p_2^{\prime } p_3^{\prime }\}| {V} | \{ p_1 p_2 p_3 \} \rangle = 2 \big \{\nonumber \\&\quad {V}_{p^{\prime }_{2} p^{\prime }_3, p_2p_3} \delta _{p^{\prime }_1 p_1} + {V}_{p^{\prime }_{2} p^{\prime }_3, p_3p_1} \delta _{p^{\prime }_1 p_2} +{V}_{p^{\prime }_{2} p^{\prime }_3, p_1p_2} \delta _{p^{\prime }_1 p_3}\nonumber \\&\qquad +{V}_{p^{\prime }_{3} p^{\prime }_1, p_2p_3} \delta _{p^{\prime }_2 p_1} +{V}_{p^{\prime }_{3} p^{\prime }_1, p_3p_1} \delta _{p^{\prime }_2 p_2} + {V}_{p^{\prime }_{3} p^{\prime }_1, p_1p_2} \delta _{p^{\prime }_2 p_3}\nonumber \\&\qquad +{V}_{p^{\prime }_{1} p^{\prime }_2, p_2p_3} \delta _{p^{\prime }_3 p_1} +{V}_{p^{\prime }_{1} p^{\prime }_2 ,p_3p_1} \delta _{p^{\prime }_3 p_2} + {V}_{p^{\prime }_{1} p^{\prime }_2, p_1p_2} \delta _{p^{\prime }_3 p_3}\big \}\nonumber .\\ \end{aligned}$$

(A.22)

Now, projecting the Schrödinger equation Eq. (A.8) onto the state $|\{p^{\prime }_1 p^{\prime }_2 p^{\prime }_3\} \rangle$ and then utilizing the completeness relation in Eq. (A.20), one arrives at

$$\begin{aligned}&\sum _{p_1< p_2 <p_3}\! \! \Big \{ \langle \{ p^{\prime }_1 p^{\prime }_2 p^{\prime }_3 \} | {T} | \{ p_1 p_2 p_3 \} \rangle \langle \{p_1 p_2 p_3 \} | \varPsi \rangle \nonumber \\&\qquad + \langle \{ p^{\prime }_1 p^{\prime }_2 p^{\prime }_3 \} | {V} | \{ p_1 p_2 p_3 \} \rangle \langle \{p_1 p_2 p_3 \} | \varPsi \rangle \!\Big \} \!= {E} \underbrace{\langle \{p_1^{\prime } p_2^{\prime } p^{\prime }_3 \}| \varPsi \rangle }_{\equiv \varPsi (p_1^{\prime } p_2^{\prime } p_3^{\prime }) } \nonumber \\ \end{aligned}$$

(A.23)

Similar to the case of a two-body basis, here it will be sufficient to consider only those components of $\varPsi (p_1^{\prime } p_2^{\prime } p_3^{\prime })$ with $p^{\prime }_1< p^{\prime }_2 < p^{\prime }_3$. With this condition, only three of the 18 kinetic terms in Eq. (A.21) survive. Hence, we have

$$\begin{aligned}&\sum _{p_1< p_2 <p_3} \langle \{ p^{\prime }_1 p^{\prime }_2 p^{\prime }_3 \} | {T} | \{ p_1 p_2 p_3 \} \rangle \langle \{p_1 p_2 p_3 \} | \varPsi \rangle \nonumber \\&\quad = ({t}_{p^{\prime }_1} + {t}_{p^{\prime }_2} + {t}_{p^{\prime }_3}) \varPsi (p^{\prime }_1 p^{\prime }_2 p^{\prime }_3). \end{aligned}$$

(A.24)

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

$$\begin{aligned}&\!\!\! \sum _{p_1< p_2<p_3} \!\!\Big \{ {V}_{p^{\prime }_{2}p^{\prime }_{3}, p_2 p_3} \delta _{p^{\prime }_1 p_1} \varPsi (p_1 p_2 p_3) \nonumber \\&\qquad + \underset{p_2 \leftrightarrow p_1}{{V}_{p^{\prime }_{2}p^{\prime }_{3}, p_3 p_1} \delta _{p^{\prime }_1 p_2} \varPsi (p_1 p_2 p_3)} +\underset{p_3 \leftrightarrow p_1}{{V}_{p^{\prime }_{2}p^{\prime }_{3}, p_1 p_2} \delta _{p^{\prime }_1 p_3} \varPsi (p_1 p_2 p_3)}\Big \}\nonumber \\&\quad = \sum _{p_1< p_2<p_3} {V}_{p^{\prime }_{2}p^{\prime }_{3}, p_2 p_3} \delta _{p^{\prime }_1 p_1} \varPsi (p_1 p_2 p_3) \nonumber \\&\qquad + \sum _{p_2<p_1<p_3} {{V}_{p^{\prime }_{2}p^{\prime }_{3}, p_3 p_2} \delta _{p^{\prime }_1 p_1} \varPsi (p_2 p_1 p_3)}\nonumber \\&\qquad +\sum _{p_3< p_2< p_1} \underset{p_2 \leftrightarrow p_3}{{{V}_{p^{\prime }_{2}p^{\prime }_{3}, p_3 p_2} \delta _{p^{\prime }_1 p_1}}} \varPsi (p_3 p_2 p_1)\nonumber \\&\quad = \sum _{p_1< p_2<p_3} {V}_{p^{\prime }_{2}p^{\prime }_{3}, p_2 p_3} \delta _{p^{\prime }_1 p_1} \varPsi (p_1 p_2 p_3) \nonumber \\&\qquad + \sum _{p_2<p_1<p_3} {{V}_{p^{\prime }_{2}p^{\prime }_{3}, p_2 p_3} \delta _{p^{\prime }_1 p_1} \varPsi (p_1 p_2 p_3)} \nonumber \\&\qquad + \sum _{p_2< p_3< p_1} {{V}_{p^{\prime }_{2}p^{\prime }_{3}, p_2 p_3} \delta _{p^{\prime }_1 p_1} \varPsi (p_1 p_2 p_3)}\nonumber \\&\quad = \sum _{p_1} \sum _{p_2< p_3 } {V}_{p^{\prime }_{2}p^{\prime }_{3}, p_2 p_3} \delta _{p^{\prime }_1 p_1} \varPsi (p_1 p_2 p_3)\nonumber \\&\quad = \sum _{p_2 < p_3 } {V}_{p^{\prime }_{2}p^{\prime }_{3}, p_2 p_3} \varPsi (p^{\prime }_1 p_2 p_3). \end{aligned}$$

(A.25)

Analogously, the next three terms in Eq. (A.22) yield

$$\begin{aligned}&\!\!\! \sum _{p_1< p_2<p_3} \!\! \Big \{ \underset{p_1 \leftrightarrow p_2}{{V}_{p^{\prime }_{3}p^{\prime }_{1}, p_2 p_3} \delta _{p^{\prime }_2 p_1} \varPsi (p_1 p_2 p_3)} \nonumber \\&\qquad + \,\, {{V}_{p^{\prime }_{3}p^{\prime }_{1}, p_3 p_1} \delta _{p^{\prime }_2 p_2} \varPsi (p_1 p_2 p_3)} +\underset{p_2 \leftrightarrow p_3}{{V}_{p^{\prime }_{3}p^{\prime }_{1}, p_1 p_2} \delta _{p^{\prime }_2 p_3} \varPsi (p_1 p_2 p_3)}\Big \}\nonumber \\&\quad \,\, = \cdots = \sum _{p_1 < p_3 } {V}_{p^{\prime }_{1}p^{\prime }_{3}, p_1 p_3} \varPsi (p_1 p^{\prime }_2 p_3), \end{aligned}$$

(A.26)

and, the three remaining terms result in

$$\begin{aligned}&\sum _{p_1< p_2<p_3} \!\! \! \Big \{ \underset{p_1 \leftrightarrow p_3}{{V}_{p^{\prime }_{1}p^{\prime }_{2}, p_2 p_3} \delta _{p^{\prime }_3 p_1} \varPsi (p_1 p_2 p_3)} \nonumber \\&\qquad + \underset{p_2 \leftrightarrow p_3}{{V}_{p^{\prime }_{1}p^{\prime }_{2}, p_3 p_1} \delta _{p^{\prime }_3 p_2} \varPsi (p_1 p_2 p_3)} + {{V}_{p^{\prime }_{1}p^{\prime }_{2}, p_1 p_2} \delta _{p^{\prime }_3 p_3} \varPsi (p_1 p_2 p_3)}\nonumber \\&\quad = \cdots = \sum _{p_1 < p_2 } {V}_{p^{\prime }_{1}p^{\prime }_{2}, p_1 p_2} \varPsi (p_1 p_2 p^{\prime }_3). \end{aligned}$$

(A.27)

Taking into account Eqs. (A.24-A.27), the Schrödinger equation Eq. (A.23) in the three-body basis Eq. (A.20) can be written as

$$\begin{aligned}&\!\! \big ({t}_{p^{\prime }_1 } + {t}_{p^{\prime }_2 } + {t}_{ p^{\prime }_3}\big ) \varPsi ({p^{\prime }_1 p^{\prime }_2 p^{\prime }_3}) + \sum _{p_1< p_2} \!\!2 {V}_{{p^{\prime }_1 p^{\prime }_2, p_1 p_2}} \varPsi ({p_1 p_2 p^{\prime }_3})\nonumber \\&\qquad + \sum _{p_1< p_3} \!\! 2 {V}_{{p^{\prime }_1 p^{\prime }_3, p_1 p_3}} \varPsi ({p_1 p^{\prime }_2 p_3}) + \!\!\sum _{p_2 < p_3} \!\! 2 {V}_{{p^{\prime }_2 p^{\prime }_3, p_2 p_3}} \varPsi ({p^{\prime }_1 p_2p_3})\nonumber \\&\quad \,\, = E \varPsi ({p^{\prime }_1 p^{\prime }_2 p^{\prime }_3}), \end{aligned}$$

(A.28)

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 states ${\tilde{p}}_1, {\tilde{p}}_2$ etc. For illustration purposes, let us consider Eq. (A.28) in an explicit basis consisting of four states, $|\{N \varLambda \varLambda \}\rangle$, $|\{N \varSigma \varSigma \}\rangle$, $|\{N \varLambda \varSigma \}\rangle$ and $|\{N N \varXi \}\rangle$. The norm of the wave function in this four-particle-state basis can be calculated as follows

$$\begin{aligned} \langle \varPsi | \varPsi \rangle&= \int d^3{{\tilde{p}}_1} d^3 {\tilde{p}}_2 d^3{{\tilde{p}}_3} \Big \{ |\frac{1}{\sqrt{2}}\varPsi _{N \varLambda \varLambda }({\tilde{p}}{_1} {\tilde{p}}_2 {\tilde{p}}_3)|^2\nonumber \\&\quad + \,|\frac{1}{\sqrt{2}}\varPsi _{N \varSigma \varSigma }({\tilde{p}}_1 {\tilde{p}}_2 {\tilde{p}}_3)|^2 +\, \, |\varPsi _{N \varLambda \varSigma }({\tilde{p}}_1 {\tilde{p}}_2 {\tilde{p}}_3)|^2\nonumber \\&\quad + \, \, |\frac{1}{\sqrt{2}}\varPsi _{N N \varXi }({\tilde{p}}_1 {\tilde{p}}_2 {\tilde{p}}_3)|^2\Big \}. \end{aligned}$$

(A.29)

Based on Eq. (A.29), we define a new set of wavefunction components

$$\begin{aligned}&\varPhi _{N \varLambda \varLambda } = \frac{1}{\sqrt{2}} \varPsi _{N \varLambda \varLambda }\,; \quad \varPhi _{N \varSigma \varSigma } = \frac{1}{\sqrt{2}} \varPsi _{N \varSigma \varSigma }\,; \quad \nonumber \\&\varPhi _{N \varLambda \varSigma } = \varPsi _{N \varLambda \varSigma }\,; \quad \varPhi _{N N \varXi } = \frac{1}{\sqrt{2}}\varPsi _{N N \varXi }.\nonumber \\ \end{aligned}$$

(A.30)

The Schrödinger equation Eq. (A.28), applying the wave function components in Eq. (A.30), now has a symmetric form,

$$\begin{aligned}&{ (T + V )} \left( \begin{array}{c} \ \varPhi _{N \varLambda \varLambda } \\ \varPhi _{N\varSigma \varSigma } \\ \varPhi _{N \varLambda \varSigma } \\ \varPhi _{ N N \varXi } \end{array} \right) = E \left( \begin{array}{c} \varPhi _{N \varLambda \varLambda } \\ \varPhi _{N \varSigma \varSigma } \\ \varPhi _{N \varLambda \varSigma } \\ \varPhi _{N N \varXi } \end{array} \right) , \end{aligned}$$

(A.31)

with T being a diagonal matrix

$$\begin{aligned} \!\! {T}= \left( \begin{array}{c c c c } 2 {t}_{\varLambda } + {t}_N &{} 0 &{} 0 &{} 0 \\ 0 &{} 2{t}_{ \varSigma } + {t}_{N} &{} 0 &{} 0 \\ 0 &{} 0 &{} {t}_{\varLambda } + {t}_{\varSigma } +{t}_N &{} 0 \\ 0 &{} 0 &{} 0 &{} 2 {t}_{N} + {t}_{\varXi } \end{array} \right) , \end{aligned}$$

(A.32)

and the symmetric potential matrix

$$\begin{aligned} {V}= \left( \begin{array}{c c c c } 2 {{\tilde{V}}}_{N \varLambda , N \varLambda } + {{\tilde{V}}}_{\varLambda \varLambda , \varLambda \varLambda } &{} {{\tilde{V}}}_{\varLambda \varLambda , \varSigma \varSigma } &{} {{\tilde{V}}}_{\varLambda \varLambda , \varLambda \varSigma } + \sqrt{2} {{\tilde{V}}}_{N \varLambda , N \varSigma } &{} \sqrt{2}{{\tilde{V}}}_{\varLambda \varLambda ,N \varXi } \\ {{\tilde{V}}}_{\varSigma \varSigma , \varLambda \varLambda } &{} 2{{\tilde{V}}}_{N \varSigma , N \varSigma } + {{\tilde{V}}}_{\varSigma \varSigma , \varSigma \varSigma } &{} {{\tilde{V}}}_{\varSigma \varSigma , \varLambda \varSigma } + \sqrt{2}{{\tilde{V}}}_{N \varSigma , N \varLambda } &{} \sqrt{2}{{\tilde{V}}}_{\varSigma \varSigma , N \varXi } \\ {{\tilde{V}}}_{\varLambda \varSigma , \varLambda \varLambda } + \sqrt{2} {{\tilde{V}}}_{N \varSigma , N \varLambda } &{} {{\tilde{V}}}_{\varLambda \varSigma , \varSigma \varSigma } + \sqrt{2}{{\tilde{V}}}_{N \varLambda , N \varSigma } &{} {{\tilde{V}}}_{N \varLambda , N \varLambda } + {{\tilde{V}}} _{N\varSigma , N \varSigma } + {{\tilde{V}}}_{\varLambda \varSigma , \varLambda \varSigma } &{} \sqrt{2} {{\tilde{V}}}_{\varLambda \varSigma , N \varXi } \\ \sqrt{2} {{\tilde{V}}}_{N \varXi , \varLambda \varLambda } &{} \sqrt{2} {{\tilde{V}}}_{N \varXi , \varSigma \varSigma } &{} \sqrt{2} {{\tilde{V}}}_{N \varXi , \varLambda \varSigma } &{} {{\tilde{V}}}_{NN,NN}+ 2 \tilde{{V}}_{N \varXi , N \varXi } \end{array} \right) .\end{aligned}$$

(A.33)

In the last step, we have expressed the potential matrix elements in terms of ${{\tilde{V}}}$ 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 the following, we want to generalize this result to an A-body system.

Table 3 Combinatorial factors of the two-body YN interactions embedded in the A-body space with strangeness $S=-2$

Table 4 Combinatorial factors of the two-body YY interactions embedded in the A-body space with strangeness $S=-2$

1.3 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 $\langle \{p_1' \ldots p_A' \} | V |\{p_1 \ldots p_A \} \rangle$ of the n-particle interaction

$$\begin{aligned} V = \frac{1}{n!} \sum _{\begin{array}{c} k_1,\ldots ,k_n \\ k_1',\ldots ,k_n' \end{array}} V_{ k_1'\ldots k_n', k_1 \ldots k_n } \, a^{\dagger }_{k_1^{\prime }} a^{\dagger }_{k_2^{\prime }} \cdots a^{\dagger }_{k_n^{\prime }} a_{k_n} \cdots a_{k_2} a_{k_{1}} \end{aligned}$$

(A.34)

a total of

$$\begin{aligned} \frac{1}{n!} \, \left( {\begin{array}{c}A\\ n\end{array}}\right) \, \left( {\begin{array}{c}A\\ n\end{array}}\right) \, n! \, n! (A-n)! \end{aligned}$$

(A.35)

different permutations of $V_{ k_1'\ldots k_n', k_1 \ldots k_n }$ contribute. Therein, the first $\frac{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 $\varPsi$ can be written as

$$\begin{aligned} \langle \{p_1' \ldots p_A' \} | \varPsi ' \rangle= & {} \sum _{i_1< i_2 \ldots<i_n} \sum _{ p_{i_1}<\ldots <p_{i_n}} \nonumber \\&\times n! \ V_{ p_{i_1}',\ldots ,p_{i_n}', p_{i_1},\ldots ,p_{i_n} } \nonumber \\&\times \langle \{p_1' \ldots p_{i_1} \ldots p_{i_n} \ldots p_A' \} | \varPsi \rangle \ . \end{aligned}$$

(A.36)

For this form, we assume that $\varPsi '$ is represented using the ordered states $p_{1}'<\ldots <p_{A}'$. Then only one of the $(A-n)!$ different spectator permutations contributes. One of the $\left( {\begin{array}{c}A\\ n\end{array}}\right)$ 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 \ldots <i_n$.

If the interacting particles are (partly) identical, we will again replace the sum over $p_{i_1}<\ldots <p_{i_n}$ 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 $\varPsi$ 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,\ldots , n_p$ particles of each species are $\sqrt{n_1! \ldots 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 $\varPhi$ states. We note that the potential matrix in terms of these states is symmetrical. In summary, the potential matrix elements then reads

$$\begin{aligned} \frac{n!\sqrt{n_1! \ldots n_p!}}{\sqrt{n_1'! \ldots n_p'!}} \ V_{ p_{i_1}',\ldots ,p_{i_n}', p_{i_1},\ldots ,p_{i_n}} \ . \end{aligned}$$

(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 n-particles that contribute identically to Eq. (A.36). The sum over $i_1< i_2 \ldots <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.

As an example, we now consider some selected matrix elements of the $S=-2$ 2-body interaction for $(A-2)$-$\varLambda \varLambda$, $(A-2)$-$\varLambda \varSigma$, $(A-2)$-$\varSigma \varSigma$, and $(A-1)$-$\varXi$ states. For the diagonal matrix elements in particle space, the square root factors in Eq. (A.37) cancel. In this case, the prefactor is just $2 \times$ the number of pairs contributing in the outgoing channel $\times 1/2$ if the active pair in the incoming channel consists of identical particles. Therefore, for our example, we find

$$\begin{aligned}&\langle (A-2)-\varLambda \varLambda | V | (A-2)-\varLambda \varLambda \rangle \nonumber \\&\quad = 2 \cdot 1 \cdot \frac{1}{2} \cdot V_{ p_{\varLambda ,1}'p_{\varLambda ,2}', p_{\varLambda ,1} p_{\varLambda ,2}} = {{\tilde{V}}}_{ p_{\varLambda ,1}' p_{\varLambda ,2}', p_{\varLambda ,1} p_{\varLambda ,2}} \nonumber \\&\langle (A-2)-\varLambda \varSigma | V | (A-2)-\varLambda \varSigma \rangle \nonumber \\&\quad = 2 \cdot 1 \cdot V_{ p_{\varLambda }' p_{\varSigma }', p_{\varLambda } p_{\varSigma }} = {{\tilde{V}}}_{ p_{\varLambda }'p_{\varSigma }', p_{\varLambda } p_{\varSigma }} \nonumber \\&\langle (A-2)-\varSigma \varSigma | V | (A-2)-\varSigma \varSigma \rangle \nonumber \\&\quad = 2 \cdot 1 \cdot \frac{1}{2} \cdot V_{ p_{\varSigma ,1}'p_{\varSigma ,2}', p_{\varSigma ,1} p_{\varSigma ,2}} = {{\tilde{V}}}_{ p_{\varSigma ,1}' p_{\varSigma ,2}', p_{\varSigma ,1} p_{\varSigma ,2}} \nonumber \\&\langle (A-1) - \varXi | V | (A-1) - \varXi \rangle \nonumber \\&\quad = 2 \cdot (A-1) \cdot V_{ p_{N}' p_{\varXi }', p_{N} p_{\varXi }} = ( A - 1 ) {{{\tilde{V}}}}_{ p_{N}'p_{\varXi }', p_{N} p_{\varXi }} \ .\nonumber \\ \end{aligned}$$

(A.38)

In the last step, we have exploited the results of the two-body system that relate V to ${{\tilde{V}}}$ in Eq. (A.19). The resulting combinatorial factors agree with the expectation that the interaction just has to be multiplied by the number of pairs contributing. More interesting is the case of transitions. Here, we first look at transitions between $\varSigma \varSigma$ and $\varLambda \varSigma$ states. Because the identity of the particles changes, we now have additionally the contribution of the square root factors. They are also important to guarantee that the interaction matrix is symmetric.

The result for the two matrix elements is

$$\begin{aligned}&\langle (A-2)-\varSigma \varSigma | V | (A-2)-\varLambda \varSigma \rangle \nonumber \\&\quad = 2 \cdot 1 \cdot 1 \cdot \frac{\sqrt{(A-2)!}}{\sqrt{2!(A-2)!}} \ \cdot V_{ p_{\varSigma ,1}'p_{\varSigma ,2}', p_{\varLambda } p_{\varSigma }} \nonumber \\&\quad = {{\tilde{V}}}_{ p_{\varSigma ,1}' p_{\varSigma ,2}', p_{\varLambda } p_{\varSigma }} \nonumber \\&\langle (A-2)-\varLambda \varSigma | V | (A-2)-\varSigma \varSigma \rangle \nonumber \\&\quad = 2 \cdot 1 \cdot \frac{1}{2} \cdot \frac{\sqrt{2!(A-2)!}}{\sqrt{(A-2)!}} \ \cdot V_{ p_{\varLambda }'p_{\varSigma }', p_{\varSigma ,1} p_{\varSigma ,2}} \nonumber \\&\quad = {{\tilde{V}}}_{ p_{\varLambda }' p_{\varSigma }', p_{\varSigma ,1} p_{\varSigma ,2}} \end{aligned}$$

(A.39)

and just reflects the number of YY pairs in the A-body state. The final examples are transitions of $\varLambda \varLambda$ and $N \varXi$. For these matrix elements, several pairs contribute and the identity of the particles changes. It is reassuring that we also find in this case symmetry of the potential matrix elements

$$\begin{aligned}&\langle (A-2)-\varLambda \varLambda | V | (A-1)- \varXi \rangle \nonumber \\&\quad = 2 \cdot 1 \cdot 1 \cdot \frac{\sqrt{(A-1)!}}{\sqrt{2!(A-2)!}} \cdot V_{ p_{\varLambda ,1}'p_{\varLambda ,2}', p_{N} p_{\varXi }} \nonumber \\&\quad = \sqrt{A-1} \ {{\tilde{V}}}_{ p_{\varLambda ,1}' p_{\varLambda ,2}', p_{N} p_{\varXi }} \nonumber \\&\langle (A-1)-\varXi | V | (A-2)- \varLambda \varLambda \rangle \nonumber \\&\quad = 2 \cdot ( A - 1 ) \cdot \frac{1}{2} \cdot \frac{\sqrt{2!(A-2)!}}{\sqrt{(A-1)!}} \cdot V_{ p_{N}'p_{\varXi }', p_{\varLambda ,1} p_{\varLambda ,2}} \nonumber \\&\quad = \sqrt{A-1} \ {{\tilde{V}}}_{ p_{N}' p_{\varXi }', p_{\varLambda ,1} p_{\varLambda ,2}} \ . \end{aligned}$$

(A.40)

In this way, it is straightforward to identify all relevant combinatorial factors for our calculations. For the $S=0$ interactions, where particle transitions do not occur, it is simply given by the number of NN pairs in the state, i.e. $\left( {\begin{array}{c}A-2\\ 2\end{array}}\right)$ and $\left( {\begin{array}{c}A-1\\ 2\end{array}}\right)$ for $(A-2)$-YY and $(A-1)$-$\varXi$ states, respectively. The factors for $S=-1$ ($S=-2$) interactions are summarized in Table 3 (Table 4). To shorten the presentation, we only include particle transitions in one direction. The other one is given by the symmetry of the potentials.