1 Introduction

The prospect of an ample production of baryons with charm offered by facilities such as the LHC at CERN [1,2,3], RHIC at BNL [4], J-PARC and KEK in Japan [5, 6], or FAIR in Germany [7,8,9] has led to a renewed interest in the in-medium properties of such baryons [10,11,12,13] and also in the question whether they, and notably the lightest charmed baryon, the \({\varLambda }_c\)(2286), could form bound states with ordinary matter [14,15,16,17,18,19]. In fact, there is a long history of speculations about possible bound systems involving the \({\varLambda }_c\) [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36] that started soon after the discovery of charmed baryons [37, 38] (see also the recent reviews [39, 40]). In principle, charmed nuclei could be produced by means of charm exchange or associate charm production reactions [41, 42], in analogy to the ones widely used in hypernuclear physics. However, the experimental production of charmed nuclei is difficult due to the short lifetimes of D-meson beams which makes it necessary to place the target as close as possible to the D-meson production point, and due to the kinematics of the reactions: the charmed particles are formed with large momentum making their capture by a target-nucleus improbable. Because of these difficulties, up to now, only three albeit controversial candidates have been reported by an emulsion experiment carried out in Dubna in the mid-1970s [43,44,45,46,47].

The \({\varLambda }_cN\) forces employed in the past investigations were predominantly derived within the meson-exchange framework, see Refs. [14, 17] for recent examples, often utilizing SU(4) flavor symmetry in one form or the other. Lately, also the constituent quark model [48] or a combination of meson-exchange and quark model [16] have been considered. Independently of that, in general, the resulting potentials turned out to be fairly attractive. Interestingly, a rather different picture emerged from recent lattice QCD (LQCD) simulations by the HAL QCD collaboration [49, 50]. Those studies, based on unphysical quark masses corresponding to pion masses of \(m_\pi \) = 410–700 MeV, suggest that the \({\varLambda }_cN\) and \({\varSigma }_cN\) interactions could be significantly less attractive than what had been proposed in the phenomenological studies mentioned above. In Ref. [51], an extrapolation of the HAL QCD results to the physical point was presented, using chiral effective field theory (EFT) [52, 53] as guideline. It revealed that the \({\varLambda }_cN\) interaction at the physical point is expected to be somewhat stronger than for large pion masses, however, still only moderately attractive and, specifically, considerably less attractive than most of the phenomenological predictions for the \({\varLambda }_cN\) interaction.

In the present work, we use the \({\varLambda }_cN\) interaction from Ref. [51] to explore the binding energies of charmed nuclei. In the literature different conventions for naming \({\varLambda }_c\) nuclei have been used in the past. We adopt here the standard nomenclature for nuclei with the proper generalization [30], spelled out for hypernuclei in Sect. I.B of Ref. [54]. It takes into account that the characterizing letter(s) for the nucleus depends on its total charge and not just on the number of protons. For example, when adding the (uncharged) \({\varLambda }\) to hydrogen (\(^2\hbox {H}\)) one gets the hypertriton (\(^3_{{\varLambda }}\hbox {H}\)) but adding the positively charged \({\varLambda }_c\) leads to \(^{\ 3}_{{\varLambda }_c}\hbox {He}\). Similarly, based on this convention, the bound state of \(\varLambda _c\) and \(^{208}\hbox {Pb}\) is \(^{209}_{{\varLambda }_c}\hbox {Bi}\).

The lightest nuclei considered, the \(A=3\) and \(A=4\) systems \(^{\ 3}_{{\varLambda }_c}\hbox {He}\) and \(^{\ 4}_{{\varLambda }_c}\hbox {He}\), are investigated by solving corresponding Faddeev–Yakubovsky equations. Indeed, for the three-body system, bound states have been reported in Ref. [16] (with total angular momentum \(J=1/2\) and 3/2) and Ref. [15] (for \(J=3/2\)). Note, however, that some of the \({\varLambda }_cN\) interactions employed in Ref. [16] are so strongly attractive that they even predict two-body bound states (in the \(^1S_0\) as well as in the \(^3S_1\) partial wave) with binding energies comparable to that of the deuteron. For calculating heavier charmed nuclei, namely from \(^{\ 5}_{\varLambda _c}\hbox {Li}\) onward to \(^{209}_{\varLambda _c}\hbox {Bi}\), a perturbative many-body approach is utilized that allows one to obtain the \(\varLambda _c\) single-particle bound states in the different nuclei from the corresponding \(\varLambda _c\) self-energy. Results for charmed nuclei computed within this framework have been reported recently in Ref. [17], for a set of \(Y_c N\) interactions deduced from an early version [55] of the hyperon–nucleon (YN) meson-exchange potential of the Jülich Group [56] via SU(4) symmetry arguments. For these interactions, the \(^{\ 5}_{\varLambda _c}\hbox {Li}\) system (\(^{\ 5}_{\varLambda _c}\hbox {He}\) in the nomenclature used in Ref. [17]) turned out to be already bound. In this context, let us mention that the HAL QCD collaboration has likewise reported results for charmed nuclei [49]. The calculations were performed with the \({\varLambda }_cN\) potentials extracted from the lattice simulations at pion masses 410–700 MeV, but using the physical masses of the \({\varLambda }_c\) and the considered nuclei. Binding energies for \({\varLambda }_c\) bound to \(^{12}\hbox {C}\), \(^{28}\hbox {Si}\), \(^{40}\hbox {Ca}\), \(^{58}\hbox {Na}\), \(^{90}\hbox {Zr}\), and \(^{208}\hbox {Pb}\) were reported.

