Jacobi no-core shell model for p-shell hypernuclei

Abstract

We extend the recently developed Jacobi no-core shell model to hypernuclei. Based on the coefficients of fractional parentage for ordinary nuclei, we define a basis where the hyperon is the spectator particle. We then formulate transition coefficients to states that single out a hyperon–nucleon pair which allow us to implement a hypernuclear many-baryon Hamiltonian for p-shell hypernuclei. As a first application, we use the basis states and the transition coefficients to calculate the ground states of $^{4}_{\varLambda }\hbox {He}$, $^{4}_{\varLambda }\hbox {H}$, $^{5}_{\varLambda }\hbox {He}$, $^{6}_{\varLambda }\hbox {He}$, $^{6}_{\varLambda }\hbox {Li}$, and $^{7}_{\varLambda }\hbox {Li}$ and, additionally, the first excited states of $^{4}_{\varLambda }\hbox {He}$, $^{4}_{\varLambda }\hbox {H}$, and $^{7}_{\varLambda }\hbox {Li}$. In order to obtain converged results, we employ the similarity renormalization group (SRG) to soften the nucleon–nucleon and hyperon-nucleon interactions. Although the dependence on this evolution of the Hamiltonian is significant, we show that a strong correlation of the results can be used to identify preferred SRG parameters. This allows for meaningful predictions of hypernuclear binding and excitation energies. The transition coefficients will be made publicly available as HDF5 data files.

1 Introduction

After more than 65 years of research on hypernuclei, our knowledge of the interaction of hyperons with nucleons or with other hyperons still remains on a modest level. This situation is rather unsatisfactory given the important role hyperons play for various aspects of nuclear physics as well as for astrophysics [1,2,3,4,5]. For example, as extensively discussed in recent years, the hyperon interaction could have a significant impact on the properties of neutron stars [3,4,5]. The reason for the large uncertainty is the tremendous difficulty to perform scattering experiments involving hyperons and the fact that no two-baryon bound state has been found so far, except for the well known deuteron. An important source of information has been the spectroscopy of hypernuclei [6]. New experiments are planned at facilities like J-PARC, FAIR, MAMI and JLab [7,8,9,10,11,12], some to study the scattering of hyperons on nucleons, but mostly measurements of bound states of ordinary nuclei with hyperons. Such new and very probably more precise data will not only be phenomenologically interesting, but also enable us to explore the underlying interactions in more detail. The latter is now possible because even fairly complex systems can be treated theoretically on a microscopic level, thanks to improved algorithms and increasing computational resources. Indeed, nowadays, one can solve the Schrödinger equation for hypernuclei up to the p-shell based on realistic and rather elaborate baryon-baryon interactions [13,14,15]. Thus, it has become feasible to study detailed features of the baryonic forces, like the spin-dependence of hypernuclear interactions, which are inaccessible in direct scattering experiments. With these theoretical advances, the new data on hypernuclei will definitely provide valuable input to pin down the underlying interactions. Eventually, the hypernuclear data could be directly utilized in fits of interaction parameters.

However, a direct use of hypernuclear data requires solving the hypernuclear many-body problem many times and, therefore, calls for a very efficient calculation scheme. Several methods have been employed in the past to study hypernuclei. For local interactions, configuration space methods, e.g. hyperspherical harmonics, Green’s function Monte Carlo, expansion in Gaussians or stochastic variational method (SVM), have been successfully used to predict properties of light hypernuclei [16,17,18,19,20]. For very light systems, that goal can be likewise achieved by solving the Faddeev- or Yakubovsky equations in momentum space [15, 21,22,23,24,25]. Those methods allow one also to deal with non-local two-body interactions, but it is difficult to extend the approaches to larger systems. Alternatively, shell model calculations have been a quite successful tool to understand properties of hypernuclei, in particular the energy level splittings [26,27,28,29]. However, that approach requires specific effective interactions that are not easily related to free-space baryon-baryon interactions. The same disadvantage also holds for density functional approaches, which have been applied to rather complex hypernuclei [30, 31]. Recently, nuclear lattice effective field theory (NLEFT) has been extended to hypernuclei using the impurity lattice Monte Carlo technique [32]. Although this first study has been performed with somewhat simplified (spin-independent) interactions, that method promises the application of free-space interactions up to medium-heavy hypernuclei.

One specifically interesting approach to tackle bound baryon systems is the no-core shell model (NCSM) [33]. An essential tool is here the representation in terms of a harmonic oscillator (HO) basis. There are several variants of the approach. In most applications so far, a single-particle Slater-determinant basis has been chosen. This realization has been very successfully employed for studying ordinary nuclei and even hypernuclei [14, 34,35,36], especially, when the so-called importance truncation is implemented [14, 35, 36]. Highly accurate results for binding energies, excitation energies and even radii have been obtained. Generally, the problem becomes very high dimensional, not least because the center-of-mass (CM) motion cannot be separated off and because angular momentum and isospin conservation cannot be exploited to limit the basis size.

Such a complication can be avoided by using a Jacobi relative coordinate basis. This, however, requires a very tedious antisymmetrization for the nucleonic states [37, 38]. Nevertheless, the method can be advantageous when many calculations are required for variations of the underlying interactions, e.g. in fitting procedures, since the antisymmetrization and other preparatory steps can be accomplished independently of the interactions. The final step of the calculation itself can then be much more efficiently performed than in the standard NCSM so that it becomes feasible to solve the problem hundreds or even thousands of times or with limited computational resources. The work of Gazda et al. [14, 34] has already been employing this Jacobi NCSM (J-NCSM) for s-shell hypernuclei. It is the main aim of the present work to extend the J-NCSM approach to p-shell hypernuclei. The new approach is then used to study in more detail the $^{4}_{\varLambda }\hbox {He}$, $^{5}_{\varLambda }\hbox {He}$, $^{6}_{\varLambda }\hbox {Li}$ and $^{7}_{\varLambda }\hbox {Li}$ systems based on the next-to-leading order (NLO) hyperon-nucleon (YN) interaction derived within chiral effective field theory (EFT) [25, 39, 40]. For interactions from chiral EFT, it is possible to obtain reliable uncertainty estimates of the results [15, 25], utilizing different orders of the chiral expansion and/or by exploiting the regulator (cutoff) dependence of these interactions (where the latter method provides only a lower limit for the error). For ordinary nuclei, such estimates are now regularly performed [41, 42].

As usual, the NCSM requires a further softening of the nucleon-nucleon (NN) and YN interactions. To this aim, we apply the similarity renormalization group (SRG) to the NN and YN potentials [43, 44]. This method has the advantage that an effective interaction can be systematically derived from the starting NN and YN interactions, which can then be equally well employed in momentum space and HO space. The SRG evolution gives rise also to so-called induced three-body and many-body forces. In the present study, we will not take into account such induced many-body forces (for the application of the SRG induced YNN forces see [35, 36, 45]). Therefore, a part of this work is devoted to study the SRG dependence of the binding energies, excitation energies and $\varLambda$-separation energies.

In Sect. 2, we start with a definition of our basis states based on the totally antisymmetrized nucleonic states defined in [38]. Practical calculations can only be performed when the transition matrix elements to states that single out NN or YN pairs are known. The calculation of these matrix elements is explained in detail in Sect. 3. This already concludes the description of the Jacobi NCSM. As mentioned above, for explicit calculations, we, however, also need soft interactions. In Sect. 4, we therefore discuss the basic features of chiral interactions and their SRG evolvement including the impact on the binding energy for $^{3}_{\varLambda }\hbox {H}$ when the SRG-induced three-baryon force (3BF) is neglected. For this study, we will make use of solutions based on the Faddeev equations. The application of the Jacobi NCSM then follows in Sect. 5. We first present a detailed benchmark for $^{4}_{\varLambda }\hbox {H}$/$^{4}_{\varLambda }\hbox {He}$ to Yakubovsky results and then continue towards A = 5–7 hypernuclei. Our conclusions are finally given in Sect. 6. Some technicalities are relegated to the appendices.

2 NCSM basis in Jacobi coordinates

The translationally invariant many-body Hamiltonian of a system consisting of $(A-1)$ nucleons and a single-strangeness hyperon Y ($Y=\varLambda$ or $\varSigma$) in Jacobi relative coordinates can be written as follows

$$\begin{aligned} H= & {} H^{S=0} + H^{S=-1} \nonumber \\= & {} \sum _{i < j=1}^{A-1} \Big ( \frac{2p^{2}_{ij}}{M(t_Y)} + V^{NN}_{ij} \Big )\nonumber \\&+ \sum _{i=1}^{A-1} \Big ( \frac{m_N + m(t_Y)}{M(t_Y)} \frac{p^2_{iY}}{2\mu _{NY}} + V^{YN}_{iY}\nonumber \\&+\frac{1}{A-1} \big (m(t_Y) - m_{\varLambda }\big ) \Big ). \end{aligned}$$

(1)

Here, $m_N$, $m(t_Y)$ and $\mu _{NY}$ are nucleon-, hyperon-, and their reduced masses, respectively, which we define by $m_N =2 m_n m_p/(m_n + m_p)$, $m(t_Y =0) = m_{\varLambda }$, and $m(t_Y =1) = (m_{\varSigma ^+} + m_{\varSigma ^{-}} + m_{\varSigma ^{0}})/3$. For simplicity, we assume isospin symmetry. A generalization to unequal masses within the isospin multiplet of nucleons and of $\varSigma$’s is straightforward but will not be considered here. The total rest mass of the system, $M(t_Y)=(A-1)m_N + m(t_Y)$, depends explicitly on the hyperon isospin $t_Y$ because an explicit $\varLambda$-$\varSigma$ conversion is allowed. The term $m(t_Y) - m_{\varLambda }$ then accounts for the difference in the rest masses of the two hyperons. The relative Jacobi momenta of NN and YN pairs,

$$\begin{aligned} p_{ij} = \frac{1}{2}(k_i - k_j), \end{aligned}$$

(2)

and

$$\begin{aligned} p_{iY } = \frac{m(t_Y)}{m_N + m(t_Y)} k_{i} - \frac{m_N}{m_N + m(t_Y)} k_Y \end{aligned}$$

(3)

are linear combinations of the momenta $k_i$ and $k_Y$ of the i-th nucleon and the hyperon, respectively. $V^{NN}_{ij}$ and $V^{YN}_{iY}$ are the corresponding NN and YN potentials.

Since hyperons ($\varLambda$, $\varSigma$) and nucleons are distinguishable, hypernuclear basis functions, denoted as $|\alpha ^{* (Y)}\rangle$, can be formed by coupling the hyperon HO states $|Y\rangle$, which describe the relative motion of a single hyperon Y with respect to the CM of the $(A-1)$N core, to the fully antisymmetrized states of the core $|\alpha _{(A-1)N}\rangle$

(4)

