Doubly magic nuclei from Lattice QCD forces at $M_{PS}=$469 MeV/c$^2$

TL;DR

This study uses ab initio Green's function calculations with Lattice QCD derived potentials at a specific pion mass to explore the structure and stability of closed shell nuclei, revealing improved energies and potential stability islands.

Contribution

First application of self-consistent Green's function methods with Lattice QCD potentials to calculate properties of medium-mass nuclei at unphysical pion mass.

Findings

Ground state energies are improved over previous methods.

Nuclei may be unbound at high pion masses and stable at lower masses.

Spectral distributions resemble those of real nuclei.

Abstract

We perform ab initio self-consistent Green's function calculations of the closed shell nuclei $^{\rm 4}$He, $^{\rm 16}$O and $^{\rm 40}$Ca, based on two-nucleon potentials derived from Lattice QCD simulations, in the flavor SU(3) limit and at the pseudo-scalar meson mass of 469~MeV/c$^{\rm 2}$. The nucleon-nucleon interaction is obtained using the HAL QCD method and its short-distance repulsion is treated by means of ladder resummations outside the model space. Our results show that this approach diagonalises ultraviolet degrees of freedom correctly. Therefore, ground state energies can be obtained from infrared extrapolations even for the relatively hard potentials of HAL QCD. Comparing to previous Brueckner Hartree-Fock calculations, the total binding energies are sensibly improved by the full account of many-body correlations. The results suggest an interesting possible behaviour in which nuclei are unbound at very large pion masses and islands of stability appear at first around the traditional doubly-magic numbers when the pion mass is lowered toward its physical value. The calculated one-nucleon spectral distributions are qualitatively close to those of real nuclei even for the pseudo-scalar meson mass considered here.

Tables

Table 1. Table 1: Ground state energies of 4 He, 16 O and 40 Ca at M P ​ S subscript 𝑀 𝑃 𝑆 M_{PS} =469 MeV/c 2 obtained from the HAL469 interaction. The ‘ T B ​ G ​ E ​ ( ω ) superscript 𝑇 𝐵 𝐺 𝐸 𝜔 T^{BGE}(\omega) +ADC(3)’ results of the present work are compared to BHF and to the exact solution. The last line is the breakup energy for splitting the system in 4 He clusters (of total energy A / 4 × A/4\times 5.09 MeV).

$E_0^A$ [MeV]		4He		16O		40Ca
BHF Inoue et al. (2015)		-8.2		-34.7		-112.7
$T^{BGE}\left(\omega \right)$+ADC(3)		-4.80(0.03)		-17.9 (0.3) (1.8)		-75.4 (6.7) (7.5)
Exact calc. Nemura (2014)		-5.09		–		–
Experiment		-28.3		-127.7		-342.0
Separation into 4He clusters:				-2.46 (0.3) (1.8)		24.5 (6.7) (7.5)

Table 2. Table 2: Computed matter and charge radii of 16 O and 40 Ca using M PS =469 MeV for N max subscript 𝑁 max N_{\mbox{\footnotesize{max}}} =11 and ℏ ​ Ω Planck-constant-over-2-pi Ω \hbar\Omega{} =11 MeV. Results are given for different levels of approximations and the charge radii from the full T B ​ G ​ E ​ ( ω ) superscript 𝑇 𝐵 𝐺 𝐸 𝜔 T^{BGE}(\omega) plus ADC(3) are compared to the experimental values. For charge radii, we assumed the physical charge distributions of the proton and the neutron (see Ref.  Cipollone  et al. ( 2015 ) for details).

					4He		16O		40Ca
	$r_{pt-matter}$[fm]:		BHF Inoue et al. (2015)		2.09		2.35		2.78
			HF		1.62		2.39		2.78
			$T^{BGE}\left(\omega \right)$ + ADC(3)		1.67		2.64		2.97
	$r_{charge}$[fm]:		$T^{BGE}\left(\omega \right)$ + ADC(3)		1.89		2.79		3.10
			Experiment Vries et al. (1987); Angeli and Marinova (2013)		1.67		2.73		3.48