Since after the publication of the \({\varLambda }_cN\) interaction [51] new results from the HAL QCD collaboration became available that include now the \({\varSigma }_cN\) channel [50], we also revisit the \({\varLambda }_cN\) interaction in order to explore in how far the inclusion of a direct interaction in the \({\varSigma }_cN\) channel modifies the extrapolation of the \({\varLambda }_cN\) interaction from the results/masses of the lattice simulations to the physical point, and in how far it influences the predictions for charmed nuclei. It turns out that adding/considering the interaction in the \({\varSigma }_cN\) channel has very little influence on the \({\varLambda }_cN\) scattering results and also not on those for \({\varLambda }_c\) nuclei. May be this is not too surprising in view of the fact that the thresholds of the two channels are separated by almost 170 MeV. In any case, for completeness, we report predictions for the \({\varSigma }_cN\)S-wave phase shifts based on our extrapolation of the HAL QCD results [50].

Another extrapolation of the lattice results for \({\varSigma }_cN\), performed in heavy baryon chiral perturbation theory and taking into account heavy quark spin symmetry, has been performed recently [57]. Let us mention already now that neither in that work nor in our calculation any signal for resonances or bound states in the \({\varSigma }_cN\) channel are found. The existence of such resonances has been suggested in Ref. [58], where the \({\varLambda }_cN\)\({\varSigma }_cN\)\({\varSigma }_c^*N\) interaction was investigated within the framework of the quark delocalization color screening model. Bound states (and resonances) in the \({\varSigma }_cN\)\(^3S_1\) partial wave around the \({\varSigma }_cN\) and \({\varSigma }_c^*N\) thresholds were also predicted in Ref. [59] based on a \(Y_c N\) interaction from meson exchange supplemented by short-range repulsion from a quark exchange model.

The paper is structured in the following way: a summary of the main characteristics of the \(\varLambda _c N\) and \(\varSigma _c N\) potentials is presented in Sect. 2. Results of the properties of \(\varLambda _c\) in infinite nuclear matter, and light and heavier nuclei are reported in Sect. 3. Finally, a brief summary and some concluding remarks are given in Sect. 4. The appendix summarizes results for \({\varSigma }_cN\) scattering.

2 The \({\varLambda }_cN\) and \({\varSigma }_cN\) potentials

The \({\varLambda }_cN\) and \({\varSigma }_cN\) interactions are constructed by using chiral EFT as guideline, following the scheme employed in our studies of the \(\varLambda N\) and \(\varSigma N\) systems [60,61,62,63]. We summarize the essentials below. More details of the approach can be found in Ref. [51]. The \(Y_c N\) potential consists of contact terms and contributions from pion exchange. The former are given by