where $\alpha _{(A-1)N}$ stands for a complete set of all necessary quantum numbers characterizing the fully antisymmetrized states of an $(A-1)$N system: the total HO energy quantum number ${\mathcal {N}}_{{A-1}}$, total angular momentum ${J}_{{A-1}}$, isospin ${T}_{{A-1}}$, and the state indices $\zeta _{A-1}$ (that distinguish different $|\alpha _{(A-1)N}\rangle$ states with the same set of ${\mathcal {N}}_{A-1}, J_{A-1}$ and $T_{A-1}$). These antisymmetrized states for $A \ge 4$ systems are computed iteratively starting from the naturally antisymmetrized basis for two nucleons, for more detail we refer to Ref. [38]. Here, the superscript $(*Y)$ represents the separation of the hyperon Y from the $(A-1)$N core. The hyperon states $|Y\rangle$ are described by a similar set of quantum numbers: the HO energy quanta $n_Y$, the orbital angular momentum $l_Y$ and spin $s_Y$ which are coupled to the relative angular momentum $I_Y$, and the isospin $t_Y$ as well. The last line in Eq. (4) defines the ordering in which the quantum numbers of the two subclusters are combined to form the total angular momentum and total isospin of the system, J and T, respectively, whose values are given by the physical state of interest. Also, for practical realization, the total HO quantum numbers ${\mathcal {N}}$ of the basis states are constrained by the maximum number of the single-particle oscillators ${\mathcal {N}}_{max}$ (also referred to as the model space size), i.e. ${\mathcal {N}} = {\mathcal {N}}_{A-1} +2n_Y +l_Y \le {\mathcal {N}}_{max}$. The state index $\zeta$ that distinguishes different basis states $|\alpha ^{* (Y)}\rangle$ with the same ${\mathcal {N}}, J$ and T is omitted for simplifying the notation. Finally, on the right-hand side of Eq. (4), the graphical representation of the basis is shown. The small red circle denotes a hyperon spectator while the big black circle represents the system of $(A-1)$N.

3 Separation of NN and YN pairs

With the basis defined in Eq. (4), the matrix elements of the Hamiltonian in Eq. (1) now read

$$\begin{aligned} \langle \alpha ^{*(Y)} | H | \alpha ^{\prime *(Y)} \rangle= & {} \, \langle \alpha ^{*(Y)} | H^{S=0} | \alpha ^{\prime *(Y)} \rangle \nonumber \\&+\, \langle \alpha ^{*(Y)} | H^{S=-1} | \alpha ^{ \prime *(Y)} \rangle . \end{aligned}$$

(5)

The basis states $|\alpha ^{*(Y)}\rangle$ are however not suitable for evaluating the $\hbox {H}^{S=0}$ and $H^{S=-1}$ matrix elements as they do not depend explicitly on the relative coordinates of the involved NN or YN pairs. To facilitate the evaluation of Eq. (5), we expand $|\alpha ^{*(Y)}\rangle$ in two additional bases of intermediate states $|\big (\alpha ^{*(2)}\big )^{*(Y)}\rangle$ and $|\alpha ^{*(YN)} \rangle$ that explicitly single out the active NN or a YN pairs, respectively. Also the superscripts represent subsystems that are separated out. Clearly, the former states $|\big (\alpha ^{*(2)}\big )^{*(Y)}\rangle$ are needed for evaluating the first part in Eq. (5) involving $H^{S=0}$, while the latter ones are necessary for the evaluation of the second part that involves $H^{S=-1}$.

The first set of auxiliary states $|\big (\alpha ^{*(2)}\big )^{*(Y)}\rangle$ can be directly constructed by coupling the hyperon states $|Y\rangle$, depending on Jacobi coordinates of Y relative to the CM of (A-1)N, to the $(A-1)$N states that consist of antisymmetrized subclusters of $(A-3)$N and 2N. In the notation of Ref. [38], this reads

(6)

Here, $\alpha ^{*(2)}_{(A-1)}$ stands for the total HO energy quantum number ${\mathcal {N}}_{\alpha ^{*(2)}_{(A-1)}},$ the total angular momentum $J^{*(2)}_{A-1}$, isospin $T^{*(2)}_{A-1}$ and state index $\zeta ^{*(2)}_{A-1}$, as introduced in [38]. Naturally, the total HO energy quantum number $\tilde{{\mathcal {N}}}$ in Eq. (6) is also restricted by $\tilde{{\mathcal {N}}} \le {\mathcal {N}}_{max}$. With the graphical representations of $|\big (\alpha ^{*(2)}\big )^{*(Y)} \rangle$ and $|\alpha ^{*(Y)}\rangle$, one can quickly relate the expansion coefficients $\big \langle \alpha ^{*(Y)} | \big (\alpha ^{*(2)}\big )^{*(Y)} \big \rangle$ to the transition coefficients of the ${(A-1)}$N system

(7)

for which an explicit expressions has been derived in [38, 46]. The Kronecker symbol $\delta _{spectator}$ is to ensure the conservation of the quantum numbers of the hyperon and the $(A-1)$N system,

$$\begin{aligned} \begin{array}{rl} \delta _{spectator} &{}= \delta _{{\mathcal {N}} \mathcal {\tilde{N}}}\delta _{Y} \delta _{core},\\ \delta _{Y } &{} = \delta _{n_Y \tilde{n}_Y} \delta _{l_Y \tilde{l}_Y} \delta _{I_Y \tilde{I}_Y} \delta _{t_Y \tilde{t}_Y},\\ \delta _{core} &{}= \delta _{{\mathcal {N}}_{A-1} {\mathcal {N}}^{*(2)}_{A-1} } \delta _{{J_{A-1}} {J^{*(2)}_{A-1}} } \delta _{{T_{A-1}} {T^{*(2)}_{A-1}} }. \end{array} \end{aligned}$$

Hence, the matrix elements of the nucleonic Hamiltonian $\langle \alpha ^{*(Y)} | H^{S=0} | \alpha ^{\prime *(Y)} \rangle$ now become

(8)

with summations over intermediate states being implied. The remaining unknown term in Eq. (8) is simply the matrix elements of $H^{S=0}$ in the basis of $(A-1)$ nucleons.

Similarly, in order to construct the intermediate states $|\alpha ^{*(YN)} \rangle$, one combines the states describing a YN pair, $|YN \rangle$, with the antisymmetrized basis of an $(A-2)$N system, $|\alpha _{(A-2)}\rangle$

(9)

Again, $|\alpha _{YN}\rangle$ and $|\alpha _{A-2}\rangle$ represent the complete sets of quantum numbers characterizing the states of the two-body hyperon-nucleon and the $(A-2)$N subsystems. Note that, in contrast to two-nucleon states, there is no antisymmetry requirement for $|\alpha _{YN}\rangle$. The relative motion of the $(A-2)$N cluster with respect to the separated out YN pair is specified by the HO energy number $n_{\lambda }$ and the orbital angular momentum ${\lambda }$. For evaluating the overlap $\langle \alpha ^{*(Y)} | \alpha ^{*(YN)}\rangle$, we need to exploit another set of auxiliary states $|\big (\alpha ^{*(1)}\big )^{*(Y)}\rangle$ in which a hyperon and a nucleon are explicitly singled out

(10)

With the help of Eq. (10), the transition coefficients $\langle \alpha ^{*(Y)} | \alpha ^{*(YN)}\rangle$ can be computed in two steps as follows

(11)

Here also an explicit summation over the auxiliary states is assumed. Clearly, the first overlap is essentially given by the coefficients of fractional parentage (cfp) of an $(A-1)$N system, which basically determine the antisymmetrized basis of $(A-1)$ nucleons in terms of the $| \alpha ^{*(1)}_{(A-1)}\rangle$ states [38], and is therefore well known. Hence, only the second transition in Eq. (11) needs to be taken care of. This transition is a transformation between different Jacobi coordinates and therefore given by the general coordinate transformation formula derived in [38]. We skip the detailed derivation but provide the final expression in Appendix A. Finally, a summation over the intermediate states is carried out. Let us again stress that both, the transition coefficients $\langle \alpha ^{*(Y)} | \big (\alpha ^{*(2)}\big )^{*(Y)}\rangle$ and $\langle \alpha ^{*(Y)} | \alpha ^{*(YN)} \rangle$, are independent of the HO frequency (HO-$\omega$) as well as of the interactions employed. They can therefore be prepared in advance and stored in the machine-independent HDF5 format so that the parallel input and output can be performed most efficiently. The corresponding files can be found at [47].

Once the transition coefficients $\langle \alpha ^{*(Y)} | \alpha ^{*(YN)}\rangle$ are known, the single-strangeness Hamiltonian matrix elements $\langle \alpha ^{*(Y)} | H^{S=-1}| \alpha ^{\prime *(Y)}\rangle$ are computed similarly as in Eq. (8):

(12)

Thus, the evaluation of the matrix elements$\langle \alpha ^{*(Y)} | H^{S=0} | \alpha ^{\prime *(Y)} \rangle$ and $\langle \alpha ^{*(Y)} | H^{S=-1} | \alpha ^{\prime *(Y)} \rangle$ can be traced back to multiplications of very large but sparse matrices. As usual, we solve the eigenvalue problem using the Lanczos method so that these matrix multiplications must be computed again and again. Therefore, an efficient method to evaluate such product matrices is extremely important. More details on the technical realization are given in Ref. [48].

4 SRG evolution for chiral NN and YN interactions

We follow the formalism initially applied by Wegner [43] to solid state physics and later employed by Bogner, Furnstahl and Perry [49] to nuclear interactions, which defines the SRG evolution in terms of a unitary transformation depending on a flow parameter s

$$\begin{aligned} H_{s}^{} = U_{s}^{} H_0 U^{\dagger }_s \equiv T_{rel}^{} + V_s^{}. \end{aligned}$$

(13)

Here $H_0 = H_{s=0}$ is the initial (bare) Hamiltonian and $T_{rel}$ is the intrinsic relative kinetic energy operator that also includes the mass difference term when one allows for particle conversions in the Hamiltonian. The parameter s has the unit of energy^-2 and varies continuously from zero to $\infty$. Note that, although the flow equation is solved with respect to s, for characterizing the SRG-evolved potentials, we will utilize a more intuitive variable

$$\begin{aligned} \lambda =\left( \frac{4 \mu ^2}{s}\right) ^{1/4}~, \end{aligned}$$

(14)

with $\mu ={m_N\, m_\varLambda }/{(m_N+m_\varLambda )}$ for YN interactions and $\mu = {m_N}/{2}$ for NN forces. A similar definition for $\lambda$ was introduced in [49]. $\lambda$ can be (to some approximation) identified with the width of the band for which the SRG evolved matrix elements of the interaction are non-zero. By differentiating the transformation Eq. (13), one obtains the evolution equation for the Hamiltonian

$$\begin{aligned} \frac{d H_s}{ds} = \frac{d V_s}{ds} = [\eta _{s}, H_s] \end{aligned}$$

(15)

where the generator

$$\begin{aligned} \eta _{s}^{} = \frac{dU_s^{}}{ds} U^{\dagger }_s = -\eta ^{\dagger }_{s} \end{aligned}$$

(16)

is an anti-hermitian operator. Usually, $\eta _{s}$ is taken as a commutator of an hermitian operator $G_s$ with the Hamiltonian, $\eta _s = [ G_s, H_s]$. The operator $G_s$ is often chosen such that the evolved Hamiltonian $H_s$ possesses a desired form. For our purpose of decoupling the low- and high-momentum components, the simplest, but yet very useful generator, is the relative kinetic energy excluding the mass shift. We take

$$\begin{aligned} G_s = \frac{p^2}{2 \mu } \end{aligned}$$