Equations

$$g_{\alpha \beta }(\omega )=$$

$${\Sigma }_{\alpha \beta }^{\star }(\omega )=$$

$${\Sigma }_{\alpha \beta }^{(\infty )}(\omega )=$$

$$=$$

$$\sum _{\beta }\left\{⟨\alpha |\frac{p^2}{2m_N}|\beta ⟩+{\Sigma }_{\alpha \beta }^{(\infty )}(\omega ={\epsilon }_r^{HF})\right\}{\psi }_{\beta }^r=$$

$$V_{rs,pq}=$$

$$L_2=$$

$$E_0^A[N_{\text{mbox}max},\hslash \Omega ]=$$

Full text

Introduction. Quantum chromodynamics (QCD) is expected to ultimately explain the structure and the interactions of all hadronic systems, together with small corrections of electroweak origin. At the moment, systematic and non-perturbative calculations of QCD can be carried out only by Lattice QCD (LQCD). Indeed, high-precision studies have been shown to be possible e.g., in the case of single-hadron masses Borsanyi et al. (2015). “Direct” calculations of multi-baryon systems have also been attempted on a lattice by several groups Beane et al. (2012); Yamazaki et al. (2012, 2015); Orginos et al. (2015); Berkowitz et al. (2015). However, the typical excitation energy $\Delta E$ for multi-baryons is one to two orders of magnitude smaller than ${\cal O}(\Lambda_{\rm QCD})$. Accordingly, high statistic data with very large Euclidean times t\mathrel{\mathchoice{\lower 2.58334pt\vbox{\halign{\m@th\displaystyle\hfill#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.58334pt\vbox{\halign{\m@th\textstyle\hfill#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.58334pt\vbox{\halign{\m@th\scriptstyle\hfill#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.58334pt\vbox{\halign{\m@th\scriptscriptstyle\hfill#\hfil\cr>\crcr\sim\crcr}}}}\hbar/\Delta E\sim 10-100 fm/c is required. This is still far beyond reach due to exponentially increasing errors in $t$ and $A$ (the atomic number), as demonstrated theoretically and numerically in recent studies  Iritani et al. (2016); Aoki et al. (2016); Iritani (2016).

In this work, we follow a different route and perform a first ab initio study for medium mass atomic nuclei directly based on QCD by taking a two-step strategy. In the first step, we extract the nuclear force from ab initio LQCD calculations with the HAL QCD method. This procedure generates consistent two-, three- and many-nucleon forces in a systematic way Ishii et al. (2007); Aoki et al. (2010, 2012, 2013). Note that the HAL QCD interactions obtained directly from spatiotemporal hadron-hadron correlations on the lattice are faithful to the scattering phase shifts and binding-energies by construction. This is done systematically without the fitting procedures needed by phenomenological potentials. Furthermore, explorative studies of three-nucleon potential show that this is weaker than the corresponding nucleon-nucleon (NN) force, in agreement with the empirically observed hierarchy of nuclear forces Doi et al. (2012). Thus, the interaction can be applied to larger nuclei. In the second step, we calculate the properties of nuclei with ab initio many-body methods using the LQCD potentials as input. The advantage of the HAL QCD approach is that one can extract the potential dictating all the elastic scattering states below the inelastic threshold from the lattice data for t\mathrel{\mathchoice{\lower 2.58334pt\vbox{\halign{\m@th\displaystyle\hfill#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.58334pt\vbox{\halign{\m@th\textstyle\hfill#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.58334pt\vbox{\halign{\m@th\scriptstyle\hfill#\hfil\cr>\crcr\sim\crcr}}}{\lower 2.58334pt\vbox{\halign{\m@th\scriptscriptstyle\hfill#\hfil\cr>\crcr\sim\crcr}}}}1 fm/c Ishii et al. (2012). This makes the LQCD calculation of potentials affordable with reasonable statistics, together with the help of advanced computational algorithms Doi and Endres (2013); Detmold and Orginos (2013); Günther et al. (2013); Nemura (2016). We can then take the advantage of the recent developments in nuclear many-body theories to calculate various information on nuclei such as binding energies and spectral distributions. Note that a similar two-step strategy has also been taken in Ref. Barnea et al. (2015); Contessi et al. (2017) where, however, effective field theories have been used to fit the LQCD data.

Past LQCD studies in the flavor SU(3) limit by the HAL QCD collaboration have led to interactions in both the nucleon and hyperon sectors with masses of the pseudo-scalar meson (which corresponds to the pion) as low as M_{PS}$$={}469 MeV/c2. In these cases, potentials have been obtained for the 1S0 and the coupled 3S1-3D1 channels Inoue et al. (2011, 2012). Exploratory calculations based on these HAL QCD potentials were performed in the Brueckner Hartree-Fock (BHF) approach Inoue et al. (2013, 2015). This is quantitative enough to give the essential underlying physics for infinite matter but it is less reliable in finite systems. More sophisticated calculations are needed in order to go beyond the mean-field level, which is necessary to properly predict binding energies and to describe the truly complex structures of nuclei at low energy. BHF becomes even more questionable for finite nuclei due to assumptions with the unperturbed single particle spectrum where there is a problem in the choice between a continuous or a gap form, neither of which is completely satisfactory.

Ab initio theories for medium mass nuclei have advanced greatly in the past decade and methods such as coupled cluster and self-consistent Green’s function (SCGF) are now routinely employed to study the structure of full isotopic chains up to Ca and Ni Somà et al. (2014a); Hergert et al. (2014); Hagen et al. (2014); Dickhoff and Barbieri (2004), with the inclusion of three-nucleon forces Carbone et al. (2013). Their use with soft interactions from chiral effective field theory have led to first principle predictions of experimental total binding energies Hebeler et al. (2015); Cipollone et al. (2013) and nuclear radii Ekström et al. (2015); Lapoux et al. (2016); Garcia Ruiz et al. (2016) with unprecedented accuracies. Interactions like the HAL QCD potentials pose a bigger challenge for calculations of medium mass isotopes due to their short-distance inter-nucleon repulsion. However, the SCGF method has the advantage that two-nucleon scattering (ladder) diagrams at large momenta—needed to resolve the short-range repulsion—can be dealt with explicitly by solving the Bethe-Goldstone equation (BGE) in the excluded space Barbieri and Dickhoff (2002); Barbieri (2006); Barbieri and Hjorth-Jensen (2009). This route was exploited in the past to study spectral strength distribution but we extend it to binding energies in the present work and find that this is accurate enough to make statements on the performance of the present HAL QCD potentials. Thus, this work is also a first step toward advancing many-body approaches that can handle hard interactions for large atomic masses. This is also important since new LQCD calculations at nearly the physical pion mass are currently in progress  Doi et al. (2016); Ishii et al. (2016); Sasaki et al. (2016); Nemura et al. (2016).

The BHF study of Ref. Inoue et al. (2013) showed that the HAL QCD interactions in the SU(3) limit do not bind at very large pion masses except for the lowest available value M_{PS}$$={}469 MeV/c2. The saturation of nuclear matter in this case is also confirmed by later SCGF calculations Carbone (2015). Thus, this is a suitable choice to investigate possible self-bound nuclei at large pion masses. In this work, we will focus on this potential and refer to it as the ‘HAL469${}_{\hbox{\small\em SU(3)}}$’ interaction, or ‘HAL469’ for simplicity.

Formalism. We follow closely the approach of Ref. Barbieri and Hjorth-Jensen (2009) and focus on the single-particle propagator given by Fetter and Walecka (2003); Dickhoff and Van Neck (2008):

[TABLE]

where $\mathcal{X}_{\alpha}^{n}\equiv\langle\Psi_{n}^{A+1}|c_{\alpha}^{\dagger}|\Psi_{0}^{A}\rangle$ ($\mathcal{Y}_{\alpha}^{k}\equiv\langle\Psi_{k}^{A-1}|c_{\alpha}|\Psi_{0}^{A}\rangle$) are the spectroscopic amplitudes, $\varepsilon_{n}^{+}\equiv E_{n}^{A+1}-E_{0}^{A}$ ($\varepsilon_{k}^{-}\equiv E_{0}^{A}-E_{k}^{A-1})$ are their quasiparticle energies and $c_{\alpha}^{\dagger}(c_{\alpha})$ are the second quantisation creation (annihilation) operators. In these definitions, $|\Psi_{n}^{A+1}\rangle$ and $|\Psi_{k}^{A-1}\rangle$ represent the exact eigenstates of the $(A\pm 1)$-nucleon system, while $E_{n}^{A+1}$ and $E_{k}^{A-1}$ are the corresponding energies. We perform calculations within a spherical harmonic oscillator model space, indicated as $P$, and use Greek indices, $\alpha$, $\beta$, … to label its basis functions. Within this space, the one-body propagator is obtained by solving the Dyson equation with an irreducible self-energy given by

[TABLE]

This expression is the sum of a mean-field (MF) term, $\Sigma^{(\infty)}$, and the contributions from dynamical correlations. The coupling (D) and interactions matrices (K and C) are computed in the third-order algebraic diagrammatic construction [ADC(3)] approximation that generates non-perturbative all-order summations of 2p1h and 2h1p configurations, as shown in Fig. 1 Schirmer et al. (1983); Barbieri et al. (2007); Barbieri and Carbone (2017). We follow the sc0 approximation of Refs. Barbieri and Hjorth-Jensen (2009); Somà et al. (2014b), in which $\Sigma^{(\infty)}$ is calculated exactly from the fully-dressed propagator $g(\omega)$ while D, K and C are written in terms of a simplified MF reference propagator, $g^{REF}(\omega)$, that is chosen to best approximate $g(\omega)$ through its first two moments at the Fermi energy (see Ref. Barbieri and Hjorth-Jensen (2009) for details).

For forces with a sizeable short-range repulsion, like the HAL QCD interactions, usual truncations of the oscillator space (of up to 12 shells in this case) are not sufficient and a resummation of ladder diagrams in the excluded Hilbert space, $Q\equiv{\bf 1}-P$, is required. We do this by solving the BGE in $Q$ according to Refs. Hjorth-Jensen et al. (1995); Engeland et al. and add the corresponding diagrams to the MF self-energy, which becomes energy dependent Barbieri and Hjorth-Jensen (2009):

[TABLE]

where $T^{BGE}_{\alpha\gamma,\,\beta\delta}(\omega)$ are the elements of the scattering t-matrix in the excluded space. We then extract a static effective interaction that we use to calculate the ADC(3) self-energy (the last term of Eq. (2)) within the model space $P$. To do this, we solve the Hartree-Fock (HF) equations with the MF potential of Eq. (3):

[TABLE]

where latin indices label HF states, and define a static interaction in this HF basis similarly to Refs. Gad and Müther (2002); Barbieri and Hjorth-Jensen (2009):

[TABLE]

The $V_{rs,pq}$ matrix elements are then transformed back to the harmonic oscillator space to be used in the computations.

Note that the BGE used to generate $T^{BGE}(\omega)$ resums scattering states where at least one nucleon is outside the whole model space. Hence, it does not depend on the isotope being computed, neither it suffers from ambiguities with the choice of the single-particle spectrum at the Fermi surface encountered with the usual G-matrix used in BHF calculations. Instead, our approach is very similar to the Bloch-Horowitz method of Ref. Luu et al. (2004), of which our $T^{BGE}(\omega)$ is also a first approximation. The ADC(3) computation accounts for all types of correlations inside the model space, including all remaining ladder diagrams. This ensures a complete many-body calculation that accounts for short distance repulsion.

It is well known that short-range repulsion, that is accounted for by Eq. (3), has the double effect of reducing the spectral strength for dominant quasiparticle peaks and of relocating it to large momenta and large quasiparticle energies Dickhoff and Barbieri (2004). Since we cannot currently calculate the location of strength at high momenta, it is not possible to quantify the magnitude of these two effects. However, they contribute to the Koltun sum rule for total binding energy with opposite signs and must cancel to a large extent. Thus, we chose to neglect both contributions and maintain a static $\Sigma^{(\infty)}$ to solve the Dyson equation. This is currently the major approximation in our calculations and its uncertainty is best estimated from the benchmark on 4He below. Resolving this requires a proper extension of the present SCGF formalism and it will be the subject of future work.

Infrared convergence. The one-body propagators of 4He, 16O and 40Ca are calculated in spherical harmonic oscillator spaces of different frequencies, $\hbar\Omega$, and increasing sizes up to $N_{\mbox{\footnotesize{max}}}$=\max\{2n+\ell\}$$=11 (and $N_{\mbox{\footnotesize{max}}}\leq$ 9 for 40Ca). The scattering matrix $T^{BGE}(\omega)$ is calculated for each frequency and model space and then used to derive the interactions of Eq. (5). We subtract the kinetic energy of the center of mass according to Ref. Hergert and Roth (2009) and calculate the intrinsic ground state energy from $g(\omega)$ using the Koltun sum rule. The same lattice simulation setup used to generate the HAL469 interaction gives a nucleon mass of $m_{N}$=1161.1 MeV/c2 in addition to the pseudo-scalar mass of $M_{PS}$=469 MeV/c2. Thus, we employ this value of $m_{N}$ in all the kinetic energy terms. Fig. 2 displays the ground state energies obtained with our $T^{BGE}(\omega)$ plus ADC(3) method. As expected, the complete resummation of ladder diagrams outside the model space tames ultraviolet corrections and results in a rather flat behaviour of the total energies for $\hbar\Omega{}\approx{}$5-20 MeV. Still, there remain some hints of oscillations with respect to $\hbar\Omega$ that could be linked to the HO truncation and to the neglect of spectral strength at high momenta as explained above.

The effectiveness of the ladder resummation is better recognised by considering the infrared (IR) convergence of the total binding energies, from where one can also extrapolate to a complete set. Ref. More et al. (2013) established that a harmonic oscillator model space, of frequency $\hbar\Omega$ and truncated to the first $N_{\mbox{\footnotesize{max}}}+1$ shells, behaves as a hard wall spherical box of radius

[TABLE]

where $b\equiv\sqrt{\hbar c^{2}/m_{N}\Omega}$ is the oscillator length. Given a soft interaction that is independent of the model space, if the frequency is large enough (i.e., $b$ is small) then ultraviolet (UV) degrees of freedom are converged. In this case, the calculated ground state energies are expected to converge exponentially when increasing the effective radius $L_{2}$:

[TABLE]

For the bare HAL469 interaction, if we use the SCGF without ladder sums outside the model space, the extrapolation according to Eq. (7) will fail because the short-distance repulsion requires extremely large $N_{\mbox{\footnotesize{max}}}$ ($\gg$ 20) to reach UV convergence, while our many-body space $P$ is limited to $N_{\mbox{\footnotesize{max}}}$=11. This is shown by the left panel of Fig. 3. However, our complete calculations use the BGE to resum all missing ladder diagrams within $Q$. Adding these to the two-particle ladders that are generated by the ADC(3) leads to a complete diagonalization of short-distance degrees of freedom, independently of the choice of $P$. The resulting dependence of the ground state energy on $L_{2}$ is shown by the middle panel and it now follows the behaviour dictated by Eq. (7). Note that single particle energies are still needed to calculate $T^{BGE}(\omega)$ but these can be identified with the free particle spectrum in the space $Q$, at large momenta. Accordingly, the flat region in Fig. 3 becomes broader as we increase $N_{\mbox{\footnotesize{max}}}$ (that is, when the boundary between the $P$ and $Q$ spaces moves away form the Fermi surface).

The binding energy of 4He for HAL469 was found to be -5.09 MeV with the exact Stochastic Variational calculations Nemura (2014), which we will use to benchmark our approach. The solid lines in Fig. 3 are the result of a nonlinear least-squares fit to Eq. (7). The points diverging from the exponential behaviour at large $L_{2}$ are assumed not to be UV converged due to the above-mentioned approximations and are excluded from the fit but are still shown in the figure. From calculations up to $\hbar\Omega$=50 MeV and the IR extrapolation, we estimate a converged binding energy of 4.80(3) MeV for 4He, where the error corresponds to the uncertainties in the extrapolation. The calculations for the other isotopes converge similarly to 4He, and we show the IR extrapolation of 40Ca in the right panel of Fig. 3, for completeness. The figure is also indicative that this methodology can be successfully applied to heavier nuclei. These results clearly show that SCGF can handle relatively hard potentials by resumming of ladders through the Bethe-Goldstone equation and they give confidence that the short-range repulsion of HAL469 is accounted for accurately.

Results. The calculated ground state energies of 4He, 16O and 40Ca are summarised in Tab. 1, together with BHF results obtained with the same gap choice and methods of Ref. Inoue et al. (2015). For 4He, the complete $T^{BGE}(\omega)$+ADC(3) result deviates from the exact solution for 4He by less than 10%. Since the SCGF approach resums linked diagrams, and thus it is size extensive, one should expect that similar errors will apply for larger isotopes Shavitt and Bartlett (2009). Thus, Tab. 1 shows both the uncertainties in the IR extrapolation More et al. (2013) and an error for the many-body truncations, for we which we take a conservative estimate of 10% based on our finding for 4He. The SCGF results are sensibly less bound than our previous BHF results Inoue et al. (2015) and we interpret this as a limitation of BHF theory.

A key feature of our calculations is the use of an harmonic oscillator space, which effectively confines all nucleons. The last line of Tab. 1 reports the deduced breakup energies for separating the computed ground states into infinitely distant 4He clusters. The 16O is unstable with respect to 4-$\alpha$ break up, by $\approx{}$2.5 MeV. Allowing an error in our binding energies of more than 10% could make oxygen bound but only very weakly. This is in contrast to the experimental results, at the physical quarks masses, where the 4-$\alpha$ breakup requires 14.4 MeV. On the other hand, 40Ca is stable with respect to breakup in $\alpha$ particles by $\approx{}$24 MeV. We expect that these observations are rather robust even when we consider the (LQCD) statistical errors in the HAL469 interaction. While such statistical fluctuations introduce additional $\sim$10% errors on binding energies Inoue et al. (2015), they are expected to be strongly correlated among 4He, 16O and 40Ca. Hence, for QCD in the SU(3) limit at $M_{PS}$=469 MeV/c2, we find that the deuteron is unbound Inoue et al. (2012) and 16O is only just slightly above the threshold for $\alpha$ breakup, while 4He and 40Ca are instead bound. The HAL469 interaction has the lowest $M_{PS}$ value among those considered in Refs. Inoue et al. (2011, 2012), while from Ref. Inoue et al. (2013) we know that it is the only one saturating nuclear matter (although not at the physical saturation point). Moreover, we have tested that SCGF attempts at calculating asymmetric isotopes, like 28O, predict strongly unbound systems even for HAL469. All these results together suggest that, when lowering of the pion mass toward its physical value, closed shell isotopes are created at first around the traditional magic numbers. This hypothesis should also be seen in the light of the limitations in the present HAL469 Hamiltonian, which was built to include only the 1S0, 3S1 and 3D1 partial waves and therefore neglects the three-body forces and spin-orbit interactions. The missing $P$ waves and Coulomb force are repulsive but could be compensated by an attracting three-body force.

Figure 4 demonstrates the spectral strength distribution of 16O obtained for $N_{\mbox{\footnotesize{max}}}$=11 and $\hbar\Omega{}$=11 MeV. Quasiparticle fragments corresponding to spin-orbit partners do not split due to the absence of a spin-orbit term in HAL469. Otherwise, all the remaining qualitative features of the experimental spectral distribution are seen also for the $M_{PS}$=469 MeV/c2. The root mean square radii are given Tab. 2 for the same model space and oscillator frequency. Although the total binding energies are 15-20% of the experimental value (Tab. 1), the computed charge radii are about the same as the experiment. This is due to the fact that the heavy nucleon mass ($m_{N}$=1161.1 MeV/${\rm c}^{2}$) used here reduces the motion of the nucleons inside the nuclei and counterbalances the effect of weak attraction of the HAL469 potential. We have also checked the dependence of the computed radii on the effective model space size, Eq. (6), and found a rather flat converged region for all three isotopes, although the values still oscillate by about 0.05 fm with changing $L_{2}$. The values in Tab. 2 are all calculated in the middle of this plateau. The HF approach of Eq. (4) and the standard BHF give similar radii in spite that they predict very different binding energies. The final radii are then increased by many-body correlations and, for all nuclei, the full $T^{BGE}(\omega)$ plus ADC(3) calculations pushes the matter distribution to larger discances. However, we note that accounting for the neglected high-momentum components—as discussed below Eq. (5)—tends to enhance the central density and would slightly reduce the calculated radii Müther and Sick (2004).

Summary. We investigated the use of the BGE in an ab initio approach and used it to resum missing two-nucleon scattering diagrams outside the usual truncations of the many-body space, while the full ADC(3) method has been retained within the model space itself. A benchmark on 4He shows that the present implementation works relatively well and it allows to solve the self-consistent Green’s function for the HAL QCD potentials derived from Lattice QCD. An investigation of the IR convergence of the ground state energies, following the work of Ref. More et al. (2013), indicates that SCGF can handle relatively hard potentials such as the HAL469, even masses as large as A=40. This opens a new path that allows full ab initio calculations of large nuclei even with hard nuclear interactions.

The present accuracy is sufficient to make quantitative statements on doubly magic nuclei, which are less bound compared to earlier BHF estimates for the HAL QCD potentials. Here, we have found that the behaviour when lowering the pion mass towards its physical value is consistent with the idea that nuclei near to the traditional magic numbers are formed at first. At $M_{PS}$=469 MeV/c2, in the SU(3) limit of QCD, both 4He and 40Ca have bound ground states while the deuteron is unbound and 16O is likely to decay into four separate alpha particles. However, 16O is already close to become bound. This suggests that the region of $M_{PS}\sim{}$500 MeV/c2 marks a transition between an unbound nuclear chart and the emergence of bound isotopes. Further studies at lower pion masses will be pivotal to test these findings and should be possible in the near future since LQCD simulations for nuclear and hyperon forces down to physical quark masses are currently underway Doi et al. (2016); Ishii et al. (2016); Sasaki et al. (2016).

Important future work will also be the inclusion the spin-orbit as well as three-nucleon forces. Proof of principle LQCD calculations for these interactions are available Murano et al. (2014); Doi et al. (2012) and follow the hierarchy of nuclear forces with three-nucleon terms smaller than the NN contributions.

Acknowledgements. We thank the HAL QCD Collaboration for providing the HAL469${}_{\hbox{\small\em SU(3)}}$ interaction. Computations of the $T^{BGE}(\omega)$ were performed using the CENS codes from Ref. Engeland et al. . This work was supported in part by the United Kingdom Science and Technology Facilities Council (STFC) under Grants No. ST/L005743/1 and ST/L005816/1, by the Royal Society International Exchanges Grant No. IE150305, by Japanese Grant-in-Aid for Scientific Research (JP15K17667and (C)26400281), by MEXT Strategic Program for Innovative Research (SPIRE) Field 5, by a priority issue (Elucidation of the fundamental laws and evolution of the universe) to be tackled by using Post “K” Computer, and by the Joint Institute for Computational Fundamental Science (JICFuS). T.D and T.H. were supported partially by the RIKEN iTHEMS Program and iTHES Project. Calculations were performed at the DiRAC Complexity system at the University of Leicester (BIS National E- infrastructure capital Grant No. ST/K000373/1 and STFC Grant No. ST/K0003259/1).