$$\begin{aligned}&V(^1S_0) = {\tilde{C}}_{^1S_0} + \tilde{D}_{^1S_0} m^2_\pi + ({C}_{^1S_0} + {D}_{^1S_0} m^2_\pi )\, ({p}^2+{p}'^2), \nonumber \\&V(^3S_1) = {\tilde{C}}_{^3S_1} + \tilde{D}_{^3S_1} m^2_\pi + ({C}_{^3S_1} + {D}_{^3S_1} m^2_\pi )\, ({p}^2+{p}'^2), \nonumber \\&V(^3D_1 -\, ^3S_1) = {C}_{\varepsilon _1}\, {p'}^2, \nonumber \\&V(^3S_1 -\, ^3D_1) = {C}_{\varepsilon _1}\, {p}^2, \end{aligned}$$
(1)

for the partial waves considered in the present study. Here \(p = |\mathbf{p}\,|\) and \({p}' = |\mathbf{p}\,'|\) are the initial and final center-of-mass (c.m.) momenta in the \({\varLambda }_cN\) or \({\varSigma }_cN\) systems. The quantities \({\tilde{C}}_{i}\), \({\tilde{D}}_{i}\), \({C}_{i}\), \({D}_{i}\) are low-energy constants (LECs) that are fixed by a fit to lattice data (phase shifts) by the HAL QCD collaboration at \(m_\pi =410\) MeV and 570 MeV. The \(m_\pi \) dependence in Eq. (1) is motivated by the corresponding expression in the standard Weinberg counting up to next-to-leading order (NLO) [52, 53] but differs from it by the term proportional to \(m^2_\pi \, ({p}^2+{p}'^2)\) which is nominally of higher order. Nevertheless, we included that term in Ref. [51] because it allowed us to obtain an optimal description of the LQCD results at \(m_\pi =410\) MeV as well as 570 MeV. We consider this as a prerequisite for a well constrained extrapolation to lower pion masses. Without such a term the phase shifts by HAL QCD for \(m_\pi =410\) MeV would be underestimated at low energies, as exemplified by results shown in Figs. 6 and 7 of Ref. [57].

The contribution of pion exchange to the \({Y}_cN\) potential is given by

$$\begin{aligned} V^{OPE}_{{Y_c N\rightarrow Y'_c N}} =-f_{{Y_c Y'_c\pi }}f_{{NN\pi }} \frac{\left( \mathbf{\sigma }_1\cdot \mathbf{q} \right) \left( \mathbf{\sigma }_2\cdot \mathbf{q} \right) }{ \mathbf{q}^{\,2}+m_{\pi }^2}, \end{aligned}$$
(2)

where \(\mathbf{q}\) is the transferred momentum, \(\mathbf{q} = \mathbf{p'} - \mathbf{p}\). The coupling constants \(f_{BB'\pi }\) are related to the axial-vector strength via \(f_{BB'\pi } = g_A^{BB'}/2F_\pi \) with \(F_\pi \) being the pion decay constant (\(F_\pi \approx 93\) MeV). The coupling constant for the \({\varLambda }_c{\varSigma }_c\pi \) vertex can be determined from the experimentally known \({\varSigma }_c\rightarrow {\varLambda }_c\pi \) decay rate, see Refs. [64, 65]. For the \({\varSigma }_c{\varSigma }_c\pi \) coupling constant lattice QCD results [66] are employed. To be concrete, \(g_A^{{\varSigma }_c{\varSigma }_c}=0.71\) [66] and \(g_A^{{\varLambda }_c{\varSigma }_c}=0.74\) [64, 65] are used, together with \(g_A^{NN}=1.27\) [67]. Note that \(f_{{\varLambda }_c{\varLambda }_c\pi } \equiv 0\) under the assumption that isospin is conserved.

Besides the coupling constants at the physical point, one needs also their \(m_\pi \) dependence:

$$\begin{aligned} f_{BB'\pi } (m^2_\pi ) = \frac{g_A^{BB'}(m^2_\pi )}{2\,F_\pi (m^2_\pi ) }. \end{aligned}$$
(3)

Results for the dependence of \(F_\pi \) on \(m^2_\pi \) are available from lattice simulations [68]. Based on that work the values \(F_\pi \approx 112\) MeV at \(m_\pi = 410\) MeV and \(F_\pi \approx 129\) MeV at \(m_\pi = 570\) MeV were deduced and employed in Ref. [51]. With regard to the dependence of \(g^{BB'}_A\) on \(m^2_\pi \) lattice QCD simulations indicate a rather small variation, at least for \(g^{{\varSigma }_c{\varSigma }_c}_A\) and \(g^{NN}_A\) where concrete results are available [66, 69]. Because of that the dependence of the \(g_A\)’s on \(m_\pi \) was neglected in Ref. [51] and the values at the physical point were used throughout.

There was no information from LQCD on the \({\varSigma }_cN\) interaction at the time when the study in Ref. [51] was performed, and, thus, the interaction in the \({\varSigma }_cN\) channel was not considered. Nonetheless, the coupling of \({\varLambda }_cN\) to \({\varSigma }_cN\) via pion exchange was already included. Due to its long-range nature, this coupling plays an important role in case of the \(\varLambda N\) and \(\varSigma N\) systems [61, 62]. Since \(f_{{\varLambda }_c{\varSigma }_c\pi }\) is empirically known and only slightly smaller than \(f_{{\varLambda }{\varSigma }\pi }\) [61] the coupling between \({\varLambda }_cN\) and \({\varSigma }_cN\) via pion exchange should still be of relevance. Indeed, the effective contribution to the \({\varLambda }_cN\) interaction, \(V_{{\varLambda }_cN \rightarrow {\varLambda }_cN} \sim V^{OPE}_{{\varLambda }_cN \rightarrow {\varSigma }_cN} \, G_{{\varSigma }_cN}\, V^{OPE}_{{\varSigma }_cN \rightarrow {\varLambda }_cN}\), is expected to be smaller for energies around the \({\varLambda }_cN\) threshold as compared to the situation for \({\varLambda }N\), but just by a factor 2-3. The reduction is due to the larger mass difference, \(M_{{\varSigma }_c}-M_{{\varLambda }_c} \approx 167\) MeV versus \(M_{\varSigma }-M_{\varLambda } \approx 78\) MeV, that enters in the corresponding Green’s functions \(G_{{\varSigma }_cN}\) or \(G_{{\varSigma }N}\). In any case, formally two-pion exchange contributions to the \({\varLambda }_cN\) interaction involving the \({\varSigma }_c\) do arise at NLO [61], and it can be expected that the piece with a \({\varSigma }_cN\) intermediate state provides the dominant contribution [70]. In the actual calculation, this (reducible) two-pion exchange contribution to the \({\varLambda }_cN\) potential is generated by solving a coupled-channel Lippmann–Schwinger (LS) equation (see below). Further (irreducible) NLO contributions from two-pion exchange [61] have been omitted in [51] for simplicity reasons. It was assumed that those can be effectively absorbed into the contact terms.

For the present study, we add a direct \({\varSigma }_cN\) interaction. The potential is determined in the same way as the one for \({\varLambda }_cN\), by using results of the HAL QCD collaboration for the corresponding \(^3S_1\) phase shift at \(m_\pi =410\) MeV and 570 MeV [50]. The main goal of this extension is to explore in how far the inclusion of a direct \({\varSigma }_cN\) interaction modifies the \({\varLambda }_cN\) results achieved earlier. Of particular interest is the question, whether it influences the results for charmed nuclei that we are concerned with here. Results for \({\varSigma }_cN\) scattering itself are discussed and summarized in the appendix.

The reaction amplitudes are obtained from the solution of a coupled-channel LS equation for the interaction potentials. After partial-wave projection [60], the equation is given by

$$\begin{aligned}&T^{\nu '\nu ,J}_{\rho '\rho }(p',p;\sqrt{s})\quad =V^{\nu '\nu ,J}_{\rho '\rho }(p',p) \nonumber \\&\qquad + \sum _{{\rho ''},\nu ''}\int _0^\infty \frac{dp''p''^2}{(2\pi )^3} \, V^{\nu '\nu '',J}_{\rho '\rho ''}(p',p'') \nonumber \\&\qquad \times \frac{2\mu _{\rho ''}}{p^2_{\rho ''}-p''^2+i\eta } T^{\nu ''\nu ,J}_{\rho ''\rho }(p'',p;\sqrt{s}) \ , \end{aligned}$$
(4)

where the label \(\rho \) indicates the channels (\({\varLambda }_cN\), \({\varSigma }_cN\)) and the label \(\nu \) the partial wave. The quantity \(\mu _\rho \) signifies the pertinent reduced mass. The on-shell momentum in the intermediate state, \(p_{\rho }\), is defined by \(\sqrt{s}=\sqrt{M^2_{B_{1,\rho }}+p_{\rho }^2}+\sqrt{M^2_{B_{2,\rho }}+p_{\rho }^2}\).

Since the integral in the LS equation (4) is divergent for the chiral potentials specified above, a regularization scheme needs to be introduced [71, 72]. For that purpose the potentials in the LS equation are cut off in momentum space by multiplication with a regulator function [60, 61], \(f(p',p) = \exp \left[ -\left( p'^4+p^4\right) /\varLambda ^4\right] \). In Ref. [51], cutoff values \(\varLambda =500\)–600 MeV were employed, in line with the range that yielded the best results in NLO studies of the \({\varLambda }N\) and \({\varSigma }N\) interactions [61, 62] and in comparable investigations of NN scattering [73]. For a more detailed discussion on this topic in the context of the employed \({\varLambda }_cN\) interaction see Ref. [51], where also further references on the issue of regularization can be found. The variations of the results with the cutoff, reflecting uncertainties due to the regularization, will be indicated by bands.

The baryon masses corresponding to the LQCD simulations at \(m_\pi = 410\) and 570 MeV are taken from Ref. [49]. For the calculation at the physical point, we use the masses \(M_{N} = 938.92\) MeV, \(M_{{\varLambda }_c} = 2286.5\) MeV, and \(M_{{\varSigma }_c} = 2455\) MeV.

3 \({\varLambda }_c\) nuclei and matter properties

In this section, we report results on the properties of the \({\varLambda }_c\) in infinite nuclear matter and on charmed nuclei. In the pertinent calculations, we employ \(Y_c N\) interactions extrapolated from results of lattice simulations by the HAL QCD collaboration [49, 50] to the physical point. In particular, we use the \({\varLambda }_cN\) potential from Ref. [51], where there is no direct \({\varSigma }_cN\) interaction, and the two potentials \(Y_c N\)-A and \(Y_c N\)-B, introduced and described in the appendix, which include a direct \({\varSigma }_cN\) interaction.

Fig. 2
figure 1

\({\varLambda }_c\) in symmetric nuclear matter. Top: the \({\varLambda }_c\) s.p. potential \(U_{{\varLambda }_c}(k_{{\varLambda }_c}=0)\) as a function of the Fermi momentum \(k_F\) in comparison to \(U_{\varLambda }(k_{\varLambda }=0)\) [62] (NLO13 with dash-dotted border; NLO19 with dashed border) for cutoff variations \({\varLambda }=500\)–600 MeV. Bottom: momentum dependence of \(U_{{\varLambda }_c}\) (solid band) and \(\varepsilon _{{\varLambda }_c}\) (hatched band) at the Fermi momentum \(k_F=1.35\hbox { fm}^{-1}\)

3.1 \({\varLambda }_c\) in infinite nuclear matter

In order to investigate the properties of the \({\varLambda }_cN\) interaction in nuclear matter, we perform a Brueckner–Hartree–Fock calculation where we adopt the so-called discontinuous prescription when solving the Bethe–Goldstone equation. We follow closely our corresponding calculation for the \({\varLambda }N\) interaction [74]. In that work and similar ones (see, e.g., Refs. [75, 76]), the reader can find details how to solve the Bethe–Goldstone equation and how the single-particle (s.p.) potential \(U_{{\varLambda }_c}\) is determined self-consistently together with the G-matrices for a specific nuclear matter density \(\rho \) (or Fermi momentum \(k_F\)).

In Fig. 1a, we present results for the dependence of \(U_{{\varLambda }_c}(k_{{\varLambda }_c}=0)\) on the Fermi momentum, in comparison to those for the \({\varLambda }\) hyperon obtained with the NLO interaction from Refs. [61, 62]. In Fig. 1b, we display the dependence of \(U_{{\varLambda }_c}(k_{{\varLambda }_c})\) and the s.p. energy, \(\varepsilon _{{\varLambda }_c}(k_{{\varLambda }_c})= k^2_{{\varLambda }_c}/2m_{{\varLambda }_c} + U_{{\varLambda }_c}(k_{{\varLambda }_c})\), on the \({\varLambda }_c\) momentum at the Fermi momentum \(k_F=1.35\hbox { fm}^{-1}\), i.e., at nuclear matter saturation density. The in-medium predictions for \({\varLambda }_c\) are based on the \(Y_c N\)-A potential (cf. appendix). Results for the \({\varLambda }_c\) properties for the alternative fit \(Y_c N\)-B, considered in the appendix, and for the \({\varLambda }_cN\) potential from Ref. [51] practically coincide with the ones for \(Y_c N\)-A and are therefore not shown.

Table 1 \({\varLambda }_cN\) scattering lengths (in fm) and partial-wave contributions to the s.p. potential \( U_{{\varLambda }_c} (k_{{\varLambda }_c} = 0)\) (in MeV) at \(k_F = 1.35 \ \mathrm{fm}^{-1}\). Results are shown for the \(Y_c N\)-A potential which includes a direct \({\varSigma }_cN\) interaction (cf. appendix), and for the \({\varLambda }_cN\) interaction from Ref. [51]. The cutoff values used (\(\varLambda = 500,\,600\) MeV) are indicated in brackets. For comparison corresponding results for the \(\varLambda N\) interactions NLO13 [61] and NLO19 [62] are given

Partial-wave contributions to \( U_{{\varLambda }_c} (k_{{\varLambda }_c} = 0)\) at \(k_F=1.35\hbox { fm}^{-1}\) are listed in Table 1. Note that the contributions of the P waves come solely from two-pion exchange involving the \({\varSigma }_cN\) intermediate state. The total potential depth amounts to around \(-20\div -18\) MeV and is quite insensitive to whether a direct \({\varSigma }_cN\) interaction is included or not. As a reminder, the “empirical” value in case of the \({\varLambda }\) hyperon is \(-30\div -27\) MeV [54]. Comparing the \({\varLambda }_c\) results with the \({\varLambda }\) case in detail, one can see that the contribution in the \(^1S_0\) partial wave is reduced by roughly a factor three. This is well in line with the corresponding interaction strengths; the \({\varLambda }_cN\) scattering length is also about a factor three smaller than that for \({\varLambda }N\), see Table 1. For the \(^3S_1\) partial wave the \({\varLambda }_cN\) and \({\varLambda }N\) contributions (for NLO13 [61]) are of comparable magnitude, despite of the fact that the \({\varLambda }_cN\) interaction is less attractive as reflected in the corresponding scattering lengths which is about 30–50 % smaller than that for \({\varLambda }N\) scattering. Obviously, for \({\varLambda }_cN\) the dispersive effects, which play an important role for the contribution of that partial wave in case of the \({\varLambda }\) [62], are smaller because the \({\varLambda }_cN\)\({\varSigma }_cN\) coupling is weaker due to a weaker transition potential and/or due to the larger threshold separation. Apparently, that reduced effect compensates for the somewhat less attractive \({\varLambda }_cN\) interaction. When comparing with the results for the NLO interaction from 2019 [62], where the \({\varLambda }N\)\({\varSigma }N\) transition potential is noticeably weaker, one sees a clear correlation between the smaller \(^3S_1\) scattering length and the reduced contribution to \(U_{{\varLambda }_c}\) (cf. Table 1).

Finally, let us compare our nuclear matter results with other predictions for \(U_{{\varLambda }_c}\) found in the literature. Reference [31] contains some results for \(U_{{\varLambda }_c}\) based on an \({Y}_cN\) potential that is adapted from one of the YN meson-exchange potentials by the Nijmegen Group by imposing SU(4) flavor symmetry. In that work, a value of \( U_{{\varLambda }_c} (k_{{\varLambda }_c} = 0)\approx -25\) MeV at nuclear matter saturation density has been found. However, note the large contributions from P waves in that study. The two S wave alone yield only around \(-3.5\) MeV. A study utilizing parity-projected QCD sum rules [10] reports a potential depth of \(\sim -20\) MeV for \({\varLambda }_c\) at nuclear matter saturation density. Yasui, in a perturbative approach based on a heavy-quark effective theory, finds a \({\varLambda }_c\) binding energy of around \(-25\div -20\) MeV in nuclear matter [13]. Though to some extent surprising, it is interesting to see that the achieved results are all fairly similar, despite of the different interactions and approaches employed.

3.2 \(^{\ 3}_{{\varLambda }_c}\)He and \(^{\ 4}_{{\varLambda }_c}\)He systems

For \(A=3\) and 4 charmed nuclei, Faddeev–Yakubovsky calculations are performed in the same way as in former studies of hypernuclei [62, 77, 78]. As discussed in Ref. [51], based on the results for \({\varLambda }\)-hypernuclei and the relative strengths of the \({\varLambda }_cN\) and \({\varLambda }N\) interactions, one can guess which light charmed nuclei could be bound. For that the pertinent mixtures of the spin-singlet and spin-triplet \({\varLambda }_cN\) (\({\varLambda }N\)) interaction for s-shell nuclei are relevant [62, 79] and, of course, the reduction of the kinetic energy associated with the \({\varLambda }_c\) as a consequence of its larger mass [30]. In view of the fact that, for the considered \(Y_c N\) interactions, the \({\varLambda }_cN\)\(^1S_0\) scattering length is only one third of the one for \({\varLambda }N\), while there is somewhat less difference in the \(^3S_1\) state (cf. Table 1), binding of light systems is expected to be mainly possible for charmed nuclei with a dominating spin-triplet \({\varLambda }_cN\) contribution, i.e. \(^{\, 3}_{{\varLambda }_c}\hbox {He}\) (\(J=\frac{3}{2}^+\)), \(^{\, 4}_{{\varLambda }_c}\hbox {He}\) (\(1^+\)), and \(^{\, 5}_{{\varLambda }_c}\hbox {Li}\) [62, 77, 79].

Additionally, whereas the Coulomb interaction is of less importance for \(\varLambda \) separation energies of hypernuclei, its contribution to \({\varLambda }_c\) separation energies of charmed nuclei is often decisive for binding [30]. For the solution of the Faddeev–Yakubovsky equations here, we take the Coulomb interactions fully into account as described in Ref. [80].

First, in order to benchmark our few-body calculations, we devised \({\varLambda }_cN\) interactions that mimic the effective range parameters predicted by a \({\varLambda }_cN\) potential obtained in the constituent quark model by Garcilazo et al. [48]. In that model, the triplet interaction is much stronger (\(a_{^3S_1} = -2.31\) fm) than the one in the singlet channel (\(a_{^1S_0} = -0.86\) fm), cf. Table 3 in that reference. Consequently, and in line with the above arguments and the explorations in Ref. [77], a bound state for the \(J=\frac{3}{2}\) state of \(^{\ 3}_{{\varLambda }_c}\)He has been reported, with a \({\varLambda }_c\) separation energy of approximately 140 keV including Coulomb [15]. Since our interactions do not reproduce the phase shifts of the quark model perfectly over a larger energy region, we use two different realizations with cutoff 600–700 MeV. We find separation energies between 60 keV and 264 keV. This includes a variation of 60 keV due to different NN interactions. The uncertainty due to different cutoffs of the \({\varLambda }_cN\) interaction is larger than the one due to different NN interactions. We also found that no bound state exists for that interaction for \(^{\ 3}_{{\varLambda }_c}\)He with \(J=\frac{1}{2}\). This result confirms the earlier calculations of Garcilazo et al. [48].

We then performed calculations for the \({\varLambda }_cN\) potentials from Ref. [51], and the interactions \(Y_c N\)-A and \(Y_c N\)-B of this work. Since all of these potentials predict a considerably weaker interaction in the \(^3\hbox {S}_1\) partial wave, none of the charmed \(A=3\) nuclei are found to be bound. This remains even true when the Coulomb interaction is not taken into account.

Note that model 4 employed by Gibson et al. [30] has \({\varLambda }_cN\)\(^1S_0\) and \(^3S_1\) scattering lengths close to those of the \(Y_c N\) interactions considered by us. No bound state for \(A=3\) was found in that work, but a bound \(^{\, 4}_{{\varLambda }_c}\)He with a \({\varLambda }_c\) separation energy of approximately 1.25 MeV was predicted, after including an estimate for the contribution of the Coulomb interaction.

Table 2 \({\varLambda }_c\) separation energies \(E_{{\varLambda }_c}\) and binding energies with respect to breakup up into four baryons, E, for the \(J=1^+\) state of \(^{\ 4}_{{\varLambda }_c}\)He. The NN interaction at order N\(^4\)LO+ with a cutoff of 450 MeV of Ref. [81] was used leading to \(E(^3\mathrm{H})=-8.141\) MeV. Expectation values for the kinetic energy \(\langle T \rangle \), the NN potential energy \(\langle V_{NN} \rangle \), and the \(Y_c N\) potential energy \(\langle V_{Y_c N} \rangle \) are also given. The probability that the four-baryon state has orbital angular momentum zero and two ( P(S) and P(D) ) is listed together with the probability \(P({\varSigma }_c)\) for the \({\varSigma }_c\) component. Energies are given in MeV and probabilities in %
Table 3 Same as Table 2 for the \(J=0^+\) state of \(^{\ 4}_{{\varLambda }_c}\hbox {He}\)

The \(A=4\) results for the interactions considered in the present work are summarized in Tables 2 and 3. In all of our calculations, we found that the \(J=1^+\) state is more bound than the \(J=0^+\) state. This differs from the situation for strangeness (\(^{\, 4}_{{\varLambda }}\)He) where the \(J=0^+\) state is more strongly bound [62]. The opposite behavior is due to differences in the relative strength of the spin-singlet and triplet interactions. Specifically, since the spin-singlet interaction is much weaker in the \({\varLambda }_cN\) potential inferred from LQCD [51] as compared to that for the \({\varLambda }N\) system [62], nuclei which are dominated by the spin-triplet interaction, like \(^{\ 4}_{{\varLambda }_c}\)He with \(J=1^+\), are more likely to be bound. For a general discussion of the role of the spin dependence for s-shell hypernuclei see Refs. [62, 79].

As obvious from Table 2, the results are very independent on whether the direct \({\varSigma }_cN\) interaction has been included or omitted. The cutoff of the \(Y_c N\) interactions has a much larger effect on the energies than the inclusion of the direct interaction. Even more striking is the observation that the results for \(A=4\) are identical for the potentials \(Y_c N\)-A and \(Y_c N\)-B. This shows unmistakably the insensitivity of the predicted bound-state properties on the \({\varSigma }_cN\) channel. The separation energies for the \(J=1^+\) state are between 100 keV and 370 keV, a clear evidence that \(A=4\) charmed nuclei could be bound. The binding energies are, however, somewhat smaller than predicted in Ref. [30]. We believe that this is partly due to omitting tensor interactions in the \(Y_c N\) in Ref. [30] which is fully taken into account in our calculations. For comparison, we have also used the interactions that simulate the quark-model potential of Ref. [48]. This interaction clearly provides stronger binding leading to 1.2–2.1 MeV separation energy depending on the cutoff used.

A few properties of the resulting wave functions are summarized in Tables 2 and 3, too. First of all, it is interesting to compare the expectation value of the Hamiltonian to the energy. For the numerical calculations, we need to restrict the number of partial waves. The most significant restriction is that the algebraic sum of all orbital angular momenta is less or equal 8. We checked that the solution of the Yakubovsky equations is converged such that the energy E is accurate to approximately 10 keV. The expection values differ at most by 20 keV. The slightly larger differences is due to a slower convergence of the wave functions compared to the Yakubovsky components. The good agreement of both numbers is a confirmation of the consistency of the numerical calculation. We also show separate expectation values of the kinetic energy, the NN potential energy and the \(Y_c N\) potential energy. It sticks out that the NN potential energy is very similar for all considered bound states. Clearly, the nuclear core is not very much distorted by the presence of the charmed hyperon. The expectation value of the \(Y_c N\) interaction is mainly dependent on the cutoff of the interaction and less sensitive to the \(\varSigma _c N\) contribution as can be seen from the similarity of the results of \(\varLambda _c N\) [51], \(Y_c N\)-A and \(Y_c N\)-B.

Finally, we give probabilities for the total orbital angular momentum of 0 and 2 P(S) and P(D). P-waves and F-waves only give a negligible contribution. Obviously, the tensor components of the NN and \(Y_c N\) interactions induce a D-wave contribution of approximately 7%. The probability to find a \(\varSigma _c\) is, similar to the \(\varSigma \) in ordinary hypernuclei, small and depends strongly on the cutoff of the \(Y_c N\) interaction.

The outcome for the \(J=0^+\) state is compiled in Table 3. In this case, we do not find a bound \(^4_{\varLambda _c}\hbox {He}\) for the interactions with a cutoff of 500 MeV. Such a state is however close to being bound as can be seen from the expectation values based on an approximate solution of the Yakubovsky equation, cf. Table 3. Also for the larger cutoff, the \(\varLambda _c\) separation energy is only 100–180 keV. Therefore, for our interactions, we can neither confirm nor exclude the existence of a \(0^+\) bound state. We note that we do find a bound \(0^+\) state for the interactions that simulate the quark model potential of Ref. [48]. In that case the separation energy is between 130 and 470 keV depending on the cutoff used.

Given that the dependence on the NN interaction was smaller than the variation with the employed regulator in the \(Y_c N\) interaction for \(A=3\), we refrain from repeating the computationally very expensive \(A=4\) calculations for different NN interactions. We do not expect that the results will be significantly different for other choices.

3.3 Heavier charmed nuclei

Now we consider the energy of the \(\varLambda _c\) single-particle bound states in heavier nuclei. To such end, we follow a perturbative many-body approach whose starting point is a nuclear matter G-matrix derived from the bare \(Y_c N\) interactions described in Sect. 2 and the appendix. This G-matrix is then used to calculate the self-energy of the \(\varLambda _c\) in the finite nucleus. Solving the Schrödinger equation with this self-energy, finally, we are able to determine the energies of all the single-particle bound states of the \(\varLambda _c\) in the nucleus. This approach also provides the real and imaginary parts of the \(\varLambda _c\) optical potential at positive energies and, therefore, allows one to study the \(\varLambda _c\)-nucleus scattering properties. This method was already used to study the properties of the nucleon [82], the \(\varDelta \) isobar [83] and the \(\varLambda \) and \(\varSigma \) hyperons [84,85,86] in finite nuclei, and very recently also those of the \(\varLambda _c\) using a meson-exchange \(Y_cN\) interaction [17]. A comprehensive description of the method can be found in these works and the interested reader is referred to any of them for details.

Results for \(^{\ 5}_{{\varLambda }_c}\hbox {Li}\), \(^{13}_{{\varLambda }_c}\hbox {N}\) , \(^{17}_{{\varLambda }_c}\hbox {F}\), \(^{41}_{{\varLambda }_c}\hbox {Sc}\), \(^{91}_{{\varLambda }_c}\hbox {Nb}\) and \(^{209}_{{\varLambda }_c}\hbox {Bi}\) are summarized in Table 4 for the \({\varLambda }_cN\) interaction from Ref. [51] and the two potentials \(Y_c N\)-A and \(Y_cN\)-B with inclusion of a direct \({\varSigma }_cN\) interaction. We note that all charmed nuclei considered consist of a closed-shell nuclear core plus a \(\varLambda _c\) sitting in a single-particle state. We note also that, although the \(\varLambda _cN\) interaction with the cutoff 600 MeV is more attractive, the \(\varLambda _c\) single-particle bound states predicted in this case are actually less bound. This is due to dispersive effects [87,88,89] in the calculation of the G-matrix which suppress the contribution from the \(\varLambda _cN\rightarrow \varSigma _cN\rightarrow \varLambda _cN\) coupling. That coupling is significantly larger for the 600 MeV cutoff and, accordingly, likewise the reduction of the overall attraction. Before analyzing the results, we would like to point out that, as discussed in Ref. [90], the approach followed tends to underestimate the energies of the \(\varLambda \) hyperon single-particle bound states for light hypernuclei such as \(^{\,5}_{\varLambda }\hbox {He}\). Accordingly, we expect \(^{\ 5}_{{\varLambda }_c}\hbox {Li}\) to be somewhat more strongly bound than what is suggested by the values given in Table 4.

Table 4 Energy of \(\varLambda _c\) single-particle bound states (in MeV) of several charmed nuclei from \(^{\ 5}_{\varLambda _c}\hbox {Li}\) to \(^{209}_{\varLambda _c}\hbox {Bi}\). For convenience the corresponding core nucleus is indicated in brackets. Results are shown for two values of the cutoff, \(\varLambda =500,\, 600\) MeV

It is interesting to observe that, contrary to single-\(\varLambda \) hypernuclei where the \(\varLambda \) is more and more bound when going from light to heavy nuclei, the binding energy of the \(\varLambda _c\) increases from \(^{\ 5}_{\varLambda _c}\hbox {Li}\) to \(^{41}_{\varLambda _c}\hbox {Sc}\) and then it decreases. This is due to the Coulomb repulsion between the \(\varLambda _c\) and the protons of the nuclear core, which together with the kinetic energy of the \(\varLambda _c\), compensates most of the attraction of the \(\varLambda _cN\) interaction. The possible existence of \(\varLambda _c\) nuclei is, therefore, subject to a delicate balance between the \(\varLambda _cN\) interaction, the kinetic energy and the Coulomb force as it has been already pointed out in Refs. [17, 33, 34, 49, 91]. In particular, in Ref. [49] it was suggested that only light- or medium-mass \(\varLambda _c\) nuclei could really exists whereas, for instance, in Ref. [17] it was found that even the heavier \(\varLambda _c\) nucleus considered in that work, namely \(^{209}_{\varLambda _c}\hbox {Bi}\), could exist, as in the present work. A small spin-orbit splitting of the \(p-\) and \(d-\)wave states of the order of a few tenths of MeV is observed in all \(\varLambda _c\) nuclei in agreement with the results obtained in Refs. [17, 33,34,35, 91]. In addition, we note also that the level spacing of the \(\varLambda _c\) single-particle states is smaller than those for the corresponding hypernuclei (see e.g. Table I of Ref. [86]). This is simply due to the fact that the mass of the \(\varLambda _c\) is larger than that of the \(\varLambda \) hyperon.

Fig. 3
figure 2

Contributions of the kinetic energy, the \(Y_cN\) interaction and the Coulomb potential to the energy of the \(\varLambda _c\) single-particle bound state \(1s_{1/2}\) as a function of the mass number of the \(\varLambda _c\) nuclei considered. Results are shown for the potential \(Y_c N\)-A with a cutoff of 500 MeV

To understand better the role of the Coulomb force in our calculation, in Fig. 2 we show the separate contributions of the kinetic energy, of the \(Y_c N\) interaction, and of the Coulomb potential to the energy of the \(\varLambda _c\) single-particle bound state \(1s_{1/2}\) for the different charmed nuclei considered in this work as function of the mass number (\(A=N+Z\), with N and Z being the neutron and atomic numbers, respectively, of the specific nucleus). When going from light to heavy \(\varLambda _c\) nuclei, the Coulomb contribution increases because of the increase of the atomic number whereas those of the kinetic energy and of the \(Y_cN\) interaction decrease. The contribution of the kinetic energy decreases with the mass number because the wave function of the \(1s_{1/2}\) state becomes more and more spread due to the larger extension of the nuclear density over which the \(\varLambda _c\) wants to be distributed (see Fig. 3). The increase of the mass number leads to a more attractive \(\varLambda _c\) self-energy (see, e.g., Figs. 2 and 3 of Ref. [86] for a detailed discussion in the case of single-\(\varLambda \) hypernuclei) that translates into a more negative contribution of the \(Y_cN\) interaction. Note that, when adding the three contributions they compensate in such a way that the energy of the \(1s_{1/2}\) decreases only by about 5 MeV from \(^{\ 5}_{\varLambda _c}\hbox {Li}\) to \(^{17}_{\varLambda _c}\hbox {F}\) and then it increases very smoothly from \(^{41}_{\varLambda _c}\hbox {Sc}\) to \(^{209}_{\varLambda _c}\hbox {Bi}\).

Fig. 4
figure 3

\({\varLambda }_c\) probability density distribution for the \(1s_{1/2}\) state in the six \({\varLambda }_c\) nuclei considered. Results are presented for the \(Y_c N\)-A potential. The variation with the cutoff is indicated by bands. The red bands show the results when the Coulomb interaction is artificially switched off. Note that different scales are used for the heavier nuclei

To end this section, we display in Fig. 3 the probability density distribution (i.e., the square of the radial wave function) of the \(\varLambda _c\) in the \(1s_{1/2}\) state for the six \(\varLambda _c\) nuclei considered, for the \(Y_c N\)-A potential. The variation with the cutoff is indicated by bands. Results for the \({\varLambda }_cN\) interaction from Ref. [51] and for \(Y_c N\)-B are not shown since the differences in the probability density are so small that they cannot be resolved in the plot. Note that, when moving from light to heavy \(\varLambda _c\) nuclei, due to the increase of the size of the nuclear core, the probability of finding the \(\varLambda _c\) close to the center of the nucleus decreases (notice the different scales of the panels), and it becomes more and more distributed over the whole nucleus. The probability density distribution when the Coulomb interaction is artificially switched off is also shown for comparison. Obviously, and as expected, the Coulomb repulsion pushes the \(\varLambda _c\) away from the center of the nuclei. A similar effect is observed for the probability densities of the other \(\varLambda _c\) single-particle bound states.

4 Summary and conclusions

In the present work, we have investigated the binding energies of charmed nuclei. As input we used \({\varLambda }_cN\) and \({\varSigma }_cN\) interactions that have been extrapolated from lattice QCD simulations by the HAL QCD collaboration [49, 50] at quark masses corresponding to \(m_\pi \) = 410–570 MeV to the physical point. For this extrapolation, we used a framework based on chiral effective field theory [51,52,53]. The \({\varLambda }_cN\) interaction established in this way is significantly weaker than what has been employed in most of the studies of charmed nuclei in the literature so far. The bound state calculations for light charmed nuclei have been carried out within the Faddeev–Yakubovsky framework. The results for heavier nuclei are from calculations of the energies of \({\varLambda }_c\) single-particle bound states, performed within a perturbative many-body approach, which allows one to determine the finite nuclei \({\varLambda }_c\) self-energy from which the energies of the different bound states can be obtained.

Our results indicate that even for a weak \({\varLambda }_cN\) interaction as suggested by the lattice simulations of the HAL QCD collaboration already \(A=4\) charmed nuclei are likely to exist. Only the lightest nucleus considered, a charmed helium \(^{\ 3}_{{\varLambda }_c}\hbox {He}\), turned out to be unbound, in contrast to conjectures reported in Refs. [15, 16].

An additional aspect considered in the present work is the effect of the \({\varSigma }_cN\) interaction. Some results from lattice simulations for this channel have become available recently [50]. There is admittedly a sizable uncertainty in the extrapolation of the HAL QCD results to the physical point, not least due to missing information on the behavior in the \({\varSigma }_c^* N\) channel, closely connected to the former by heavy quark spin symmetry. This makes reliable predictions for \({\varSigma }_cN\) observables rather difficult at the moment. On the other hand, we found that the uncertainties due to the present situation in the \({\varSigma }_cN\) channel do not affect the conclusions on the properties of the \({\varLambda }_cN\) interaction at low energies, relevant for the quest of charmed nuclei. Specifically, taking into account the coupling of \({\varLambda }_cN\) to \({\varSigma }_cN\) and the direct \({\varSigma }_cN\) interaction as suggested by the HAL QCD results has very little influence on the existence of such bound states. Indeed, for the \({\varLambda }_cN\) interaction, the extrapolation of the lattice results to the physical point seems to be fairly reliable and stable and, therefore, we believe that robust predictions for the properties of the \({\varLambda }_c\) in finite and infinite nuclear matter can be given based on the \({\varLambda }_cN\) potentials established in Ref. [51] and in this work.

Prospects for detecting charmed nuclei at J-PARC have been discussed at various occasions, see, e.g., Ref. [92]. Corresponding opportunities by the CBM experiment at FAIR are considered in Ref. [93]. The option for discovering charmed nuclei with neutrino beams is addressed in Ref. [94]. An alternative on a different scope is offered by high-energy experiments such as the pp- and/or heavy-ion collisions [95] presently pursued by the ALICE collaboration [96, 97] at the LHC/CERN or the STAR collaboration at RHIC/BNL [98, 99]. Here “exotic” nuclei such as the anti-hypertriton or the \(^4\overline{\mathrm{He}}\) were already produced and detected, and the lightest charmed nuclei might be within reach—now or in the near future—should they indeed exist. That said, one should be aware that there are tremendous experimental challenges for producing and detecting charmed nuclei, as has been summarized in Ref. [17] but also indicated in the introduction to the present work.