(17)

with p being the particles relative momentum. The flow equation Eq. (15) now becomes an operator equation

$$\begin{aligned} \frac{d V_s}{ds} = \Big [\Big [\frac{p^2}{2\mu }, V_s\Big ], H_s\Big ] \ . \end{aligned}$$

(18)

This is then solved in a partial-wave relative momentum basis

$$\begin{aligned} | p \,(ls)J; \, t_1 m_{t_1}S_1\, \,t_2 m_{t_2}S_2 \rangle \equiv |p \alpha \rangle , \end{aligned}$$

(19)

where l is the orbital angular momentum that combines with the total spin s to form the total angular momentum J. Further, $(t_i, m_{t_i}, S_i)_{i=1,2}$ are sets of the intrinsic quantum numbers that distinguish different particle states: isospin, isospin projection and strangeness. The normalization of the basis states Eq. (19) simply reads

$$\begin{aligned} \sum _{\alpha } \int dp p^2 \,|p \alpha \rangle \langle p \alpha | = 1. \end{aligned}$$

(20)

After projecting Eq. (18) onto the basis Eq. (19), one obtains the flow equation in form of an integro-differential equation

$$\begin{aligned} \frac{dV^{\alpha \alpha ^{\prime }}_{s}(p p^{\prime }) }{ds}= & {} \Big [ T^{\alpha }_{rel}(p) \frac{p^{\prime 2}}{2\mu ^{\alpha ^{\prime }} } + T^{\alpha ^{\prime }}_{rel}(p^{\prime }) \frac{p^{ 2}}{2\mu ^{\alpha } } \nonumber \\&-T^{\alpha }_{rel}(p) \frac{p^{2}}{2\mu ^{\alpha } } - T^{\alpha ^{\prime }}_{rel}(p^{\prime }) \frac{p^{\prime 2}}{2\mu ^{\alpha ^{\prime }} } \Big ] V^{\alpha \alpha ^{\prime }}_{s}(p p^{\prime })\nonumber \\&+ \sum _{{\tilde{\alpha }}}\int _{0}^{\infty } dk k^2 \Big [ \frac{p^{2}}{2\mu ^{\alpha }} + \frac{p^{\prime 2}}{2\mu ^{\alpha ^{\prime }}} - \frac{k^2}{\mu ^{{\tilde{\alpha }}}} \Big ] \nonumber \\&\times V^{\alpha {\tilde{\alpha }}}_s(pk) V^{{\tilde{\alpha }} \alpha ^{\prime }}_s(k p^{\prime }). \end{aligned}$$

(21)

Here, the reduced mass $\mu$ and $T_{rel}$ depend explicitly on the particle states $\alpha$ since physical masses are employed for the SRG evolution. We solve the flow equation Eq. (21) numerically using a non-equidistant momentum grid characterized by the ultraviolet momentum cutoff $p_{max}$ and N Gauss-Legendre integration points $p_n$ with corresponding weights $w_n (n=1,\cdots N)$. Since the initial potentials often vary at low momenta faster than at high momenta, it is useful to define the grid such that it is sparse at high momenta but denser at the low-momentum region.

Fig. 1

Contour plot of the YN potential matrix elements for all possible particle channels with charge $Q=0$ and in the $^1S_0$ partial wave. The potentials are evolved to four different values of the YN flow parameter: $\lambda _{YN} =98$ fm$^{-1}$ (first column, almost non-evolved), $\lambda _{YN}= 3$ fm$^{-1}$ (second column, slightly evolved), $\lambda _{YN} =1.6$ fm$^{-1}$ (third column) and $\lambda _{YN} =0.868$ fm$^{-1}$ (last column). The initial potential is the YN NLO interaction with a regulator of $\varLambda _{Y}= 650$ MeV

Discretizing the flow equation leads to a set of coupled differential equations which is then solved using the advanced multi-step Adams PECE (Predict Estimation Correct Estimation) method [50]. The SRG-evolution of the YN interaction NLO19 with a regulator of $\varLambda _{Y}=650$ MeV is illustrated in Fig. 1. The contour plots are the potentials for all the particle channels with zero charge and in the $^1S_0$ partial wave. The initial potential NLO19(650) is evolved to four different values of the YN flow parameter: $\lambda _{YN} = 98$ fm$^{-1}$ (almost non-evolved, bare interaction), $\lambda _{YN} = 3$ fm$^{-1}$ (slightly evolved), $\lambda _{YN} = 1.6$ fm$^{-1}$ (commonly used) and the extreme case $\lambda _{YN} = 0.868$ fm$^{-1}$. As expected, the SRG evolution steadily drives the potentials toward a diagonal form decoupling the low- and higher-momentum states. While the bare NLO19 shows a strong repulsive behavior for almost all particle channels over the entire momentum range, the SRG-evolved potentials become slightly attractive at low momenta but remain repulsive at high momenta.

We explicitly checked that NN and YN scattering observables remain unchanged by this unitary transformation. At this point, we neglect induced three-baryon forces (3BFs). In this approximation, the evolution of NN and YN forces is not linked to each other and we can choose $\lambda _{NN}$ and $\lambda _{YN}$ independently.

As a first application, we apply the SRG transformed interactions to obtain binding energies $E(^3_\varLambda \mathrm{H})$ and the $\varLambda$ separation energies $B_\varLambda (^3_\varLambda \mathrm{H}) = E(^2 \mathrm{H}) - E(^3_\varLambda \mathrm{H})$ of $^3_\varLambda \hbox {H}$.

Since the $^3_\varLambda \hbox {H}$ is predominantely a weakly bound $\varLambda$ to a significantly stronger bound deuteron, it is very difficult to obtain converged results for the binding energies using the NCSM. Therefore, for this study, we use solutions based on Faddeev equations (see Appendix B). With this method, an accuracy of 1 keV for these energies is routinely achieved.

Fig. 2

Dependence of $B_\varLambda (^3_\varLambda \mathrm{H})$ on $\lambda _{YN}$ for $\lambda _{NN} = 1.6$ (blue +) and 2.4 $\hbox {fm}^{-1}$ (orange x). Starting point of the NN SRG evolution is the Idaho-$\hbox {N}^3$LO(500) interaction [51]. For YN, the NLO19(600) interaction [25] is used. The black solid horizontal line and cyan band indicates the experimental value [52] and its uncertainty. The blue dashed and orange dash-dotted lines are results for the bare YN interaction and for $\lambda _{NN} = 1.6$ and 2.4 $\hbox {fm}^{-1}$, respectively

In Fig. 2, $B_\varLambda (^3_\varLambda \mathrm{H})$ is shown for one typical choice of the NN and YN starting interactions. It can be seen that the dependence on the flow parameter of the NN interaction is of the order of 20 keV. But, unfortunately, it is also clear that the dependence on $\lambda _{YN}$ is rather significant, indicating a non-negligible contribution of SRG induced three-baryon interactions. We will discuss later in Sect. 5.5 how this issue could be possibly resolved without explicitly taking the induced 3BFs into account. Note that, for $\lambda _{YN} \lesssim 1.0$ $\hbox {fm}^{-1}$, the separation energy is in fair agreement with experiment and the result of calculations based on the bare YN interaction.

5 Results

As first application of the Jacobi NCSM, we employ the approach to investigate some interesting hypernuclear systems up to the p-shell. Since 3BFs are not included in the current study, our primary focus will be the impact of different chiral NN and YN interactions as well as their SRG evolution on the separation energies. For the NN interaction we consider the next-to-next-to-next-to-leading order potential from the Idaho group with a regulator of $\varLambda _{N} = 500$ MeV (Idaho-N³LO(500)) [51], and the high-order semilocal momentum-space (SMS) potential regularized with $\varLambda _{N} = 450$ MeV

(SMS N⁴LO+(450)) [53]. Two chiral potentials at next-to-leading order, namely NLO13 and NLO19 [25, 40] with the range of regulators $\varLambda _{Y}$ = 550–650 MeV, are chosen for the YN interaction. In all calculations, contributions of the NN and YN potentials in partial waves higher than $J = 6$ are left out. The high partial waves affect the energies only by a few keV. For simplicity, the electromagnetic part of the NN interaction [54] as well as the Coulomb point-like contribution in some YN channels are not included in the SRG evolution, but only added afterwards. We observed that evolving these interactions changes hypernuclear binding energies only by few keV.

5.1 Extrapolation of the binding energies

Due to the finite truncation in the single-particle Hilbert space, results from the NCSM calculations are dependent on the HO frequency $\omega$ as well as on the model space size ${\mathcal {N}}$. Both parameters can be understood in terms of an ultraviolet and infrared cutoff. Based on this insight, theoretically founded extrapolations can be performed with respect to the infrared cutoff [55,56,57,58]. This is especially interesting for the calculation of expectation values of long range operators, like radii, because the infrared dependence is pronounced in this case. Since we will be most concerned about the ultraviolet dependence, we follow here a simple, but practical approach.

In order to obtain converged binding energies, and, at the same time, to be able to systematically estimate the numerical uncertainties, we follow a two-step procedure as employed in [38]. The first step is to minimize (eliminate) the HO-$\omega$ dependence. For each model space size ${\mathcal {N}}$, we first calculate the binding energies, $E(\omega , {\mathcal {N}})$, for a range of HO-$\omega$ and then utilize the following ansatz,

$$\begin{aligned} E(\omega , {\mathcal {N}}) = E_{{\mathcal {N}}} + \kappa (\text {log}(\omega ) - \text {log}(\omega _{\text {opt}}))^2, \end{aligned}$$

(22)

to extract the lowest binding energy $E_{{\mathcal {N}}}$ for the considered model space ${\mathcal {N}}$ and the corresponding optimal HO frequency $\omega _{\text {opt}}$. As an example, we show in Fig. 3 the HO-$\omega$ dependence of $E({^{4}_{\varLambda } \text {He},0^{+}})$ for model space ${\mathcal {N}}$ varying from 10 to 22. We notice that the optimal frequency $\omega _{\text {opt}}$ shifts to lower values as the model space size ${\mathcal {N}}$ increases, and the $\omega$-dependence of $E(\omega , {\mathcal {N}})$ flattens out as we move forward to the largest model space ${\mathcal {N}}_\text {max}$.

Fig. 3

$E{(^4_{\varLambda }\text {He},0^+)}$ as a function of HO $\omega$. 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). The calculations are based on the Idaho-N³LO(500) (NN) and NLO19(600) (YN) interactions, SRG-evolved to $\lambda _{NN} =1.6$ fm$^{-1}$ and $\lambda _{YN} = 2.00$ fm$^{-1}$, respectively

In the second step, the binding energies with the minimal $\omega$-dependence, $E_{{\mathcal {N}}}$, are used for extrapolating to a converged result in infinite model space assuming an exponential ansatz

$$\begin{aligned} E_{{\mathcal {N}}} = E_{\infty } + A e^{-B{\mathcal {N}}}~. \end{aligned}$$

(23)

