Predictions for charmed nuclei based on $Y_c N$ forces inferred from lattice QCD simulations

Charmed nuclei are investigated utilizing $\Lambda_c N$ and $\Sigma_c N$ interactions that have been extrapolated from lattice QCD simulations at unphysical masses of $m_\pi = 410$--$570$ MeV to the physical point using chiral effective field theory as guideline. Calculations of the energies of $\Lambda_c$ single-particle bound states for various charmed nuclei from $^{\ 5}_{\Lambda_c}$Li to $^{209}_{\Lambda_c}$Bi are performed using a perturbative many-body approach. This approach allows one to determine the finite nuclei $\Lambda_c$ self-energy from which the energies of the different bound states can be obtained. Though the $\Lambda_c N$ interaction inferred from the lattice results is only moderately attractive, it supports the existence of charmed nuclei. Already the lightest nucleus considered is found to be bound. The spin-orbit splitting of the p- and d-wave states turns out to be small, as in the case of single $\Lambda$ hypernuclei. Additional calculations based on the Faddeev-Yakubovsky equations suggest that also $A=4$ systems involving a $\Lambda_c$ baryon are likely to be bound, but exclude a bound $^{\, 3}_{\Lambda_c}$He state.


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 Λ 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 Λ c [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35] that started soon after the discovery of charmed baryons [36,37] (see also the recent reviews [38,39]). In principle, charmed nuclei could be produced by means of charm exchange or associate charm production reactions [40,41], 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 [42][43][44][45].
The Λ c N 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 [46] 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 [47,48]. Those studies, based on unphysical quark masses corresponding to pion masses of m π = 410 − 700 MeV, suggest that the Λ c N and Σ c N interactions could be significantly less attractive than what had been proposed in the phenomenological studies mentioned above. In Ref. [49], an extrapolation of the HAL QCD results to the physical point was presented, using chiral effective field theory (EFT) [50,51] as guideline. It revealed that the Λ c N 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 Λ c N interaction.
In the present work, we use the Λ c N interaction from Ref. [49] to explore the binding energies of charmed nuclei. In the literature different conventions for naming Λ c nuclei have been used in the past. We adopt here the standard nomenclature for nuclei with the proper gener-alization [29], spelled out for hypernuclei in Sect. I.B of Ref. [52]. 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) Λ to hydrogen ( 2 H) one gets the hypertriton ( 3 Λ H) but adding the positively charged Λ c leads to 3 Λc He. Similarly, based on this convention, the bound state of Λ c and 208 Pb is 209 Λc Bi. The lightest nuclei considered, the A = 3 and A = 4 systems 3 Λc He and 4 Λc 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 Λ c N interactions employed in the former works are so strongly attractive that they even predict two-body bound states (in the 1 S 0 as well as in the 3 S 1 partial wave) with binding energies comparable to that of the deuteron. For calculating heavier charmed nuclei, namely from 5 Λc Li onward to 209 Λc Bi, a perturbative many-body approach is utilized that allows one to obtain the Λ c single-particle bound states in the different nuclei from the corresponding Λ 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 [53] of the hyperon-nucleon (Y N ) meson-exchange potential of the Jülich Group [54] via SU(4) symmetry arguments. For these interactions, the 5 Λc Li system ( 5 Λc He in the nomenclature used in [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 [47]. The calculations were performed with the Λ c N potentials extracted from the lattice simulations at pion masses 410-700 MeV, but using the physical masses of the Λ c and the considered nuclei. Binding energies for Λ c bound to 12 C, 28 Si, 40 Ca, 58 Na, 90 Zr, and 208 Pb were reported.
Since after the publication of the Λ c N interaction [49] new results from the HAL QCD collaboration became available that include now the Σ c N channel [48], we also revisit the Λ c N interaction in order to explore in how far the inclusion of a direct interaction in the Σ c N channel modifies the extrapolation of the Λ c N 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 Σ c N channel has very little influence on the Λ c N scattering results and also not on those for Λ 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 Σ c N S-wave phase shifts based on our extrapolation of the HAL QCD results [48].
Another extrapolation of the lattice results for Σ c N , performed in heavy baryon chiral perturbation theory and taking into account heavy quark spin symmetry, has been performed recently [55]. Let us mention already now that neither in that work nor in our calculation any signal for resonances or bound states in the Σ c N channel are found.
The existence of such resonances has been suggested in Ref. [56], where the Λ c N -Σ c N -Σ * c N interaction was investigated within the framework of the quark delocalization color screening model. Bound states (and resonances) in the Σ c N 3 S 1 partial wave around the Σ c N and Σ * c N thresholds were also predicted in Ref. [57] based on a Y c N interaction from meson exchange supplemented by shortrange repulsion from a quark exchange model.
The paper is structured in the following way: A summary of the main characteristics of the Λ c N and Σ c N potentials is presented in Sec. 2 . Results of the properties of Λ c in infinite nuclear matter, and light and heavier nuclei are reported in Sec. 3. Finally, a brief summary and some concluding remarks are given in Sec. 4. The appendix summarizes results for Σ c N scattering.

The Λ c N and Σ c N potentials
The Λ c N and Σ c N interactions are constructed by using chiral EFT as guideline, following the scheme employed in our studies of the ΛN and ΣN systems [58][59][60][61]. We summarize the essentials below. More details of the approach can be found in Ref. [49]. The Y c N potential consists of contact terms and contributions from pion exchange. The former are given by for the partial waves considered in the present study. Here p = |p | and p = |p | are the initial and final center-ofmass (c.m.) momenta in the Λ c N or Σ c N systems. The quantitiesC i ,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 π = 410 MeV and 570 MeV. The m π dependence in Eq. (1) is motivated by the corresponding expression in the standard Weinberg counting up to next-to-leading order (NLO) [50,51] but differs from it by the term proportional to m 2 π (p 2 +p 2 ) which is nominally of higher order. Nevertheless, we included that term in Ref. [49] because it allowed us to obtain an optimal description of the LQCD results at m π = 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 π = 410 MeV would be underestimated at low energies, as exemplified by results shown in Figs. 6 and 7 of Ref. [55].
The contribution of pion exchange to the Y c N potential is given by where q is the transferred momentum, q = p − p. The coupling constants f BB π are related to the axial-vector strength via f BB π = g BB A /2F π with F π being the pion decay constant (F π ≈ 93 MeV). The coupling constant for the Λ c Σ c π vertex can be determined from the experimentally known Σ c → Λ c π decay rate, see Refs. [62,63]. For the Σ c Σ c π coupling constant lattice QCD results [64] are employed. To be concrete, g ΣcΣc A = 0.71 [64] and g ΛcΣc A = 0.74 [62,63] are used, together with g N N A = 1.27 [65]. Note that f ΛcΛcπ ≡ 0 under the assumption that isospin is conserved.
Besides the coupling constants at the physical point, one needs also their m π dependence: Results for the dependence of F π on m 2 π are available from lattice simulations [66]. Based on that work the values F π ≈ 112 MeV at m π = 410 MeV and F π ≈ 129 MeV at m π = 570 MeV were deduced and employed in Ref. [49]. With regard to the dependence of g BB A on m 2 π lattice QCD simulations indicate a rather small variation, at least for g ΣcΣc A and g N N A where concrete results are available [64,67]. Because of that the dependence of the g A 's on m π was neglected in Ref. [49] and the values at the physical point were used throughout.
There was no information from LQCD on the Σ c N interaction at the time when the study in Ref. [49] was performed, and, thus, the interaction in the Σ c N channel was not considered. Nonetheless, the coupling of Λ c N to Σ c N via pion exchange was already included. Due to its long-range nature, this coupling plays an important role in case of the ΛN and ΣN systems [59,60]. Since f ΛcΣcπ is empirically known and only slightly smaller than f ΛΣπ [59] the coupling between Λ c N and Σ c N via pion exchange should still be of relevance. Indeed, the effective contribution to the Λ c N interaction, V ΛcN →ΛcN ∼ V OP E ΛcN →ΣcN G ΣcN V OP E ΣcN →ΛcN , is expected to be smaller for energies around the Λ c N threshold as compared to the situation for ΛN , but just by a factor 2 − 3. The reduction is due to the larger mass difference, M Σc − M Λc ≈ 167 MeV versus M Σ − M Λ ≈ 78 MeV, that enters in the corresponding Green's functions G ΣcN or G ΣN . In any case, formally two-pion exchange contributions to the Λ c N interaction involving the Σ c do arise at NLO [59], and it can be expected that the piece with a Σ c N intermediate state provides the dominant contribution [68]. In the actual calculation, this (reducible) two-pion exchange contribution to the Λ c N potential is generated by solving a coupled-channel Lippmann-Schwinger (LS) equation (see below). Further (irreducible) NLO contributions from twopion exchange [59] have been omitted in [49] for simplicity reasons. It was assumed that those can be effectively absorbed into the contact terms.
For the present study, we add a direct Σ c N interaction. The potential is determined in the same way as the one for Λ c N , by using results of the HAL QCD collaboration for the corresponding 3 S 1 phase shift at m π = 410 MeV and 570 MeV [48]. The main goal of this extension is to explore in how far the inclusion of a direct Σ c N interaction modifies the Λ c N 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 Σ c N 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 [58], the equation is given by where the label ρ indicates the channels (Λ c N , Σ c N ) and the label ν the partial wave. The quantity µ ρ signifies the pertinent reduced mass. The on-shell momentum in the intermediate state, Since the integral in the LS equation (4) is divergent for the chiral potentials specified above, a regularization scheme needs to be introduced [69,70]. For that purpose the potentials in the LS equation are cut off in momentum space by multiplication with a regulator function [58,59], f (p , p) = exp − p 4 + p 4 /Λ 4 . In Ref. [49], cutoff values Λ = 500-600 MeV were employed, in line with the range that yielded the best results in NLO studies of the ΛN and ΣN interactions [59,60]. 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 π = 410 and 570 MeV are taken from Ref. [47].

Λ c nuclei and matter properties
In this section, we report results on the properties of the Λ 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 [47,48] to the physical point. In particular, we use the Λ c N potential from Ref. [49], where there is no direct Σ c N interaction, and the two potentials Y c N -A and Y c N -B, introduced and described in the appendix, which include a direct Σ c N interaction.

Λ c in infinite nuclear matter
In order to investigate the properties of the Λ c N 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 ΛN interaction [71]. In that work and similar ones (see, e.g., Refs. [72,73]), the reader can find details how to solve the Bethe-Goldstone equation and how the singleparticle (s.p.) potential U Λc is determined self-consistently together with the G-matrices for a specific nuclear matter density ρ (or Fermi momentum k F ).
In Fig. 1(a), we present results for the dependence of U Λc (k Λc = 0) on the Fermi momentum, in comparison to those for the Λ hyperon obtained with the NLO interaction from Refs. [59,60]. In Fig. 1(b), we display the dependence of U Λc (k Λc ) and the s.p. energy, ε Λc (k Λc ) = k 2 Λc /2m Λc + U Λc (k Λc ), on the Λ c momentum at the Fermi momentum k F = 1.35 fm −1 , i.e., at nuclear matter saturation density. The in-medium predictions for Λ c are based on the Y c N -A potential (cf. appendix) and the Λ c N potential from Ref. [49]. Results for the Λ c properties for the alternative fit Y c N -B, considered in the appendix, practically coincide with the ones for Y c N -A and are therefore not shown. Partial-wave contributions to U Λc (k Λc = 0) at k F = 1.35 fm −1 are listed in Table 1. Note that the contributions of the P waves come solely from two-pion exchange involving the Σ c N intermediate state. The total potential depth amounts to around −20 ÷ −18 MeV and is quite insensitive to whether a direct Σ c N interaction is included or not. As a reminder, the "empirical" value in case of the Λ hyperon is −30 ÷ −27 MeV [52]. Comparing the Λ c results with the Λ case in detail, one can see that the contribution in the 1 S 0 partial wave is reduced by roughly a factor three. This is well in line with the corresponding interaction strengths; the Λ c N scattering length is also about a factor three smaller than that for ΛN , see Table 1. For the 3 S 1 partial wave the Λ c N and ΛN contributions (for NLO13 [59]) are of comparable magnitude, despite of the fact that the Λ c N interaction is less attractive as reflected in the corresponding scattering lengths which is about 30-50 % smaller than that for ΛN scattering. Obviously, for Λ c N the dispersive effects, which play an important role for the contribution of that partial wave in case of the Λ [60], are smaller because the Λ c N -Σ c N 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 Λ c N interaction. When comparing with the results for the NLO interaction from 2019 [60], where the ΛN -ΣN transition potential is noticeably weaker, one sees a clear correlation between the smaller 3 S 1 scattering length and the reduced contribution to U Λc (cf. Table 1).
Finally, let us compare our nuclear matter results with other predictions for U Λc found in the literature. Ref. [30] contains some results for U Λc based on an Y c N potential that is adapted from one of the Y N meson-exchange potentials by the Nijmegen Group by imposing SU(4) flavor symmetry. In that work, a value of U Λc (k Λc = 0) ≈ −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 ∼ −20 MeV for Λ c at nuclear matter saturation density. Yasui, in a perturbative approach based on a heavy-quark effective theory, finds a Λ c binding energy of around −25 ÷ −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
Λc He and 4 Λ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 [60,74,75]. As discussed in Ref. [49], based on the results for Λ-hypernuclei and the relative strengths of the Λ c N and ΛN interactions, one can guess which light charmed nuclei could be bound. For that the interaction pertinent mixtures of the spin-singlet and spin-triplet Λ c N (ΛN ) interaction for s-shell nuclei are relevant [60,76] and, of course, the reduction of the kinetic energy associated with the Λ c as a consequence of its larger mass [29]. In view of the fact that, for the considered Y c N interactions, the Λ c N 1 S 0 scattering length is only one third of the one for ΛN , while there is somewhat less difference in the 3 S 1 state (cf. Additionally, whereas the Coulomb interaction is of less importance for Λ separation energies of hypernuclei, its contribution to Λ c separation energies of charmed nuclei is often decisive for binding [29]. For the solution of the Faddeev-Yakubovsky equations here, we take the Coulomb interactions fully into account as described in [77]. First, in order to benchmark our few-body calculations, we devised Λ c N interactions that mimic the effective range parameters predicted by a Λ c N potential obtained in the constituent quark model by Garcilazo et al. [46]. In that model, the triplet interaction is much stronger (a3 S1 = −2.31 fm) than the one in the singlet channel (a1 S0 = −0.86 fm), cf. Table 3 in that reference. Consequently, and in line with the above arguments and the explorations in Ref. [74], a bound state for the J = 3 2 state of 3 Λc He has been reported, with a Λ 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 N N interactions. The uncertainty due to different cutoffs of the Λ c N interaction is larger than the one due to different N N interactions. We also found that no bound state ex-ists for that interaction for 3 Λc He with J = 1 2 . This result confirms the earlier calculations of Garcilazo et al.
We then performed calculations for the Λ c N potentials from Ref. [49], 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 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. [29] has Λ c N 1 S 0 and 3 S 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 Λc He with a Λ c separation energy of approximately 1.25 MeV was predicted, after including an estimate for the contribution of the Coulomb interaction.
The A = 4 results for the interactions considered in this 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. As obvious from Table 2, the results are very independent on whether the direct Σ c N 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 Σ c N 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 [29]. We believe that this is partly due to omitting tensor interactions in the Y c N in [29] 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. [46]. This interaction clearly provides stronger binding leading to 1.2 to 2.1 MeV separation energy depending on the cutoff used. Table 2. Λc separation energies EΛ c and binding energies with respect to breakup up into four baryons E for the J = 1 + state of 4 Λc He. The N N interaction at order N 4 LO+ with a cutoff of 450 MeV of Ref. [78] was used leading to E( 3 H) = −8.141 MeV. Expectation values for the kinetic energy T , the N N potential energy VNN , and the YcN potential energy VY c N are also given. The probability that the four-baryon have orbital angular momentum zero and two ( P (S) and P (D) ) are listed together with the probability P (Σc) for the Σc component. Energies are given in MeV and probabilities in %.
ΛcN (500) [49]  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 N N potential energy and the Y c N potential energy. It sticks out that the N N 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 dependent on the Σ c N contribution as can be seen from the similarity of the results of Λ c N [49], 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 N N and Y c N interactions induce a D-wave contribution of approximately 7%. The probability to find a Σ c is, similar to the Σ 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 Λc 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 Λ 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. [46]. In that case the separation energy is between 130 keV and 470 keV depending on the cutoff used.
Given that the dependence on the N N interaction was smaller than the dependence on 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 N N interactions. We do not expect that the results will be significantly different for other choices.

Heavier charmed nuclei
Now we consider the energy of the Λ 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 Sec. 2 and the appendix. This G-matrix is then used to calculate the self-energy of the Λ 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 Λ c in the nucleus. This approach also provides the real and imaginary parts of the Λ c optical potential at positive energies and, therefore, allows one to study the Λ c -nucleus scattering properties. This method was already used to study the properties of the nucleon [79], the ∆ isobar [80] and the Λ and Σ hyperons [81][82][83] in finite nuclei, and very recently also those of the Λ c using a meson-exchange Y c N 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.  Table 4 for the Λ c N interaction from Ref. [49] and the two potentials Y c N -A and Y c N -B with inclusion of a direct Σ c N interaction. We note that all charmed nuclei considered consist of a closed-shell nuclear core plus a Λ c sitting in a single-particle state. We note also that, although the Λ c N interaction with the cutoff 600 MeV is more attractive, the Λ c single-particle bound states predicted in this case are actually less bound. This is due to dispersive effects [84][85][86] in the calculation of the G-matrix which suppress the contribution from the Λ c N → Σ c N → Λ c N 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. [88], the approach followed tends to underestimate the energies of the Λ hyperon single-particle bound states for light hypernuclei such as 5 Λ He. Accordingly, we expect 5 Λc Li to be somewhat more strongly bound than what is suggested by the values given in Table 4.
It is interesting to observe that, contrary to single-Λ hypernuclei where the Λ is more and more bound when going from light to heavy nuclei, the binding energy of the Λ c increases from 5 Λc Li to 41 Λc Sc and then it decreases. This is due to the Coulomb repulsion between the Λ c and the protons of the nuclear core, which together with the kinetic energy of the Λ c , compensates most of the attraction of the Λ c N interaction. The possible existence of Λ c nuclei is, therefore, subject to a delicate balance between the Λ c N interaction, the kinetic energy and the Coulomb force as it has been already pointed out in Refs. [17,32,33,47,87]. In particular, in Ref. [47] it was suggested that only lightor medium-mass Λ c nuclei could really exists whereas, for instance, in Ref. [17] it was found that even the heavier Λ c nucleus considered in that work, namely 209 Λc 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 Λ c nuclei in agreement with the results obtained in Refs. [17,[32][33][34]87]. In addition, we note also that the level spacing of the Λ c single-particle states is smaller than those for the corresponding hypernuclei (see e.g. Table I of Ref. [83]). This is simply due to the fact that the mass of the Λ c is larger than that of the Λ hyperon.
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 Λ c singleparticle 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 Λ c nuclei, the Coulomb contribution increases because of the increase of the atomic number whereas those of the kinetic energy and of the Y c N 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 Λ c wants to be distributed (see Fig. 3). The increase of  the mass number leads to a more attractive Λ c self-energy (see, e.g., Figs. 2 and 3 of Ref. [83] for a detailed discussion in the case of single-Λ hypernuclei) that translates into a more negative contribution of the Y c N 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 Λc Li to 17 Λc F and then it increases very smoothly from 41 Λc Sc to 209 Λc Bi. To end this section, we display in Fig. 3 the probability density distribution (i.e., the square of the radial wave function) of the Λ c in the 1s 1/2 state for the six Λ c nuclei considered, for the Y c N -A. The cut-off dependence is indicated by bands. Results for the Λ c N interaction from Ref. [49] 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 Λ c nuclei, due to the increase of the size of the nuclear core, the probability of finding the Λ c close to the center of the nucleus decreases, 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 Λ c away from the center of the nuclei. A similar effect is observed for the probability densities of the other Λ c single-particle bound states.

Summary and Conclusions
In the present work, we have investigated the binding energies of charmed nuclei. As input we used Λ c N and Σ c N interactions that have been extrapolated from lattice QCD simulations by the HAL QCD collaboration [47,48] at quark masses corresponding to m π = 410−570 MeV to the physical point. For this extrapolation, we used a framework based on chiral effective field theory [49][50][51]. The Λ c N 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 Λ c single-particle bound states, performed within a perturbative many-body approach, which allows one to determine the finite nuclei Λ c self-energy from which the energies of the different bound states can be obtained.
Our results indicate that even for a weak Λ c N 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 Λc 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 Σ c N interaction. Some results from lattice simulations for this channel have become available recently [48]. 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 Σ * c N channel, closely connected to the former by heavy quark spin symmetry. This makes reliable predictions for Σ c N observables rather difficult at the moment. On the other hand, we found that the uncertainties due to the present situation in the Σ c N channel do not affect the conclusions on the properties of the Λ c N interaction at low energies, relevant for the quest of charmed nuclei. Specifically, taking into account the coupling of Λ c N to Σ c N and the direct Σ c N interaction as suggested by the HAL QCD results has very little influence on the existence of such bound states. Indeed, for the Λ c N 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 Λ c in finite and infinite nuclear matter can be given based on the Λ c N potentials established in Ref. [49] and in this work.
Prospects for detecting charmed nuclei at J-PARC have been discussed at various occasions, see, e.g., Ref. [89]. Corresponding opportunities by the CBM experiment at FAIR are considered in Ref. [90]. The option for discovering charmed nuclei with neutrino beams is addressed in Ref. [91]. An alternative on a different scope is offered by high-energy experiments such as the pp-and/or heavy-ion collisions [92] presently pursued by the ALICE collaboration [93,94] at the LHC/CERN or the STAR collaboration at RHIC/BNL [95,96]. Here "exotic" nuclei such as the anti-hypertriton or the 4 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.
We acknowledge helpful communications with Vadim Baru, Benjamin Dönigus, and Hirokazu Tamura. This work is sup-ported in part by the DFG and the NSFC through funds provided to the Sino-German CRC 110 "Symmetries and the Emergence of Structure in QCD" (DFG grant. no. TRR 110), and the COST Action CA16214. The numerical calculations were performed on JURECA, the JURECA-Booster of the Jülich Supercomputing Centre, Jülich, Germany, and the Centro di Calcolo of the INFN Sezione di Catania, Catania, Italy.

Appendix: Σ c N scattering
In this appendix, we discuss the results for Σ c N scattering. However, let us emphasize from the beginning that these have to be interpreted with caution. For charmed baryons, heavy quark spin symmetry plays a role [55,97] and, thus, one should include not only the Σ c N channel but also Σ * c N . Indeed the Σ c N and Σ * c N thresholds are just about 65 MeV apart so that the coupling between those systems should be important. Unfortunately, there are no results for the Σ * c N interaction from lattice simulations and, therefore, it cannot be explicitly included in the analysis. For that reason, it was omitted in our earlier study [49] and it was assumed that any effect of the Σ * c N interaction can be effectively absorbed into the Λ c N LECs. After all, since M Σ * c − M Λc = 234 MeV, there should be little influence on the low-energy Λ c N amplitude anyway. In case of the Σ c N channel, such an assumption is questionable. In Ref. [55] heavy quark spin symmetry was taken into account in the derivation of the potential. But also in that work, the actual coupling of the Σ c N and Σ * c N channels was ignored in the evaluation of the scattering amplitude.
Another complicacy comes from the fact that the mass difference M Σc − M Λc is larger than the pion mass. Because of that the use of the static approximation for the pion exchange (see Eq. (2)) in the Λ c N → Σ c N transition potential is rather problematic for energies near and/or above the Σ c N threshold where the exchanged pion can go on-shell. In principle, one should take into account that there is a Λ c N π three-body cut, see Refs. [98][99][100] for a discussion on this issue for the analogous DD * -DDπ case in the context of the X(3872). As a matter of fact, for the unphysical pion masses of the LQCD calculation, this problem does not arise. However, it could have a noticeable influence on the extrapolation of the Σ c N results to the physical point. This aspect is ignored in our calculation, and it is also not addressed in the Σ c N study of Meng et al. [55]. We emphasize that, for energies around the Λ c N threshold, the static approximation is well justified.
LQCD results for the Σ c N 3 S 1 phase shifts are available for m π = 410, 570, 700 MeV [48]. We determine the LECs of the contact interaction, cf. Eq. (1), by a fit to the lattice data at the two lower pion masses. It allows us to determine the LECsC i and C i as well asD i and D i , i.e., the ones which encode the pion-mass dependence of the contact interaction. The fits are done to the central values of the phase shifts as given in Fig. 2 of Ref. [48], for energies up to 30 MeV, cf. upper panel of Fig. 4. Alternative fits with particular emphasis on the near-threshold  behavior of the HAL QCD results were performed, too, cf. the lower panel. Of course, in both cases, we made sure that we produce larger near-threshold phase shifts for m π = 410 MeV than for 570 MeV, as suggested by the lattice simulation. It should be said that trying to fit to the HAL QCD results at somewhat higher energies is not very meaningful in view of the fact that the Σ * c N channel is not explicitly included. Its threshold is around 85 MeV for the HAL QCD calculations [48] and, as said above, at roughly 65 MeV for physical masses.
A combined fit to the Λ c N and Σ c N 3 S 1 phase shifts turned out to be unnecessary because the inclusion of a direct Σ c N interaction had practically no effect on the Λ c N results reported in [49], at least for the energies considered there. We did not attempt to reproduce the inelasticity parameter, given in Ref. [48] in terms of the Table 5. ΛcN and ΣcN (I = 1/2) scattering lengths in the 3 S1 partial wave (in fm), for the YcN potentials described in the text. In case of ΣcN real and imaginary parts are given. ΛcN results for the interaction from Ref. [49] are listed as well.

ΛcN ΣcN
ΛcN (500)  Λ c N S-matrix, η = |S ΛcN,ΛcN |. In the lattice simulation it is with values around 0.995 basically compatible with 1, which suggests that there is practically no channel coupling. However, we believe that this could be an artifact of the way how the analysis by the HAL QCD collaboration is done. In the ΛN − ΣN systems, the strong channel coupling arises primarily from the tensor force mediated by pion exchange, and that leads to a strong coupling of the 3 S 1 to the 3 D 1 partial wave, see, e.g., Fig. 7 in Ref. [59]. Indeed, for 3 S 1 → 3 S 1 transitions (as well as for 1 S 0 → 1 S 0 ), the expectation value of the tensor operator is zero. However, in the analysis of the HAL QCD collaboration, D waves are not considered so that this component of the Y c N force is only effectively included in the S-wave interactions. In our calculation, we include the full one-pion exchange and that means that there is a coupling to the 3 D 1 partial wave -at the physical point and also for the pion masses of the lattice simulation. Since the pion mass is known and the Λ c Σ c π coupling constant is known, too (at least at the physical point), one can consider our result as genuine prediction for the strength of the channel coupling and, thus, we did not impose any further constraints on it. Results for the 3 S 1 phase shift are presented in Fig Fig. 4. If we require a quantitative reproduction of the low-energy behavior (Y c N -B; lower panel), then there is agreement with the lattice results only up to around 25 MeV. The cutoff dependence of the fits to the HAL QCD results is fairly small so that the (hatched) bands are barely visible. The phase shift obtained from the interaction at the physical point (solid bands) do exhibit a noticeable but still moderate cutoff dependence. That said, there is a sizable uncertainty in the extrapolation of the two fitting scenarios considered. When emphasis is put on the very low-energy results by HAL QCD, then the phase shifts at the physical point are larger and actu-  Fig. 5. Results for the ΛcN 3 S1 phase shift for the interaction YcN -A. Same description of curves and symbols as in Fig. 4. Lattice QCD results are from Ref. [47].
ally close to the lattice results for unphysical pion masses ( Fig. 4; lower panel) whereas the fit over a larger energy region yields perceptibly smaller results.
For convenience, we compiled the scattering lengths for Y c N -A and B in Table. 5. Obviously, with the Σ c N interaction included there is a small reduction in the Λ c N 3 S 1 scattering length (for A and for B) as conpared to the results from Ref. [49]. However, the variation is rather moderate and stays within the uncertainty due to the regulator dependence. Also the variation in the corresponding Λ c N 3 S 1 phase shift is very small. Actually, the Λ c N results with inclusion of a direct Σ c N interaction, cf. Fig. 5 for Y c N -A, can be hardly distinguished from those presented in Fig. 2 in Ref. [49]. For the Σ c N channel the situation is rather different. First there is a sizable difference in the scattering length obtained for the interactions A and B. In addition, and more disturbing, there is a fairly drastic regulator dependence. One can certainly say that there is no indication for a near-by Σ c N bound state, a conclusion already drawn in Ref. [55]. (The existence of Σ c N bound states or resonances has been suggested in some studies in the past [56,57]). On the other hand, more quantitative conclusions are difficult to draw. Definitely, there is a significant coupling to the Λ c N channel, mediated by one-pion exchange, which gives rise to an appreciable imaginary part of the scattering length in our calculation. Moreover, one should not forget that the static approximation is used by us for simplicity reasons. In the real world another channel is open, namely Λ c N π, which contributes likewise to the inelasticity. An appropriate treatment is desirable but technically demanding and, thus, postponed to the future when hopefully more information from lattice simulations will be available. Fortunately, these uncertainties have basically no effect on the predictions for the properties of the Λ c N interaction at low energies, cf. Table. 5 and Fig. 5, and also not on the binding energies of Λ c nuclei, as discussed in the main part of this paper.