The confidence interval for each $E_{{\mathcal {N}}}$ in Eq. (23) can be determined either from the spread of the energy in the vicinity of $\omega _{\text {opt}}$ or from the slope between two successive energies, $E_{{\mathcal {N}}}$ and $E_{{\mathcal {N}}+2}$. The latter is mostly employed in our calculations. It should however be stressed that the two ways of assigning confidence intervals are practically equivalent and lead to the same results within the numerical uncertainties. The determined intervals will serve as a weight for each $E_{{\mathcal {N}}}$ in the model-space fit using the ansatz in Eq. (23). The model-space extrapolation for $E(^4_{\varLambda }\text {He},0^+)$ is illustrated in Fig. 4. The final uncertainty (shaded area) is then taken as the difference between the extrapolated $E_{\infty }$ and $E_{ {{\mathcal {N}}} \text {max}}$.

Fig. 4

${\mathcal {N}}$-dependence of $E{(^4_{\varLambda } \text {He},0^+)}$. The symbols and uncertainties represent results extracted from Eq. (22). The black line is obtained using Eq. (23). The (red) straight line with shaded area indicates the converged result and its uncertainty. Same description of interactions as in Fig. 3

Fig. 5

${\mathcal {N}}$-dependence of $B_\varLambda {(^4_{\varLambda } \text {He},0^+)}$. Same description as in Fig. 4

Fig. 6

${\mathcal {N}}$-dependence of: a $B_{\varLambda }(^5_{\varLambda } \text {He})$, b $B_{\varLambda }(^7_{\varLambda } \text {Li}, \frac{1}{2}^+ 0)$. Same description as in Fig. 4. The Idaho-N³LO NN and NLO19(600) YN potentials are SRG evolved to $\lambda _{NN} =1.6$ fm$^{-1}$ and $\lambda _{YN} =2.6$ fm$^{-1}$, respectively

In hypernuclear physics, we are generally more interested in the so-called $\varLambda -$separation energy, $B_{\varLambda }$, which is defined as the difference between the binding energies of a hypernucleus and of the corresponding core nucleus. Hence,

$$\begin{aligned} B_{\varLambda }(^4_{\varLambda }\text {He}) = E(^3\text {He}) - E(^4_{\varLambda }\text {He})~. \end{aligned}$$

(24)

Following the definition Eq. (24), in principle, one can subtract the separation energy for each $\omega$ and ${\mathcal {N}}$,

$$\begin{aligned} B_{\varLambda }(^4_{\varLambda }\text {He},\omega , {\mathcal {N}} ) = E(^3\text {He},\omega ,{\mathcal {N}}) - E(^4_{\varLambda }\text {He}, \omega , {\mathcal {N}})~, \end{aligned}$$

(25)

and then employ the described two-step procedure to extrapolate the converged $B_{\varLambda }$. We have however observed that for each model space size ${\mathcal {N}}$, the useful ranges of $\omega$ and hence the optimal frequencies $\omega _{\text {opt}}$ for the nuclear core $^3{\text {He}}$ and hypernucleus $^4_{\varLambda }\hbox {He}$ are somewhat different. It is therefore advisable to eliminate the $\omega$-dependence of the binding energies of $^3\text {He}$ and $^4_{\varLambda }$He separately. After that, one subtracts $B_{\varLambda }({\mathcal {N}})$ for every model space ${\mathcal {N}}$

$$\begin{aligned} B_{\varLambda }(^4_{\varLambda }\text {He},{\mathcal {N}}) = E(^3\text {He},{\mathcal {N}}) - E(^4_{\varLambda }\text {He},{\mathcal {N}})~, \end{aligned}$$

(26)

and utilizes the ansatz Eq. (23) to extract the converged result for $B_{\varLambda }(^4_{\varLambda }\text {He})$ together with its uncertainty, see Fig. 5. Clearly, the $\varLambda$-separation energy $B_{\varLambda }(^4_{\varLambda }\text {He})$ exhibits a slightly faster convergence pattern as that of the binding energy $E(^4_{\varLambda }\text {He})$. This tendency is also observed for all other investigated hypernuclei. For completeness, the model-space extrapolations of $B_{\varLambda }(^5_{\varLambda }\text {He})$ and $B_{\varLambda }(^7_{\varLambda }\text {Li},1/2^+)$ are shown in Fig. 6.

It is stressed that there is no fundamental reason that the separation energies monotonically converge with increasing model space, but we observed this monotonic behavior in all systems computed so far. This motivated the use of Eq. (23) for the extrapolation of the separation energies. Note that the resulting energies are consistent with results of FY calculations and/or with a fit of the ${{{\mathcal {N}}}}$ dependence to a constant. The latter way of fitting is less preferable since it generally leads to larger uncertainties.

Let us finally emphasize that, although the described procedure is computationally rather expensive, it allows for a systematic and, most importantly, reliable extraction of the final results of the NCSM calculations. Within the Jacobi-basis formalism such a robust extrapolation is feasible and yields plausible results for light p-shell hypernuclei as one will see in the following sections.

5.2 Benchmark results for $^4_{\varLambda }$He

As mentioned above, to validate the J-NCSM we benchmark our converged results with the binding energies obtained when solving the Faddeev-Yakubovsky equations [23]. More details are given in Appendix B.

Table 1 Ground- and excited-state energies (in MeV) of $^4_\varLambda \hbox {He}$ obtained from the Faddeev-Yakubovsky (FY) and J-NCSM approaches. The calculations are based on the Idaho-N³LO(500) NN interaction, SRG-evolved to $\lambda _{NN} =1.6$ fm$^{-1}$, and the NLO19(600) YN potential, evolved to three different SRG flow values, namely $\lambda _{YN} = 1.6$, 3.0 and 14.0 fm$^{-1}$

The binding energies for the ground state $(0^+)$ and first excited state $(1^+)$ of $^4_{\varLambda }\hbox {He}$ are tabulated in Table 1. Clearly, within the numerical accuracy of better than 20 keV, the two approaches, J-NCSM and FY, agree very nicely.

5.3 Effects of NN chiral interactions on $B_{\varLambda }$

It is known that the nuclear binding energy $E(^3\text {He})$ and consequently $E(^4_{\varLambda }\text {He})$ are very sensitive to the employed NN potentials when three-nucleon (3N) and higher-body forces are not included. This is noticeable in the binding energies of the $^4_{\varLambda }\text {He}(0^+)$ state shown in Fig. 7, obtained for various NN forces: the Idaho-N³LO(500), the improved chiral N²LO and N⁴LO with a configuration-space regulator of $R = 0.9$ fm [59, 60] and the SMS N⁴LO+(450).

Fig. 7

$E(^4_{\varLambda } \text {He},0^+)$ as a function of $\lambda _{YN}$. The calculations are based on the NLO19(600) YN potential and four chiral NN interactions: the Idaho-N³LO(500) (red circles), two improved chiral-N⁴LO (blue triangles) and chiral-N²LO (green diamonds) interactions regularized in configuration space with a cutoff $R=0.9$ fm [59, 60], and the SMS N⁴LO+(450) potential (black crosses). All NN potentials are evolved to a flow parameter of $\lambda _{NN}=1.6$ fm$^{-1}$. The error bars show the estimated numerical uncertainties

Fig. 8

$B_{\varLambda }(^4_{\varLambda } \text {He},0^+)$ as a function of $\lambda _{YN}$. Same description of the curves as in Fig. 7

All NN forces are evolved to an SRG parameter of $\lambda _{NN} =1.6$ fm$^{-1}$. For that value overall the binding energies of the A = 3–6 nuclei are reasonably well described. Of course, this requirement can be fulfilled within a certain range of $\lambda _{NN}$ values so that the actual choice is to some extent arbitrary. The YN potential is evolved to a wide range of flow parameters, $1.0 \le \lambda _{YN}\le 3.0$ fm$^{-1}$. One clearly sees that the binding-energy variations due to different chiral NN forces can be as large as 270 keV. However, being evolved to the same $\lambda _{NN}=1.6$ fm$^{-1}$, these NN potentials have a rather similar impact on the $\varLambda$ removal energy, in particular for low SRG-YN flow parameters $\lambda _{YN} \le 1.6$ fm$^{-1}$ where there is practically no difference in$B_{\varLambda }(^4_{\varLambda }\hbox {He}$, $0^+$), see also Fig. 8. For higher values of $\lambda _{YN}$, the discrepancies among the computed values of $B_{\varLambda }(^4_{\varLambda }\hbox {He}$, $0^+$) somewhat increase but remain relatively small, about 50 keV at most (at $\lambda _{YN} =2.0$ fm$^{-1}$). We stress that a similar behavior is also observed for the $\varLambda$-separation energies of $^4_{\varLambda }\hbox {He} (1^+)$, $^5_{\varLambda }\text {He}$ and $^7_{\varLambda }\text {Li}(\frac{1}{2}^+, 0)$.

Fig. 9

${\varLambda }$-separation energies versus binding energies of the nuclear core: a $^4_{\varLambda } \text {He}(0^+)$ and $^3\hbox {He}$, b $^4_{\varLambda } \text {He}(1^+)$ and $^3\hbox {He}$, c $^5_{\varLambda } \text {He}$ and $^4\hbox {He}$, d $^7_{\varLambda } \text {Li}(\frac{1}{2}^+, 0)$ and $^6\hbox {Li}$. The calculations are based on the Idaho-N³LO(500) (red circles) and the SMS N⁴LO+(450) (blue asterisks) NN potentials, evolved to several values of $\lambda _{NN}$, in combination with the NLO19(600) YN interaction, SRG evolved to $\lambda _{YN} = 2.0$ fm$^{-1}$. The error bars show the estimated numerical uncertainties

Hence, in order to further explore the effect of the NN interaction on $B_{\varLambda }$, we perform calculations using the two most accurate NN potentials, namely Idaho-N³LO(500) and SMS N⁴LO+(450) evolved to several $\lambda _{NN}$ flow variables. It is remarked that, although these two NN potentials describe the available NN scattering data almost perfectly, they indeed have very different matrix elements, particularly in the high-momentum region. It is therefore of great interest to study their predictions for $B_{\varLambda }(A=4-7)$ more carefully. To speed up the convergence of the results, the NLO19(600) YN potential is evolved to a flow parameter of $\lambda _{YN} = 2.0$ fm$^{-1}$. This specific choice of $\lambda _{YN}$ is based on the above observation (cf. Fig. 8) that the largest discrepancy in $B_{\varLambda }$ is generally observed at that flow parameter. The results for the A = 4–7 hypernuclei are displayed in Fig. 9 where the $\varLambda$-separation energies are plotted against the binding energies of the corresponding core nucleus. For the chosen YN flow parameter, the hypernuclei are strongly overbound compared to experiment. We show the results here to emphasize the effect of different NN interactions. For a direct comparison with the experimental separation energies, see below. The energies obtained with the Idaho-N³LO(500) and SMS N⁴LO+(450) potentials are denoted by red squares and blue crosses, respectively. Also, the error bars are added in order to indicate the estimated numerical uncertainties, which in many cases are hardly visible. The light colored bands indicate the variation of the separation energies depending on the binding energy of the core nucleus. Evidently, there is a general trend that stronger nuclear binding energies lead to larger $\varLambda$-separation energies. Furthermore, the overall variations in the $\varLambda$-separation energies of the two states $^4_{\varLambda }\text {He}(0^+, 1^+)$ due to the change in the $^3\hbox {He}$ core binding energies are noticeable, i.e. around 400 keV (see panels (a), (b)). However, the width of the band is rather small, of the order of 80 keV only. For the $^5_{\varLambda }\hbox {He}$ system, panel (c), the variation of $B_{\varLambda }$ stemming from the SRG evolution of the individual NN interactions is roughly 600 keV while the overall discrepancy caused by these two NN potentials can be twice as large. It can be also clearly seen that the width of the band for $^5_{\varLambda }\hbox {He}$ is rather large, about 220 keV. However, given the considerable variation in $B_{\varLambda }(^5_{\varLambda }\text {He})$, the relative width (roughly $22\%$ of the 1 MeV total variation for all NN interactions employed) is of the same order of magnitude as that for the two states of $^4_{\varLambda }\text {He}$. Similarly, the effect of the SRG-NN evolution on $B_{\varLambda }(^7_{\varLambda }\text {Li}, \frac{1}{2}^+)$ for the SRG-YN flow parameter of $\lambda _{YN}=2.0$ fm$^{-1}$ is also pronounced. Here, one of the individual NN potentials, i.e. the Idaho-$\hbox {N}^3\hbox {LO}$(500), already causes a discrepancy in $B_{\varLambda }(^7_{\varLambda }\text {Li}, \frac{1}{2}^+)$ of about 0.8 MeV, which is almost twice the variation in $B_{\varLambda }(^5_{\varLambda }\text {He})$. The total variation when considering both interactions is however similar for both hypernuclei, namely 1.1 MeV. But the relative variation (i.e. the relative width of the colored band in panel (d)) is rather large, about 400 keV ($40\%$ of the 1.1 MeV). For larger $\lambda _{NN}$ ($\lambda _{NN} > 1.6$ fm$^{-1}$), the numerical uncertainties become visible for $^7_\varLambda \hbox {Li}$ and its core. Since the larger $\lambda _{NN}$ significantly increase the width of the band, its width might be further reduced when more converged calculations become available also for these flow parameters. In any case, one can expect from the correlations shown in Fig. 9 that the dependence of $B_{\varLambda }$ on the nuclear interactions can be substantially reduced once the 3N forces are properly included so that nuclear core binding energies are in fair agreement with experiment. Work in this direction is in progress.

5.4 Effects of the NLO YN interactions on $B_{\varLambda }$

We are now in the position to study the impact of the NLO13 and NLO19 YN interactions on the $\varLambda$-separation energies. The two NLO potentials are practically equivalent in terms of describing two-body YN observables. Furthermore, by construction, they reproduce the experimental binding energy of $^3_{\varLambda }\hbox {H}$ within its uncertainty (of order of 50 keV). However, as discussed in Ref. [25], the NLO19 interaction is characterized by a different (somewhat weaker) $\varLambda$-$\varSigma$ transition strength, particularly in the $^3{S}_1$ partial-wave channel, a feature that is believed to be closely related to the strength of chiral YNN forces [25, 35]. The latter is expected to manifest itself in the predictions of observables (e.g. separation energies) for $A \ge 4$ hypernuclei and in infinite nuclear matter. Indeed, it has been found that the NLO19 potential is more attractive in the medium than NLO13 [25]. In addition, in that work, the possible impact of the NLO13 and NLO19 potentials on the $A=4$ hypernucleus has been thoroughly investigated, using the Faddeev-Yakubovsky approach. We provide here again results for the spin-doublet states of $^4_{\varLambda }\text {He}$ for benchmarking. Furthermore, we extend the study to the A = 5–7 hypernuclei. For our purpose, it is sufficient to choose the SMS N⁴LO+(450) potential with $\lambda _{NN} =1.6$ fm$^{-1}$.

Fig. 10

$\varLambda$-separation energies of a $^4_{\varLambda }\text {He}(0^+)$, b $^4_{\varLambda }\text {He}(1^+)$, c $^5_{\varLambda }\text {He}(\frac{1}{2}^{+})$, d $^7_{\varLambda }\text {Li}(1/2^+)$, e $^7_{\varLambda }\text {Li}(3/2^+)$ as a function of the SRG-YN flow parameter $\lambda _{YN}$. Black lines with grey bands represent experimental value of $B_{\varLambda }$ and the uncertainties, respectively. The calculations are based on the NN interaction SMS N⁴LO+(450) with the SRG-NN evolution parameter of $\lambda _{NN}=1.6$ fm$^{-1}$ in combination with the NLO13 (red solid lines) and NLO19 (dashed blue lines) YN potentials for four regulators, $\varLambda _{Y} =$ 500 (triangles), 550 (stars), 600 (crosses) and 650 (circles) MeV

The separation energies $B_{\varLambda }$ of the ground- and first-excited states of the A = 4–7 hypernuclei evaluated for the two NLO YN potentials with various regulators $\varLambda _{Y}$ = 500–650 MeV are presented in Fig. 10. In that calculation, both YN interactions are evolved to the same range of the SRG-YN flow parameters, $0.8 \le \lambda _{YN} \le 3.0$ fm$^{-1}$. For the two states of $^7_{\varLambda }\hbox {Li}$, the calculations have only been performed up to $\lambda _{YN} \le 1.6$ fm$^{-1}$ in order to save some computational resources. Overall, the dependence of $B_{\varLambda }$ on the chiral regulator $\varLambda _{Y}$ is somewhat stronger for the NLO19 than for the NLO13 potential. This, however, does not relate to any physical reason but simply reflects the fact that, in the NLO19 realization, one has less freedom to absorb regulator artifacts into the parameters of the chiral interactions (low-energy constants, LECs) because some of the LECs are determined (and taken over) from fits to NN phase shifts in line with SU(3) flavor symmetry, see [25]. There are also noticeable differences between the $\varLambda$-separation energies obtained with the two interactions, which apparently exceed the $\varLambda _{Y}$-dependence. For all states except $^4_{\varLambda }\text {He}(0^+)$, see panels (b–e), one observes a general tendency toward larger $B_{\varLambda }$ values predicted by NLO19 than those calculated with NLO13. In other words, the interaction with a weaker $\varLambda$–$\varSigma$ conversion potential generally leads to larger $\varLambda$-separation energies. That trend is, however, not clear for the ground state of $^4_{\varLambda }\hbox {He}$ as can be seen in panel a. We remark that a similar (chiral) regulator dependence and sensitivity to the YN potential has been observed in the Faddeev–Yakubovsky results for $A=3$, 4 hypernuclei, computed directly with the bare YN interactions [25]. There, it was already found that NLO19 leads to somewhat stronger binding which might be a result of the weaker $\varLambda$–$\varSigma$ conversion of NLO19 compared to NLO13. The pronounced variations of $B_{\varLambda }$ predicted by the two interactions are a striking evidence for possible contributions of 3BFs to the $\varLambda$ separation energy. These discrepancies are expected to be largely removed once proper chiral YNN forces are taken into account explicitly [61].

Let us mention that the strong sensitivity of the $\varLambda$-separation energies of $^4_{\varLambda }\text {He}(1^+)$ and $^5_{\varLambda }\text {He}$ to the $\varLambda$-$\varSigma$ transition potential can be understood using a simple approximation for the effective spin-dependent $\varLambda N$ potential in s-shell hypernuclei, which can be written as [62, 63]

$$\begin{aligned} \begin{aligned} ^3_{\varLambda }\text {H}: \quad&\tilde{V}_{\varLambda N} \approx \frac{3}{4} V^{s}_{\varLambda N} + \frac{1}{4} V^{t}_{\varLambda N}\\ ^4_{\varLambda }\text {He}(0^+): \quad&\tilde{V}_{\varLambda N} \approx \frac{1}{2} V^{s}_{\varLambda N} + \frac{1}{2} V^{t}_{\varLambda N}\\ ^4_{\varLambda }\text {He}(1^+): \quad&\tilde{V}_{\varLambda N} \approx \frac{1}{6} V^{s}_{\varLambda N} + \frac{5}{6} V^{t}_{\varLambda N}\\ ^5_{\varLambda }\text {He}: \quad&\tilde{V}_{\varLambda N} \approx \frac{1}{4} V^{s}_{\varLambda N} + \frac{3}{4} V^{t}_{\varLambda N},\\ \end{aligned} \end{aligned}$$

(27)

where $V^{s}_{\varLambda N}$ and $V^{t}_{\varLambda N}$ are the singlet- and triplet two-body potentials, respectively. It follows clearly from Eq. (27) that the two states, $^4_{\varLambda }\text {He}(1^+)$ and $^5_{\varLambda }\text {He}$, are dominated by the spin-triplet $\varLambda N$ interaction, which is, as already mentioned, strongly influenced by the $\varLambda$-$\varSigma$ conversion. Interestingly, as can be seen in Fig. 10, the results for $^4_\varLambda \hbox {He} \left( 1^+\right)$, $^5_\varLambda \hbox {He} \left( {1}/{2}^+\right)$ and $^7_\varLambda \hbox {Li} \left( {3}/{2}^+\right)$ in panels (b), (c) and (e) are clearly different for the NLO13 and NLO19 set of interactions. To a lesser extend this can also be seen for $^7_\varLambda \hbox {Li} \left( {1}/{2}^+\right)$ in panel (d). Since $^4_{\varLambda }\text {He}(1^+)$ and $^5_{\varLambda }\text {He}$ are dominated by the $^3 \hbox {S}_1$ interaction,cf. Eq. (27), this suggests that the $^3 \hbox {S}_1$ contribution is also very important for $^7_\varLambda \hbox {Li}$, especially for the ${3}/{2}^+$ state. A future more detailed study will be necessary to validate this hypothesis.

Fig. 11

Probabilities of finding the $\varSigma$ hyperon in the wave functions of a $^4_{\varLambda }\text {He}(0^+)$, b $^4_{\varLambda }\text {He}(1^+)$, c $^5_{\varLambda }\text {He}({1}/{2}^+)$, d $^7_{\varLambda }\text {Li}(1/2^+)$, e $^7_{\varLambda }\text {Li}(3/2^+)$ as a function of SRG-YN flow parameter $\lambda _{YN}$. Same NN potential, symbols and lines as in Fig. 10

In this context, the probabilities of finding a $\varSigma$ particle in the hypernuclear wave functions ($P_{\varSigma }$) are of great interest, too. Clearly, they are an indication for the strength of the $\varLambda$-$\varSigma$ conversion of the YN interaction. Moreover, it can be expected that there are some correlations to the charge-symmetry breaking (CSB) of $\varLambda$ separation energies of mirror hypernuclei as well [23, 24]. Our calculated $\varSigma$-probabilities for A = 4–7 hypernuclei obtained with the two NLO potentials are shown in Fig. 11. It is interesting that in all systems $P_\varSigma$ decreases with decreasing $\lambda _{YN}$ for $\lambda _{YN} \ge 1$ $\hbox {fm}^{-1}$ but increases again for $\lambda _{YN} < 1$ $\hbox {fm}^{-1}$. Additionally, The results displayed in panel (a) clearly indicate a noticeable dependence of $P_{\varSigma }(^4_{\varLambda }\text {He}, 0^+)$ on the chiral cutoff $\varLambda _{Y}$. That regulator dependence, however, becomes somewhat less visible for all other states, see panels (b-e). Also, the variation of the $\varSigma$-probabilities caused by the two chiral interactions is most pronounced for $^4_{\varLambda }\text {He}(0^+)$. This is exactly opposite to the observations for the $\varLambda$-separation energies as discussed above. Moreover, there is an overall tendency toward larger $P_{\varSigma }$ predicted by the interaction with a stronger $\varLambda$-$\varSigma$ transition (i.e. NLO13) although it is somewhat blurred by the regulator dependence. We further note that, while there is a visible difference between the $\varSigma$-probabilities of the s-shell spin-doublet states (in particular for the predictions of NLO13), the p-shell doublet $P_{\varSigma }(^7_{\varLambda }\text {Li}, {1}/{2}^+)$ and $P_{\varSigma }(^7_{\varLambda }\text {Li}, {3}/{2}^+)$ are quite similar for both interactions. Clearly, one sees that the $\varLambda$-separation energies and $\varSigma$-probabilities in A = 4–7 hypernuclei are somewhat correlated. However, we do not observe a definite one-to-one correlation between the two quantities.

5.5 Correlation of $\varLambda$-separation energies

Fig. 12

Correlations of $\varLambda$-separation energies between $^5_{\varLambda }\hbox {He}$ and a $^3_{\varLambda }\hbox {H}$, b the $0^+$ state of $^4_{\varLambda }\hbox {He}$ (red) and $^4_{\varLambda }\hbox {H}$ (blue), c the $1^+$ state of $^4_{\varLambda }\hbox {He}$ (red) and $^4_{\varLambda }\hbox {H}$ (blue), d $^6_{\varLambda }\hbox {He}$ (red) and $^6_{\varLambda }\hbox {Li}$ (blue), e $^7_{\varLambda }\text {Li}({1}/{2}^{+},0)$ and f $^7_{\varLambda }\text {Li}({3}/{2}^{+},0)$, for a wide range of flow parameters $\lambda _{YN}$. The error bars represent numerical uncertainties which are small in most of the cases. The experimental $\varLambda$-separation energy for $^5_{\varLambda }\text {He}$ is from [52]. The results for other systems are taken from a [52], b–c [64] for $^4_{\varLambda }\text {He}$ (black asterisk) and $^4_{\varLambda }\text {H}$ (grey square), d [65] for $^6_{\varLambda }\text {He}$ (black asterisk) and $^6_{\varLambda }\text {Li}$ (grey square), e [52] and f [66]. The Idaho-$\hbox {N}^3$LO(500) evolved to 1.6 $\hbox {fm}^{-1}$ and NLO19(600) was used for the NN and YN interaction, respectively

In Sect. 5.4, we have observed surprisingly similar trends of the $\varLambda$-separation energies for all investigated hypernuclei with respect to the running SRG-YN flow parameter $\lambda _{YN}$. This probably hints at some intriguing correlations between the $\varLambda$-separation energies of these systems. In order to quantitatively study these correlations, we compute $B_{\varLambda }$ for all considered hypernuclei, for the same range of $\lambda _{YN}$ evolution parameters, and compare the results with each other for selected values of $\lambda _{YN}$. It is known that $^5_{\varLambda }\hbox {He}$ is the experimentally best studied hypernucleus so far. Also, our J-NCSM results for this hypernucleus are well-converged. We therefore use $^5_{\varLambda }\hbox {He}$ as a benchmark system and plot $B_{\varLambda }(^5_{\varLambda }\text {He})$ against the separation energies of other hypernuclear systems (A = 3–7), see Fig. 12. For that, we choose Idaho-N³LO(500) evolved to an SRG-NN flow variable of $\lambda _{NN}=1.6$ fm$^{-1}$ for the NN interaction and NLO19 with a regulator of $\varLambda _{Y} = 600$ MeV for the YN interaction. However, we want to emphasize that similar trends are observed for SMS N⁴LO+(450) and in combination with other YN interactions, see also [48]. Let us first look at the correlation between the $\varLambda$ removal energies of the $^5_{\varLambda }\hbox {He}$ hypernucleus and of the hypertriton. Here $B_{\varLambda }(^3_{\varLambda }\text {H})$ are computed within the Faddeev–Yakubovky approach since NCSM calculations are very difficult for this weakly bound system. The correlation plot is presented in panel (a) of Fig. 12. Here each symbol represents the numerical $B_{\varLambda }$ of the two systems calculated at the same flow parameter $\lambda _{YN}$, and it also includes the estimated uncertainties that are small in most of the cases. The straight line is obtained from a linear fit to the results, reminding one of the Tjon line between the binding energies of $^4\hbox {He}$ and $^3\hbox {He}$ [67,68,69,70,71,72]. We observe a nearly perfect linear correlation between $B_{\varLambda }(^3_{\varLambda }\text {H})$ and $B_{\varLambda }(^5_{\varLambda }\text {He})$ for flow parameters up to $\lambda _{YN}=2.0$ fm$^{-1}$ and a slight deviation from the straight line as $\lambda _{YN}$ further increases. The latter can be attributed to the possible contribution of 3BFs [48]. Interestingly, the correlation line goes through the experimental $\varLambda$-separation energies of the two systems at $\lambda _{YN}=0.836$ fm$^{-1}$. The value of $\lambda _{YN}$, at which the $^5_{\varLambda }\text {He}$ hypernucleus is properly described, will be referred to as the magic flow parameter $\lambda ^{m}_{YN}$. For that value, the separation energy of $^3_\varLambda \hbox {H}$ is 92 keV. Using the bare NLO19(600) and the same NN interaction, we found 119 keV which is in reasonable agreement with the result at $\lambda ^{m}_{YN}$. Obviously, the concrete value of $\lambda ^{m}_{YN}$ will depend on the YN interactions as well as their regulators.

The correlation plots for the ground and excited states of $^4_{\varLambda }\text {He}/^4_{\varLambda }\text {H}$ are displayed in panels (b) and (c), respectively. While there is a strictly linear correlation between the separation energies $B_{\varLambda }(^4_{\varLambda }\text {He} / ^4_{\varLambda }\text {H}, 1^+)$ and $B_{\varLambda }(^5_{\varLambda }\text {He})$, the correlation line for $B_{\varLambda }(^4_{\varLambda }\text {He} / ^4_{\varLambda }\text {H}, 0^+)$ and $B_{\varLambda }(^5_{\varLambda }\text {He})$ exhibits a small loop to the right for large values of $\lambda _{YN}$, $\lambda _{YN} \ge 2.4$ fm$^{-1}$ similar to the behavior of the correlation line for $B_{\varLambda }(^3_{\varLambda }\text {H})$ and $B_{\varLambda }(^5_{\varLambda }\text {He})$. Also, from panels (b) and (c), one easily notices almost identical results for the isospin mirrors $^4_{\varLambda }\hbox {He}$ and $^4_{\varLambda }\hbox {H}$. This is because there are no CSB terms in the employed version of the chiral YN potential. The CSB effect arising from the point Coulomb interactions is included in the calculation, but its contribution is minor [73, 74]. It is interesting that, at the magic flow parameter, $\lambda ^{m}_{YN}=0.836$ fm$^{-1}$, the experimental value of $B_{\varLambda }(^4_{\varLambda }\text {He}, 1^+)$ is exactly reproduced while the ground state is somewhat underbound. Furthermore, at this $\lambda ^{m}_{YN}$ our J-NCSM results for the spin doublet of $^4_{\varLambda }\text {He}$, $B_{\varLambda }(0^+ (1^+)) = 1.57 (0.97)$ MeV, are surprisingly close to the those obtained within the exact Faddeev-Yakubovsky method using the non-evolved bare YN interactions, $B_{\varLambda }(0^+ (1^+)) = 1.61 (1.18)$ MeV. The slight deviation between the two results is consistent with the size of 3BFs expected from the power counting of chiral EFT [25].

Similarly, almost perfectly linear correlations are also found between $B_{\varLambda }(^5_{\varLambda }\text {He})$ and the ground-state energies $E_{\varLambda }$ of the p-shell $^6_{\varLambda }\text {He}$ and $^6_{\varLambda }\text {Li}$ hypernuclei, panel (d), as well as the $\varLambda$-separation energies $B_{\varLambda }$ of the ground and first excited states in $^7_{\varLambda }\text {Li}$, panels (e) and (f), respectively. Note that the resonance energies $E_{\varLambda }(^6_{\varLambda }\text {Li}/^6_{\varLambda }\text {He})$ are computed as the difference between the hypernuclear binding energies $E(^6_{\varLambda }\text {Li}/^6_{\varLambda }\text {He})$ and the binding energy $E(^4\text {He})$. This removes most of the NN-interaction dependence. In panel (d), one notices a pronounced difference in the binding energies $E_{\varLambda }$ of $^6_{\varLambda }\text {He}$ and $^6_{\varLambda }\text {Li}$ (about 1.08 MeV), which simply results from different contributions of the Coulomb interactions of the two nuclear cores $^5{\text {He}}$ and $^5{\text {Li}}$. We remark that the NLO19(600) YN potential with the magic flow parameter $\lambda ^{m}_{YN}=0.836$ fm$^{-1}$ underbinds the $^6_{\varLambda }\text {He}/^6_{\varLambda }\text {Li}$ systems while it slightly overbinds the first excited state in $^7_{\varLambda }\hbox {Li}$. The obtained $\varLambda$-separation energy for the ground state, $B_{\varLambda }(^7_{\varLambda }\text {Li},\,{1}/{2}^+) =5.59 \pm 0.01$ MeV, is, however, in very good agreement with the result from emulsion experiments, $B_{\varLambda }(^7_{\varLambda }\text {Li}, {1}/{2}^+) = 5.58 \pm 0.03$ MeV [52]. It should be noted that counter experiments reported a somewhat larger value for $^7_{\varLambda }\text {Li}({1}/{2}^{+},0)$, namely $B_{\varLambda }(^7_{\varLambda }\text {Li}, {1}/{2}^+) = 5.85 \pm 0.13 \pm 0.1$ MeV [75].

The observed linear correlations between the separation energies of different hypernucler systems is rather striking and interesting. It will be important to examine those correlations using different YN bare interactions in order to check whether this useful property is a universal feature or just a signature of the chiral interactions. Nevertheless, our finding for the chiral forces with SRG evolution suggests that the missing SRG-induced three-body forces might be parameterized by only one adjustable parameter (effects of SRG-induced higher-body forces on $B_{\varLambda }$ are expected to be insignificant [35]). If this is the case, one is able to minimize the effects of the omitted three-body forces by tuning the SRG-YN flow parameters $\lambda _{YN}$ to the magic value for which a particular hypernucleus, for example $^5_{\varLambda }\text {He}$, is properly described. This magic flow parameter $\lambda ^{m}_{YN}$ then can serve as a good starting point for hypernuclear calculations requiring a SRG-YN evolution – which, in turn, may provide a good opportunity to study hypernuclear structure as well as the YN forces in a less expensive but realistic approach. A possible application of this finding has been considered in [15, 48]. In this context let us mention that similar linear correlations have been also observed in Ref. [19] for the double-$\varLambda$ hypernuclei $^{\ \, 5}_{\varLambda \varLambda }\hbox {H}$ and $^{\ \, 6}_{\varLambda \varLambda }\hbox {He}$.

As discussed in Ref. [25], the contribution of chiral 3BFs is comparable to the uncertainty at NLO of approximately 200–300 keV for $A=4$. The full $\lambda _{YN}$ dependence of the result is an order of magnitude larger than what is expected for 3BFs by chiral power counting. This situation is very different from that for ordinary nuclei where SRG-induced and chiral 3BFs are of comparable size. Wirth and Roth have pointed out that the size of the SRG-induced 3BFs is probably enhanced because the $\varSigma$ contribution is significantly weakened when $\lambda _{YN}$ is lowered [35]. Our observation here is that, for extreme values of $\lambda _{YN}$ below 1 $\hbox {fm}^{-1}$, the $P_\varSigma$ value increases again and the overbinding disappears. For such $\lambda _{YN}$, the contribution of 3BFs is again in line with the expectation from the chiral power counting. Especially, it seems to be neglible for $^3_\varLambda \hbox {H}$.

6 Conclusions

In this work, we have extended the nuclear J-NCSM to describe baryonic systems with strangeness $S=-1$. The inclusion of the strangeness degree of freedom significantly complicates the implementation of the approach in part because the particle conversion $\varLambda$-$\varSigma$ is explicitly taken into account. Accordingly, the Jacobi basis now consists of two orthogonal subsets, characterized by the $\varLambda$ and $\varSigma$ hyperons. For the applications of the two-body NN and YN forces, we introduced two auxiliary bases that explicitly single out the involved NN and YN pairs, respectively. Like the coefficients of fractional parentage, the expansion coefficients can also be computed in a preparatory step separately from any binding-energy calculations. Once they are known, the evaluation of the many-body Hamiltonian matrix elements in the Jacobi basis (and therefore the energy calculations) are straightforward.

As a first application of the Jacobi NCSM, we utilized the approach to investigate hypernuclear systems with A = 4–7. Here, the $\varLambda$-separation (binding) energies are extracted systematically via a two-step procedure that enables an effective removal of the HO-$\omega$ sensitivity of the final results as well as a reliable estimation of the numerical uncertainties.

We performed the energy calculations based on various SRG-evolved chiral interactions. In particular, we considered the Idaho N³LO and SMS N⁴LO NN potentials in combination with the next-to-leading order YN interactions, NLO13 and NLO19. We found that at low values of the SRG YN flow parameter, $\lambda _{YN} \le 1.4$ fm$^{-1}$, the separation energies are not very sensitive to the NN potentials. The dependence somewhat increases for higher $\lambda _{YN}$, however, the relative variations remain quite similar for all systems.

It turned out that, for some of the considered hypernuclei, there are large differences between the predictions of the two practically phase-equivalent YN potentials NLO13 and NLO19. Those can be attributed to possible (but so far neglected) contributions of chiral three-body (YNN) forces [61]. We also observed that there are almost perfect linear correlations between the $\varLambda$ separation energies of the A = 4–7 hypernuclei calculated for a wide range of the SRG-YN flow parameter. Interestingly, at the magic value $\lambda _{YN}^{m}$ that yields the empirical $B_{\varLambda }(^5_{\varLambda }\text {He})$, the separation energies of $^3_{\varLambda }\text {H}$ and $^4_{\varLambda }\text {He}(0^+,1^+)$ are in good agreement with the results for the non-evolved YN interactions (at least within the expected contributions of the chiral 3BF), while the one for $^7_{\varLambda }\text {Li}$ is surprisingly close to the experiment. This may suggest that by tuning the SRG parameter such that the $^5_{\varLambda }\hbox {He}$ hypernucleus is correctly reproduced, one can effectively minimize the effects of the missing SRG-induced 3BF. Therefore, the special flow parameter $\lambda _{YN}^{m}$ can be a good starting point for hypernuclear calculations that require an SRG evolution. Such calculations will be useful to develop improved YN interactions. Eventually, taking SRG-induced and chiral 3BF into account will be necessary. Work in this direction is in progress.

Appendices

Appendix A: Transition $\big \langle \big (\alpha ^{*(1)}\big )^{*(Y)}| \alpha ^{*(YN)} \big \rangle$

The states $| \alpha ^{*(YN)} \rangle$ and $\big |\big (\alpha ^{*(1)} \big )^{*(Y)}\big \rangle$ with directions of momenta are illustrated in Figs. 13 and 14, respectively. The explicit quantum numbers of these states are

(A.1)

with

(A.2)

and

(A.3)

Thereby, the total HO quantum number is given by $\tilde{{{\mathcal {N}}}} = {\mathcal {N}}^{*(1)}_{(A-1)} +2 n_Y +l_Y$ and ${{{\mathcal {N}}}} = {{{\mathcal {N}}}}_{YN} + {{{\mathcal {N}}}}_{A-2} + 2n_\lambda +\lambda$. The transition $\big \langle \big (\alpha ^{*(1)}\big )^{*(Y)}| \alpha ^{*(YN)} \big \rangle$ can be interpreted as a transformation between different Jacobi coordinates. We can therefore make use of the general Jacobi-coordinate transformation formula Eq. (11) in [38]. For that, we first need to specify the directions of the relative motions of particles (subclusters) in the two states $| \alpha ^{*(YN)} \rangle$ and $\big |\big (\alpha ^{*(1)} \big )^{*(Y)}\big \rangle$. These directions are depicted in Figs. 13 and 14.

Fig. 13

$|\alpha ^{*(YN)}\rangle$ state with directions of momenta

Fig. 14

$\big |\big (\alpha ^{*(1)}\big )^{*(Y)}\big \rangle$ state with directions of momenta

Comparing the definitions of our two states $| \alpha ^{*(YN)} \rangle$ and $\big |\big (\alpha ^{*(1)} \big )^{*(Y)}\big \rangle$ with the corresponding ones in Eq. (11) in [38], one notices that the directions of the relative momenta are the same, however, the ordering of the coupling of the angular momenta and isospins in $\big |\big (\alpha ^{*(1)} \big )^{*(Y)}\big \rangle$ and $|\alpha \rangle _{(13)2}$ are different. The recoupling from $(\tilde{T}_{A-2}\, t_N)T^{*(1)}_{A-1} \text{ to } ( t_N\, \tilde{T}_{A-2})T^{*(1)}_{A-1}$ requires a simple phase factor,

$$\begin{aligned} \begin{aligned} \big | \big (\tilde{T}_{A-2}\, t_N)T^{*(1)}_{A-1} \big \rangle = (-1)^{\tilde{T}_{A-2} + t_N - T^{*(1)}_{A-1}} \big |( t_N\, \tilde{T}_{A-2})T^{*(1)}_{A-1}\big \rangle . \end{aligned} \end{aligned}$$

(A.4)

And, changing the coupling $\big |(\tilde{J}_{A-2}(l_N s_N)I_N)J^{*(1)}_{A-1} \big \rangle$ to $\big |(l_N(s_N \tilde{J}_{A-2})S_{A-1})J^{*(1)}_{A-1} \big \rangle$ can be done with the help of 6j-symbols

$$\begin{aligned}&|\big (\tilde{J}_{A-2}( l_N s_N)I_N) J^{*(1)}_{A-1}\big \rangle = (-1)^{I_N + 2\tilde{J}_{A-2} + l_N + s_N}\nonumber \\&\quad \times \sum _{S_{A-1}=\tilde{J}_{A-2} + s_N} {\hat{I}}_N {\hat{S}}_{A-1} \left\{ \begin{array}{ccc} \tilde{J}_{A-2} &{} s_N &{} S_{A-1}\\ l_N &{} J^{*(1)}_{A-1} &{} I_N \end{array} \right\} \nonumber \\&\quad \times \big |(l_N(s_N \tilde{J}_{A-2})S_{A-1})J^{*(1)}_{A-1}\big \rangle , \end{aligned}$$

(A.5)

where the abbreviation ${\hat{I}}_{N} = \sqrt{2I_N +1}$, etc., is introduced. Now taking into account Eqs. (A.5) and (A.4) and then making use of the Jacobi-coordinate transformation formula in [38], one obtains

$$\begin{aligned} \begin{aligned}&\big \langle \big ( \alpha ^{*( 1)}\big )^{*(Y)} | \alpha ^{*(YN)} \big \rangle \\&\quad = \delta _{{\mathcal {N}} \tilde{{\mathcal {N}}}} \, \delta _{t_{Y} \tilde{t}_Y} \, \delta _{\tilde{T}_{A-2} T_{A-2}} \delta _{\tilde{J}_{A-2} J_{A-2}} \delta _{\mathcal {\tilde{N}}_{A-2} {\mathcal {N}}_{A-2}}\delta _{{\tilde{\zeta }}_{A-2} \zeta _{A-2}}\\&\qquad \times {\hat{I}}_{N}\,{\hat{I}}_Y \, {\hat{J}}_{YN} \, {\hat{S}}_{YN}\,{\hat{I}}_{A-2} {\hat{J}}^{*(1)}_{A-1} \, {\hat{T}}^{*(1)}_{A-1}\, {\hat{T}}_{YN}\\&\qquad \times (-1)^{3J_{A-2} + 2T_{A-2} + T_{YN} + S_{YN} +\lambda + t_Y + l_Y + t_N + l_N + I_N+ 1}\\&\qquad \times \sum _{S_{A-1}=\tilde{J}_{A-2} + s_N} (-1)^{S_{A-1}} {\hat{S}}^{2}_{A-1} \left\{ \begin{array}{ccc} {J}_{A-2} &{} s_N &{} S_{A-1}\\ l_N &{} J^{*(1)}_{A-1} &{} I_N \end{array} \right\} \\&\qquad \times \sum _{L,S} {\hat{L}}^2 {\hat{S}}^2 \left\{ \begin{array}{ccc} l_{N} &{} S_{A-1} &{} J^{*(1)}_{A-1}\\ l_{Y} &{} s_{Y} &{} I_{Y}\\ L &{} S &{}{J} \end{array}\right\} \left\{ \begin{array}{ccc} l_{YN} &{} S_{YN} &{} J_{YN}\\ \lambda &{} J_{A-2} &{} I_{A-2}\\ L &{} S &{} {J} \end{array}\right\} \\&\qquad \times \langle n_N \,l_{N}\, n_Y\, l_{Y} : L \,|\, n_{YN} \,l_{YN} \, n_{\lambda } \, \lambda : L \rangle _{d}\\&\qquad \times \left\{ \begin{array}{ccc} s_Y &{} s_N &{} S_{YN}\\ J_{A-2} &{} S &{} S_{A-1} \end{array} \right\} \left\{ \begin{array}{ccc} t_Y &{} t_N &{} T_{YN}\\ T_{A-2} &{} {T} &{} T^{*(1)}_{A-1} \end{array} \right\} , \end{aligned} \end{aligned}$$

(A.6)

where the HO bracket $\langle n_N \,l_{N}\, n_Y\, l_{Y} : L \,|\, n_{YN} \,l_{YN} \, n_{\lambda } \, \lambda : L \rangle _{d}$ follows the same convention as in [76] with the mass ratio given by

$$\begin{aligned} d=\frac{(A-2) \, m(t_Y)}{(A-1)\,m_N + m(t_Y)}. \end{aligned}$$

(A.7)

Appendix B: Faddeev-Yakubovsky equations for hypernuclei

Faddeev-Yakubovsky equations in momentum space are a well established tool to solve the Schrödinger equations for light hypernuclei with $A=3$ or 4 [21, 23]. We use $A=4$ results to benchmark the NCSM and to provide results for bare interactions. For $A=3$, momentum space is much more efficient for the representation of wave functions and, therefore, for the solution of the Schrödinger equation than HO wave functions. The weak binding of $^3_\varLambda \hbox {H}$ leads to an extremely slow convergence of the energy with respect to ${{\mathcal {N}}}$.

Our solution follows Ref. [77]. For $A=3$, we need to solve a set of coupled Faddeev equations

$$\begin{aligned} | \psi _A \rangle= & {} G_0 t_{12} (1-P_{12}) | \psi _B\rangle \nonumber \\ | \psi _B \rangle= & {} G_0 t_{31} \left( | \psi _A\rangle -P_{12} | \psi _B \rangle \right) \end{aligned}$$

(B.8)

for two Faddeev amplitudes $| \psi _A \rangle$ and $| \psi _B \rangle$. Here, we assume that particles 1 and 2 are nucleons and particle 3 is the hyperon. The permutation operator $P_{12}$ exchanges all coordinates and quantum numbers of the two nucleons. $G_0$ is the free three-baryon propagator. The two off-shell t-matrices $t_{12}$ and $t_{31}$ are solutions of the Lippmann-Schwinger equation of subsystem (12) or (31), respectively. For the solution, we use two momentum Jacobi bases:

$$\begin{aligned} | p_{12} p_3 \alpha \rangle = \left| p_{12} p_3 \!\left[ (l_{12}s_{12})j_{12}\! \left( l_3 \frac{1}{2} \right) \! I_3 \right] \! J (t_{12} t_Y) T M_T \right\rangle \nonumber \\ \end{aligned}$$

(B.9)

and

$$\begin{aligned} | p_{31} p_2 \beta \rangle = \left| p_{31} p_2 \!\left[ (l_{31}s_{31})j_{31}\! \left( l_2 \frac{1}{2} \right) \!I_2 \right] \! J (t_{31} \frac{1}{2}) T M_T \right\rangle .\nonumber \\ \end{aligned}$$

(B.10)

Here, $p_{ij}$ denotes the magnitude of the relative momentum in subsystem (ij) and $p_k$ the relative momentum of particle k relative to the other two particles. The angular dependence is expanded in corresponding orbital angular $l_{ij}$ and $l_k$. These are coupled to the spin of the two-baryon subsystem $s_{ij}$ to the total angular momentum of the subsystem $j_{ij}$. $l_k$ and the spin 1/2 of the third baryon couple to the spectator angular momentum $I_k$. Finally, the total angular momentum J of the three-body system is obtained by coupling $j_{ij}$ and $I_k$. The isospin of the pair $t_{ij}$ is coupled either with the isospin of the hyperon $t_Y=0,1$ or with the isospin of the nucleon 1/2 to the total isospin T and its third component $M_T$. For the hypertriton $M_T=0$ and $T=0$ is the by far dominant component of the wave function. For the solution, the two Faddeev amplitudes $| \psi _A \rangle$ and $| \psi _B \rangle$ are expanded in terms of their natural set of basis states $| p_{12} p_3 \alpha \rangle$ and $| p_{31} p_2 \beta \rangle$, respectivley. Because of the short-ranged hypernuclear interactions, the t-matrices converge quickly with respect to partial waves which is then also true for the Faddeev components if expressed in their natural set of basis states. Note that transitions between these states are required in order to solve the equations. In this work, we restrict the partial wave so that $j_{ij}\le 6$. This ensures that energies are converged to better than 1 keV. We also calculated the wave function

$$\begin{aligned} | \varPsi \rangle = | \psi _A \rangle + (1-P_{12})|\psi _B\rangle \end{aligned}$$

(B.11)

and checked explicitly that the expectation value of the operator agrees with the energy obtained by solving the Faddeev equations. For this check, we need to include the $T=1$ and $T=2$ contributions to reach an accuracy of the order of 1 keV.

For $A=4$ hypernuclei, we need to solve a set of five coupled Yakubovsky equations

$$\begin{aligned} | \psi _{1A} \rangle= & {} G_0 t_{12} (P_{13}P_{23}+P_{12}P_{23}) \nonumber \\&\left[ | \psi _{1A}\rangle + | \psi _{1B}\rangle + | \psi _{2A}\rangle \right] \nonumber \\ | \psi _{1B} \rangle= & {} G_0 t_{12} \left[ (1-P_{12}) (1-P_{23}) | \psi _{1C}\rangle \right. \nonumber \\&\left. + (P_{13}P_{23}+P_{12}P_{23}) | \psi _{2B}\rangle \right] \nonumber \\ | \psi _{1C} \rangle= & {} G_0 t_{14} \left[ | \psi _{1A}\rangle + | \psi _{1B}\rangle + | \psi _{2A}\rangle -P_{12} | \psi _{1C}\rangle \right. \nonumber \\&\left. +P_{13}P_{23} | \psi _{1C}\rangle +P_{12}P_{23} | \psi _{2B}\rangle \right] \nonumber \\ | \psi _{2A} \rangle= & {} G_0 t_{12} \left[ (P_{12}-1)P_{13} | \psi _{1C}\rangle + | \psi _{2B}\rangle \right] \nonumber \\ | \psi _{2B} \rangle= & {} G_0 t_{34} \left[ | \psi _{1A}\rangle + | \psi _{1B}\rangle + | \psi _{2A}\rangle \right] \end{aligned}$$

(B.12)

for the five Yakobovsky components $| \psi _{1A}\rangle$, $| \psi _{1B}\rangle$, $| \psi _{1C}\rangle$, $| \psi _{2A}\rangle$ and $| \psi _{2B}\rangle$. Each of these components is expanded in terms of its natural Jacobi coordinate, respectively, as defined below

$$\begin{aligned}&| p_{12} p_3 q_4 \alpha _A \rangle \nonumber \\&\quad = \left| p_{12} p_3 q_4 \left[ \left[ (l_{12}s_{12})j_{12} \, \left( l_3 \frac{1}{2} \right) I_3 \right] j_{123} \ \left( l_4 \frac{1}{2} \right) I_4 \right] J \quad \right. \nonumber \\&\qquad \left. \left[ (t_{12} \frac{1}{2}) \tau _{123} t_Y \right] T M_T \right\rangle \ \nonumber \\&| p_{12} p_4 q_3 \alpha _B \rangle \nonumber \\&\quad = \left| p_{12} p_4 q_3 \left[ \left[ (l_{12}s_{12})j_{12} \, \left( l_4 \frac{1}{2} \right) I_4 \right] j_{124} \ \left( l_3 \frac{1}{2} \right) I_3 \right] J \quad \right. \nonumber \\&\qquad \left. \left[ (t_{12} t_Y ) \tau _{124} \frac{1}{2} \right] T M_T \right\rangle \ \nonumber \\&| p_{14} p_2 q_3 \alpha _C \rangle \nonumber \\&\quad = \left| p_{14} p_2 q_3 \left[ \left[ (l_{14}s_{14})j_{14} \, \left( l_2 \frac{1}{2} \right) I_2 \right] j_{124} \ \left( l_3 \frac{1}{2} \right) I_3 \right] J \quad \right. \nonumber \\&\qquad \left. \left[ (t_{14} \frac{1}{2}) \tau _{124} \frac{1}{2} \right] T M_T \right\rangle \ \nonumber \\&| p_{12} p_{34} q \beta _A \rangle \nonumber \\&\quad = \left| p_{12} p_{34} q \left[ \left[ (l_{12}s_{12})j_{12} \, \lambda \right] I \ (l_{34}s_{34})j_{34} \right] J\right. \nonumber \\&\qquad \left. (t_{12} t_{34}) T M_T \right\rangle \nonumber \\&| p_{34} p_{12} q \beta _B \rangle \nonumber \\&\quad \!= \! \left| p_{34} p_{12} q \left[ \left[ (l_{34}s_{34})j_{34} \, \lambda \right] I \ (l_{12}s_{12})j_{12} \right] J \right. \nonumber \\&\qquad \left. (t_{34} t_{12}) T M_T \right\rangle \ .\!\!\!\! \end{aligned}$$

(B.13)

Here, the free propagator $G_0$ and the t-matrices $t_{ij}$ are of course embedded into the four-baryon system. The coupling scheme is much more complicated than in the three-baryon system. Now there are two types of Jacobi coordinates required. The first three basis sets are of the “3+1” type. Here, three momenta $p_{ij}$, $p_k$ and $q_l$ are required that are relative momenta within the pair ij, of particle k with respect to pair ij and of particle l with respect to the three-body subsytem ijk. Additionally to the quantum numbers of the three-body system, we have now introduced $j_{ijk}$ and $\tau _{ijk}$ for the total angular momentum and isospin of the three-body subsystem. J, T and $M_T$ are the total angular momentum, isospin and third component of isospin of the four-baryon system. We have again omitted the spins and isospins of the two baryons in the inner most subsystem since only $t_4=t_Y$ differs from 1/2. The last two basis sets are of the “2+2” type. Here, relative momenta of two two-body subsytems $p_{ij}$ and $p_{kl}$ are introduced together with angular momenta and isospins for these subsystems. Additionally, the relative momentum of the two pairs q and its angular momentum $\lambda$ is required. In order to finally define the total four-body angular momentum, an additional intermediate angular momentum I needs to be introduced as seen in the definition of the states.

For four-baryon states, it is not sufficient to constrain the two-body angular momenta in order to get a finite number of partial waves. Additional constraints on other angular momenta are necessary. For the calculations of this work, we chose $j_{ij}\le 5$, $l_i\le 6$, $\lambda \le 6$, $l_{ij}+l_k+l_l \le 10$ and $l_{ij}+l_{kl}+\lambda \le 10$. In order to save computational resources, we restrict ourselves to the by far most important isospin component $T={1}/{2}$, although the contribution to the energy of the Yakubovsky equations induces an uncertainty of 10 keV. Interestingly, the other isospin components are more important when calculating the expectation value for which they contribute approximately 20 keV. Note that also our J-NCSM results are based on the dominant isospin components only. We therefore need to take an uncertainty of approximately 20 keV in the four-baryon systems into account due to missing isospin components.

Once the Yakubovsky components are found, we obtain the wave function by

$$\begin{aligned} | \varPsi \rangle= & {} (1+P_{13}P_{23}+P_{12}P_{23}) | \psi _{1A} \rangle \nonumber \\&+ (1+P_{13}P_{23}+P_{12}P_{23}) | \psi _{1B} \rangle \rangle \nonumber \\&+ (1-P_{12}) (1+ P_{13}P_{23}+P_{12}P_{23}) | \psi _{1C} \rangle \nonumber \\&+(1+P_{13}P_{23}+P_{12}P_{23}) | \psi _{2A} \rangle \rangle \nonumber \\&+(1+P_{13}P_{23}+P_{12}P_{23}) | \psi _{2B} \rangle \ . \end{aligned}$$

(B.14)

This briefly summarizes the Faddeev-Yakubovsky approach as we used it for benchmarking our J-NCSM results.