Novel chiral Hamiltonian and observables in light and medium-mass nuclei

TL;DR

This paper introduces a new chiral effective field theory-based Hamiltonian for light and medium-mass nuclei, showing promising results in reproducing energies and spectra, and offering a competitive alternative to existing models.

Contribution

A novel Hamiltonian combining next-to-next-to-leading order three-nucleon operators with an existing two-body interaction, tested across various nuclei with encouraging results.

Findings

Accurately reproduces binding energies and spectra in light nuclei.

Performs well for medium-mass nuclei, from oxygen to nickel.

Outperforms some existing Hamiltonians in spectral quality.

Abstract

A novel parameterisation of a Hamiltonian based on chiral effective field theory is introduced. Specifically, three-nucleon operators at next-to-next-to-leading order are combined with an existing (and successful) two-body interaction containing terms up to next-to-next-to-next-to-leading order. The resulting potential is labelled $N\!N\!$+$3N\text{(lnl)}$. The objective of the present work is to investigate the performance of this new Hamiltonian across light and medium-mass nuclei. Binding energies, nuclear radii and excitation spectra are computed using no-core shell model and self-consistent Green's function approaches. Calculations with $N\!N\!$+$3N\text{(lnl)}$ are compared to two other representative Hamiltonians currently in use, namely NNLO$_{\text{sat}}$ and the older $N\!N\!$+$3N(400)$. Overall, the performance of the novel interaction is very encouraging. In light nuclei, total energies are generally in good agreement with experimental data. Known spectra are also well reproduced with a few notable exceptions. The good description of ground-state energies carries on to heavier nuclei, all the way from oxygen to nickel isotopes. Except for those involving excitation processes across the $N=20$ gap, which is overestimated by the new interaction, spectra are of very good quality, in general superior to those obtained with NNLO$_{\text{sat}}$. Although largely improving on $N\!N\!$+$3N(400)$ results, charge radii calculated with $N\!N\!$+$3N\text{(lnl)}$ still underestimate experimental values, as opposed to the ones computed with NNLO$_{\text{sat}}$ that successfully reproduce available data on nickel. On the whole, the new two- plus three-nucleon Hamiltonian introduced in the present work represents a promising alternative to existing nuclear interactions.

Tables

Table 1. Table 1: 4 He ground-state energies and mean values of the five chiral 3 ​ N 3 𝑁 3N N 2 LO terms (in MeV) for the N ​ N 𝑁 𝑁 N\!N\! + 3 ​ N ​ ( 400 ) 3 𝑁 400 3N(400)   Roth  et al. ( 2012 ) , NNLO sat sat {}_{\text{sat}}   Ekström  et al. ( 2015 ) , N ​ N 𝑁 𝑁 N\!N\! + 3 ​ N ​ (lnl) 3 𝑁 (lnl) 3N\text{(lnl)} (present work), and 2 ​ N 2 𝑁 2N + 3 ​ N ​ ( 500 ) c D ​ 0.83 3 𝑁 subscript 500 subscript 𝑐 𝐷 0.83 3N(500)_{c_{D}\!0.83} Gazit  et al. ( 2019 ) Hamiltonians. All interactions are bare, i.e. not evolved via SRG techniques. The experimental 4 He ground-state energy is -28.29 MeV.

4He:	$⟨H⟩$	$⟨c_1⟩$	$⟨c_3⟩$	$⟨c_4⟩$	$⟨c_D⟩$	$⟨c_E⟩$
$NN$+$3N\left(400\right)$	-28.28	-0.06	1.27	-3.93	-0.28	-0.66
NNLO${}_{\text{sat}}$	-28.43	-0.24	-0.73	-3.76	1.39	0.42
$NN$+$3N\text{(lnl)}$	-28.25	-0.18	-1.36	-3.27	0.74	0.43
$2N$+$3N{\left(500\right)}_{c_{D}0.83}$	-28.36	-0.26	-1.50	-3.79	0.78	0.30

Table 2. Table 2: Values of the five LECs ( c 1 subscript 𝑐 1 c_{1} , c 3 subscript 𝑐 3 c_{3} , c 4 subscript 𝑐 4 c_{4} in GeV -1 ) that define the 3 ​ N 3 𝑁 3N force at N 2 LO for the N ​ N 𝑁 𝑁 N\!N\! + 3 ​ N ​ ( 400 ) 3 𝑁 400 3N(400)   Roth  et al. ( 2012 ) , NNLO sat sat {}_{\text{sat}}   Ekström  et al. ( 2015 ) , N ​ N 𝑁 𝑁 N\!N\! + 3 ​ N ​ (lnl) 3 𝑁 (lnl) 3N\text{(lnl)} (present work), and 2 ​ N 2 𝑁 2N + 3 ​ N ​ ( 500 ) c D ​ 0.83 3 𝑁 subscript 500 subscript 𝑐 𝐷 0.83 3N(500)_{c_{D}\!0.83}   Gazit  et al. ( 2019 ) Hamiltonians. The c 1 subscript 𝑐 1 c_{1} , c 3 subscript 𝑐 3 c_{3} and c 4 subscript 𝑐 4 c_{4} constants enter both the 2 ​ N 2 𝑁 2N and the 3 ​ N 3 𝑁 3N sectors of the Hamiltonian, they are equal for N ​ N 𝑁 𝑁 N\!N\! + 3 ​ N ​ ( 400 ) 3 𝑁 400 3N(400) , N ​ N 𝑁 𝑁 N\!N\! + 3 ​ N ​ (lnl) 3 𝑁 (lnl) 3N\text{(lnl)} and 2 ​ N 2 𝑁 2N + 3 ​ N ​ ( 500 ) c D ​ 0.83 3 𝑁 subscript 500 subscript 𝑐 𝐷 0.83 3N(500)_{c_{D}\!0.83}  since these share the same bare 2 ​ N 2 𝑁 2N interaction from Ref.  Entem and Machleidt ( 2003 ) .

	$c_1$	$c_3$	$c_4$	$c_D$	$c_E$
$NN$+$3N\left(400\right)$	-0.81	-3.20	5.40	-0.20	0.098
NNLO${}_{\text{sat}}$	-1.122	-3.925	3.766	0.817	-0.040
$NN$+$3N\text{(lnl)}$	-0.81	-3.20	5.40	0.70	-0.06
$2N$+$3N{\left(500\right)}_{c_{D}0.83}$	-0.81	-3.20	5.40	0.83	-0.052

Table 3. Table 3: Comparison between the ground state energies of 16 O (in MeV) as computed with the ADC(2) and ADC(3) many-body approximations and with the NCSM. The SCGF results are obtained in a full e max = 13 subscript 𝑒 max 13 e_{\text{max}}=13 space and one could expect residual errors of a few % due to the many-body truncations beyond ADC(3), as detailed in the text. The NCSM results for NNLO sat sat {}_{\text{sat}} and N ​ N 𝑁 𝑁 N\!N\! + 3 ​ N ​ (lnl) 3 𝑁 (lnl) 3N\text{(lnl)} show uncertainties arising from N max subscript 𝑁 max N_{\rm max} extrapolation and importance-truncation. The NCSM result for N ​ N 𝑁 𝑁 N\!N\! + 3 ​ N ​ ( 400 ) 3 𝑁 400 3N(400) is taken from Ref.  Hergert  et al. ( 2013 ) .

$E_{g.s.}$ (16O):	ADC(2)	ADC(3)	NCSM
$NN$+$3N\left(400\right)$	-128.54	-130.81	-130(2)
NNLO${}_{\text{sat}}$	-124.63	-126.23	-125(5)
$NN$+$3N\text{(lnl)}$	-123.91	-127.27	-126(2)
Experiment:			-127.62

Full text

Recent advances in nuclear structure theory have led to the availability of several complementary ab initio many-body techniques applicable to light and medium-mass nuclei as well as nuclear matter. After successful benchmarks of different approaches, the focus is moving to the development of improved models of nuclear Hamiltonians, currently representing the largest source of uncertainty in ab initio calculations of nuclear systems. In particular, none of the existing two- plus three-body interactions is capable of satisfactorily reproducing all the observables of interest in medium-mass nuclei.

Purpose

A novel parameterisation of a Hamiltonian based on chiral effective field theory is introduced. Specifically, three-nucleon operators at next-to-next-to-leading order are combined with an existing (and successful) two-body interaction containing terms up to next-to-next-to-next-to-leading order. The resulting potential is labelled $N\!N\!$+$3N\text{(lnl)}$. The objective of the present work is to investigate the performance of this new Hamiltonian across light and medium-mass nuclei.

Methods

Binding energies, nuclear radii and excitation spectra are computed using state-of-the-art no-core shell model and self-consistent Green’s function approaches. Calculations with $N\!N\!$+$3N\text{(lnl)}$ are compared to two other representative Hamiltonians currently in use, namely NNLO${}_{\text{sat}}$ and the older $N\!N\!$+$3N(400)$.

Results

Overall, the performance of the novel $N\!N\!$+$3N\text{(lnl)}$ interaction is very encouraging. In light nuclei, total energies are generally in good agreement with experimental data. Known spectra are also well reproduced with a few notable exceptions. The good description of ground-state energies carries on to heavier nuclei, all the way from oxygen to nickel isotopes. Except for those involving excitation processes across the $N=20$ gap, which is overestimated by the new interaction, spectra are of very good quality, in general superior to those obtained with NNLO${}_{\text{sat}}$. Although largely improving on $N\!N\!$+$3N(400)$ results, charge radii calculated with $N\!N\!$+$3N\text{(lnl)}$ still underestimate experimental values, as opposed to the ones computed with NNLO${}_{\text{sat}}$ that successfully reproduce available data on nickel.

Conclusions

The new two- plus three-nucleon Hamiltonian introduced in the present work represents a promising alternative to existing nuclear interactions. In particular, it has the favourable features of (i) being adjusted solely on $A=2,3,4$ systems, thus complying with the ab initio strategy, (ii) yielding an excellent reproduction of experimental energies all the way from light to medium-heavy nuclei and (iii) behaving well under similarity renormalisation group transformations, with negligible four-nucleon forces being induced, thus allowing large-scale calculations up to medium-heavy systems. The problem of the underestimation of nuclear radii persists and will necessitate novel developments.

I Introduction

In the last decade, advances in many-body approaches and inter-nucleon interactions have enabled significant progress in ab initio calculations of nuclear systems. At present, several complementary methods to solve the (time-independent) many-body Schrödinger equation are available, tailored to either light systems Navrátil et al. (2016); Carlson et al. (2015), medium-mass nuclei Dickhoff and Barbieri (2004); Lähde et al. (2014); Hagen et al. (2014a); Hergert et al. (2016); Stroberg et al. (2017); Tichai et al. (2018) or extended nuclear matter Carbone et al. (2014); Lynn et al. (2016); Drischler et al. (2019). New developments, which promise to extend (most of) these methods to higher accuracy and/or heavy nuclei, are being currently proposed Tichai et al. (2019a, b).

Over the past few years, benchmarks calculations have allowed assessment of the systematic errors associated with both the use of a necessarily finite-dimensional Hilbert space and the truncation of the many-body expansion at play in each of the formalisms of interest. In state-of-the-art implementations, these errors add up to at most $5\%$, much less than the uncertainty attributable to the input nuclear Hamiltonian Hergert et al. (2014); Hebeler et al. (2015); Lapoux et al. (2016); Duguet et al. (2017); Leistenschneider et al. (2018). As a result, ab initio calculations have also acquired the role of diagnostic tools as the focus of the community is shifting towards developing improved models of nuclear interactions.

The large majority of these developments currently takes place in the context of chiral effective field theory ($\chi$EFT) Epelbaum et al. (2009); Machleidt and Entem (2011) based on Weinberg’s power counting (WPC) Weinberg (1990, 1991); Ordóñez and van Kolck (1992). Building on a low-energy expansion with nucleons and pions as explicit degrees of freedom, $\chi$EFT provides a framework in which two- and many-nucleon interactions can be systematically derived, in principle with uncertainties associated with each order of the expansion. This approach to $\chi$EFT suffers from issues related to its non renormalizability and different strategies have been proposed to resolve this problem, requiring the use of a different power counting and a mixture of non-perturbative and perturbative calculations Hammer et al. (2019). However, applications of these procedures are still in their infancy stage. From the practical point of view, the standard implementations of $\chi$EFT interactions are much more advanced and therefore in better position to provide useful predictions. The present study follows the conventional approach of seeking fully non-perturbative solutions of the many-body Schrödinger equation and focuses on the performance of ab initio calculations for isotopes up to medium masses.

Within this framework, one of the first successful nuclear Hamiltonians combined a next-to-next-to-next-to-leading order (N3LO) two-nucleon ($2N$) force Entem and Machleidt (2003) with a local N2LO three-nucleon ($3N$) interaction Navrátil (2007), whose associated cutoff was subsequently reduced to 400 MeV/c to optimise its behaviour under similarity renormalization group (SRG) transformations Roth et al. (2012). This Hamiltonian, labelled $N\!N\!$+$3N(400)$, has constituted a standard for many early applications of ab initio techniques in light as well as medium-mass nuclei. As systematic calculations beyond the light sector became available, deficiencies associated to this interaction emerged. In particular, it was shown to lead to an underestimation of nuclear radii and a substantial overbinding, i.e. an overestimation of total binding energies in medium- and heavy-mass nuclei Somà et al. (2014a); Binder et al. (2014); Lapoux et al. (2016). Recently, this behaviour was related to the poor description of saturation properties of symmetric nuclear matter Simonis et al. (2017); Drischler et al. (2019).

While the initial success of $N\!N\!$+$3N(400)$ represented an important breakthrough, its poor performance for isotopes above the oxygens prompted new work aimed at improving predictions of nuclear saturation. A successful route was put forward by Ekström and collaborators in Ref. Ekström et al. (2015). There, in contrast to the standard strategy of constraining the Hamiltonian only on few-nucleon data, the binding energies and charge radii of nuclei up to $A=25$ were employed in the simultaneous fit of the low-energy constants (LECs) that enter the $2N$ and $3N$ interactions. The resulting next-to-next-to-leading order Hamiltonian, labelled NNLO${}_{\text{sat}}$, provides a much improved description of charge radii and fair saturation properties of extended matter. However, this is obtained at the price of deteriorating certain properties of two- and few-body systems Ekström et al. (2015); Lapoux et al. (2016). Furthermore, mild underbinding is observed in medium-mass isotopes. Another family of Hamiltonians that have proven to perform well for either binding energies or radii Simonis et al. (2017) was introduced by Hebeler and collaborators in Ref. Hebeler et al. (2011). One drawback of these interactions is that $2N$ and $3N$ sectors are not consistently evolved under SRG, which introduces an additional source of uncertainty.

In the meantime, two problematic features of $N\!N\!$+$3N(400)$ were identified, both in the $3N$ sector. The first one concerns the LEC $c_{D}$ entering both the $3N$ contact operator and the contact axial current. This constant was originally fitted, via its contribution to the axial current in $\beta$-decay processes, to the triton half-life Gazit et al. (2009). Recently, an error in the derivation of the axial current was pointed out Schiavilla (2018); Gazit et al. (2019), which raises questions about the pertinence of the previous fit, as well as about subsequent many-body calculations that employed this interaction. The second point relates to possible artefacts introduced when certain types of regulator functions are used. In principle, any function suppressing high momentum modes could be employed. In practice, however, it was shown Hagen et al. (2014b); Dyhdalo et al. (2016); Lynn et al. (2017); Lu et al. (2019) that local regulators, as employed for $3N$ operators in the standard $N\!N\!$+$3N(400)$ implementation, are likely to cut off regions of the phase space that are instead relevant. The success of NNLO${}_{\text{sat}}$was also attributed in part to the use of nonlocal regulators in the $3N$ interaction.

To overcome the above difficulties, a variant of the $N\!N\!$+$3N(400)$ interaction that remedies the contact axial current fit of the LECs and employs both local and nonlocal (lnl) $3N$ regulators is introduced and tested in the present work. This new Hamiltonian is labelled $N\!N\!$+$3N\text{(lnl)}$. Low-energy constants of $2N$ and $3N$ operators are fitted to properties of $A=2,3,4$ nuclei, for which an excellent description is maintained including the triton half-life once the correct contact axial current is employed. Properties of light nuclei are investigated by means of no-core shell model (NCSM) Barrett et al. (2013) calculations. Encouragingly, results are generally in good agreement with experiment. In order to assess its quality in medium-mass nuclei, many-body calculations are performed within the self-consistent Green’s function approach, in both its closed-shell (i.e. Dyson) Cipollone et al. (2013) and open-shell (i.e. Gorkov) Somà et al. (2011) versions. Three representative isotopic chains, namely oxygen, calcium and nickel, are addressed. Total binding energies, two-neutron separation energies, charge radii and density distributions, as well as one-nucleon addition and removal spectra are systematically analysed and compared with calculations performed with other Hamiltonians. Overall, for what concerns ground- and excited-state energies, the performance of $N\!N\!$+$3N\text{(lnl)}$ is very satisfactory, in line with NNLO${}_{\text{sat}}$ or even superior in certain regions. Only for radii (and densities) does the picture change, with $N\!N\!$+$3N\text{(lnl)}$ improving with respect to $N\!N\!$+$3N(400)$ but still underestimating experimental data and NNLO${}_{\text{sat}}$ results. Possible sources of such disagreement are discussed.

The paper is organised as follows. The $2N+3N$ Hamiltonians used in the present study are described in Sec. II, with particular emphasis on the novel $N\!N\!$+$3N\text{(lnl)}$ interaction. Section III discusses the performance of the new Hamiltonian in light systems, focusing on total energies and excitation spectra across $s$- and $p$-shell nuclei. Results for medium-mass nuclei are presented in Sec. IV. First, the employed many-body method is introduced (Sec. IV.1) and convergence of the calculations with respect to model space (Sec. IV.2) and many-body truncations (Sec. IV.3) are discussed. Next, results for several observables are systematically studied: ground-state energies (Sec. IV.4), charge radii and distributions (Sec. IV.5), and spectra of odd-$A$ nuclei (Sec. IV.6). Concluding remarks are presented in Sec. V.

II Hamiltonians

Three different Hamiltonians are employed in this work. The first one, labelled $N\!N\!$+$3N(400)$, is based on the chiral N3LO nucleon-nucleon potential from Entem and Machleidt Entem and Machleidt (2003); Machleidt and Entem (2011) combined with the chiral N2LO $3N$ interaction with a local regulator Navrátil (2007). The $2N$ interaction of Ref. Entem and Machleidt (2003) was built with a cutoff of 500 MeV/c, however, a 400 MeV/c regulator was used for the $3N$ sector Roth et al. (2012). This Hamiltonian has been employed extensively in calculations of $p$- and $sd$-shell nuclei and describes well the binding energy of oxygen, nitrogen and fluorine isotopes Hergert et al. (2013); Cipollone et al. (2013). Nevertheless, as also shown in this paper, it leads to overbinding in medium-mass nuclei starting in the calcium chain and underpredicts radii even for O isotopes Binder et al. (2014); Somà et al. (2014a); Lapoux et al. (2016). In this Hamiltonian, the $3N$ low-energy constant $c_{D}$ had been set to -0.2 based on a 3H half-life fit in Ref. Gazit et al. (2009). However, recently it was pointed out that an error was present in the relationship between the two-body axial current LEC $d_{R}$ and the $c_{D}$ (specifically, a missing factor of -1/4) and that the $c_{D}$ that actually fits the 3H half-life is 0.83(24) Gazit et al. (2019). The $c_{E}=0.098$ LEC was determined by fitting the 4He binding energy Roth et al. (2012) after the cutoff was reduced to 400 MeV/c from the original 500 MeV/c Gazit et al. (2009) to mitigate four-body terms induced by the SRG evolution Roth et al. (2012).

With the main goal of improving the description of radii in medium-mass nuclei, a new chiral Hamiltonian with terms up to N2LO was developed in Ref. Ekström et al. (2015). It is characterised by a simultaneous fit of $2N$ and $3N$ LECs that does not rely solely on two-nucleon and $A$=3,4 data, but also on binding energies of 14C and 16,22,24,25O as well as charge radii of 14C and 16O. The resulting interaction, named NNLO${}_{\text{sat}}$, successfully describes the saturation of infinite nuclear matter Ekström et al. (2015), the proton radius of 48Ca Hagen et al. (2016a) and the nuclear radii of neutron-rich carbon isotopes as well as other medium mass nuclei Kanungo et al. (2016); Lapoux et al. (2016). It also performed well in several other applications, including the description of the parity inversion in 11Be Calci et al. (2016), the 10C($p,p$)10C elastic scattering Kumar et al. (2017), electron scattering Rocco and Barbieri (2018); Barbieri et al. (2019), giant dipole resonances Raimondi and Barbieri (2019), and the derivation of microscopic optical potentials Idini et al. (2019). Unlike the $N\!N\!$+$3N(400)$ interaction, NNLO${}_{\text{sat}}$ employs a non-local regulator.

Motivated by the success of NNLO${}_{\text{sat}}$, the objective of the present work is to amend the original $N\!N\!$+$3N(400)$ interaction, and in particular its $3N$ part. While the latter has been shown to be problematic, its $2N$ part is instead believed to perform relatively well and thus is kept unchanged. Being based on the N3LO potential, which provides a better description of nucleon-nucleon data compared to the lower-order NNLO${}_{\text{sat}}$, it guarantees superior features in light systems, e.g. a better reproduction of spectroscopy of natural parity states in $p$- and light $sd$-shell nuclei. A comparison of the different Hamiltonians in the calculation of the ground-state energy of 4He is shown in Table 1, where mean values of the five $3N$ N2LO terms are displayed. Curiously, for $N\!N\!$+$3N(400)$ and NNLO${}_{\text{sat}}$, $\langle c_{3}\rangle$, $\langle c_{D}\rangle$, and $\langle c_{E}\rangle$ terms contribute with opposite sign. This is particularly disturbing for the three-nucleon contact term, $c_{E}$. In fact, one might argue that this could be at the origin of the severe overbinding generated by $N\!N\!$+$3N(400)$ in medium-mass nuclei. Consequently, here the $c_{D}$ and $c_{E}$ LECs are changed to 0.7 and -0.06, respectively, to get 4He results more in line with the ones of NNLO${}_{\text{sat}}$. In addition, the regulator of the $3N$ interaction was modified by introducing a non-local regulator of the same type as that used in the NNLO${}_{\text{sat}}$ on top of the local one employed in $N\!N\!$+$3N(400)$. For technical reasons, a completeness in the three-nucleon antisymmetrized harmonic-oscillator (HO) basis was applied to evaluate the matrix elements of the new $3N$ interaction. The non-local regulator was set to 500 MeV/c to be consistent with the cutoff of the N3LO $2N$ interaction and the local regulator was increased from 400 MeV/c to 650 MeV/c to achieve a larger binding in A=3,4 systems in agreement with experiment. With these LECs, the new interaction reproduces very reasonably experimental ground-state energies of 3H, 3He and 4He, as well as the 3H half-life.

The performance of this new interaction, named $N\!N\!$+$3N\text{(lnl)}$, in 4He is shown in Table 1. One observes that all $3N$ terms contribute with the same sign as in NNLO${}_{\text{sat}}$ and the $\langle c_{E}\rangle$ is about the same. The $N\!N\!$+$3N\text{(lnl)}$ interaction was already applied in Ref. Leistenschneider et al. (2018) to the description of binding energies of neutron-rich titanium isotopes and neighbouring isotopic chains. Later on, it was employed in the calculation of beta decays of selected light and medium mass nuclei Gysbers et al. (2019), where it was denoted NN-N3LO+3Nlnl. Very recently, it was also used in a study of quasifree nucleon knockout from 54Ca Chen et al. (2019). Its performance in these preliminary applications was very promising, which is confirmed in the present systematic and extensive study. It turns out, however, that the good quality of this new interaction (and of NNLO${}_{\text{sat}}$) in medium-mass nuclei is not just a consequence of a particular choice of the LECs $c_{D}$ and/or $c_{E}$. In Table 1, results of the corrected $2N$+$3N(500)_{c_{D}\!0.83}$ interaction from Ref. Gazit et al. (2019) are also shown. As compared to $N\!N\!$+$3N\text{(lnl)}$, this Hamiltonian has the same $c_{1},c_{3},c_{4}$ LECs (i.e. those from the $2N$ N3LO interaction Entem and Machleidt (2003) also used in $N\!N\!$+$3N(400)$) and almost identical $c_{D}=0.83$ and $c_{E}=-0.052$, but it employs a local regulator as the $N\!N\!$+$3N(400)$ although with a 500 MeV/c cutoff. One can see that it gives very similar results in 4He as the $N\!N\!$+$3N\text{(lnl)}$. However, it severely overbinds medium-mass nuclei in many-body perturbation theory calculations Miyagi (2019). Consequently, the choice of the regulator appears crucial for the correct description beyond light nuclei and perhaps a hint of the reason of the superior performance of $N\!N\!$+$3N\text{(lnl)}$ manifests in the 4He results by the absolute values of $\langle c_{3}\rangle$ and $\langle c_{4}\rangle$ reduced while the $\langle c_{E}\rangle$ contribution enhanced compared to $2N$+$3N(500)_{c_{D}\!0.83}$. Table 2 summarises the LECs used to construct the 3N interactions of all the above Hamiltonians.

III Light nuclei

As an initial test of the $N\!N\!$+$3N\text{(lnl)}$ Hamiltonian, NCSM Barrett et al. (2013) calculations of energies of $s$-shell and selected $p$-shell nuclei were performed. In the NCSM, nuclei are considered to be systems of $A$ nonrelativistic point-like nucleons interacting via realistic two- and three-body interactions. Each nucleon is an active degree of freedom and the translational invariance of observables, the angular momentum, and the parity of the nucleus are conserved. The many-body wave function is expanded over a basis of antisymmetric $A$-nucleon harmonic oscillator (HO) states. The basis contains up to $N_{\rm max}$ HO excitations above the lowest possible Pauli configuration, so that the motion of the center of mass is fully decoupled and its kinetic energy can be subtracted exactly. The basis is characterised by an additional parameter $\Omega$, the frequency of the HO well, and may depend on either Jacobi relative Navrátil et al. (2000) or single-particle coordinates Navrátil et al. (2000). The convergence of the HO expansion can be greatly accelerated by applying an SRG transformation on the $2N$ and $3N$ interactions Wegner (1994); Bogner et al. (2007); Hergert and Roth (2007); Bogner et al. (2010); Jurgenson et al. (2009). Except for $A{=}3,4$ nuclei, here and in the following of the paper an SRG evolution is applied to the $N\!N\!$+$3N(400)$ and $N\!N\!$+$3N\text{(lnl)}$ interactions down to a scale of $\lambda{=}2$ fm*-1*. On the contrary, calculations with NNLO${}_{\text{sat}}$ are performed with the bare Hamiltonian.

In Figs. 1, 2 and 3 the excitation energy spectra of selected Li, Be, B, and C isotopes are displayed. A correct ordering of low-lying levels is found for all the considered lithium and beryllium isotopes, namely 6,7,9Li and 8,9Be. The $2^{+}0$ and $1^{+}_{2}0$ states in 6Li as well as some of the excited states in 7Li and 8,9Be are broad resonances. Here a more realistic description of 6Li and 9Be would require a better treatment of continuum effects, see Refs. Hupin et al. (2015) and Langhammer et al. (2015), respectively, in this regard. Let us note that all excited states of 6Li are unbound with respect to the emission of an $\alpha$ particle and that 7Li has only one excited state below the $\alpha$-separation threshold. Similarly, 8Be is never bound and even its ground state is unstable against decay into two $\alpha$. The lowest states in 10B are known to be very sensitive to the details of nuclear forces, and the $3N$ interaction in particular Navrátil et al. (2007). Here a good description is achieved by $N\!N\!$+$3N\text{(lnl)}$, with only the $1^{+}_{2}0$ state resulting incorrectly placed. The correct level ordering is also found in 11B, with the spectrum being overall too compressed as compared to the experimental one. Finally, worth-noting is the correct ordering of $T{=}1$ states in 12C, also known to be sensitive to the $3N$ interaction. On the other hand, the alpha-cluster dominated $0^{+}0$ Hoyle state in 12C cannot be reproduced in the limited NCSM basis employed here Epelbaum et al. (2011). In general, $N\!N\!$+$3N\text{(lnl)}$ yields spectra that are in good agreement with experiment. Some underestimation of level-splitting in 9Li, 11B, and 13C emerges, and could be associated with a weaker spin-orbit interaction strength. This is comparable to what has been found with earlier parameterisations of chiral $3N$ forces (see, e.g., Ref. Navrátil et al. (2007)).

Ground-state energies of 3H, 3,4He, and selected $p$-shell nuclei from 6He to 16O are shown in Fig. 4. The calculated values (red lines) obtained with the $N\!N\!$+$3N\text{(lnl)}$ interaction are compared to experiment (blue lines). Theoretical error bars represent the NCSM extrapolation uncertainty. Overall, experimental data are very reasonably reproduced, with differences of at most a few percent. The agreement is the best for $T_{z}{=}0$ and $|T_{z}|{=}1/2$ nuclei. Some deficiency of the interaction is observed with increasing $|T_{z}|$; for example, while 4He is in a perfect agreement with experiment, the 6He and 8He are barely bound. Note however that, again, a proper treatment of continuum effects, not included here, is likely to provide additional binding to systems close to the dripline Calci et al. (2016). Overall, the performance of the $N\!N\!$+$3N\text{(lnl)}$ Hamiltonian in light nuclei is very encouraging.

IV Medium-mass nuclei

IV.1 Self-consistent Green’s function theory

In standard, i.e. Dyson, self-consistent Green’s function theory (DSCGF) Dickhoff and Barbieri (2004); Barbieri and Carbone (2017), the solution of the $A$-body Schrödinger equation is achieved via its rewriting in terms of one-, two-, …, $A$-body objects named propagators or, indeed, Green’s functions (GFs). Green’s functions are expanded in a perturbative series, which in self-consistent schemes is recast in terms of the exact GFs so that a large portion of non-skeleton diagrams are implicitly resummed. One is mostly interested in the one-body Green’s function since this provides access to all one-body observables and to the ground-state energy via the so-called Galitskii-Migdal-Koltun sum rule Galitskii and Migdal (1958); Koltun (1972). The latter can be properly generalised to account for three-body forces Carbone et al. (2013). In addition, the one-body GF contains information on neighbouring nuclei. Specifically, the residues from its Lehmann representation are related to transition matrix elements for one-nucleon addition and removal, while the poles give direct access to ground and excited states of ($A\pm 1$)-nucleon systems. Note that in all calculations the intrinsic form of the Hamiltonian is employed, i.e. the center-of-mass kinetic energy is subtracted from the start. Since the latter depends on the number of nucleons at play, different calculations are performed with the Hamiltonian corresponding to mass number $A$ or $A\pm 1$ depending on whether ground-state quantities or nucleon addition or removal spectra are computed, as detailed in Ref. Cipollone et al. (2015).

The one-body GF is obtained by solving the Dyson equation that is intrinsically non-perturbative and in which the irreducible self-energy encodes all non-trivial many-body correlations arising from the interactions of a nucleon with the nuclear medium. The self-energy is particularly important since it encodes information on both the $A$-nucleon ground state and the scattering states of the $A+1$ system. Hence it provides a natural ab initio approach for consistent calculations of structure and reactions Barbieri and Jennings (2005); Idini et al. (2019). In this work the self-energy is computed in the so-called algebraic diagrammatic construction [ADC($n$)] approach up to order $n$=3 Schirmer et al. (1983); Barbieri and Carbone (2017); Raimondi and Barbieri (2018). This entails including all perturbative contribution up to $n$-th order plus any additional resummation needed to preserve its spectral representation. At first order, ADC(1) includes only mean-field terms and it is nothing else than the standard Hartree-Fock (HF) approximation. The higher orders, ADC(2) and ADC(3), add dynamical correlations in terms of 2p1h and 2h1p configurations. However, these remain minimally included and non-interacting at ADC(2) while ADC(3) includes infinite-order resummations of both particle-particle/hole-hole and particle-hole ladders. Generally speaking, ADC($n$) defines a truncation scheme that is systematically improvable up to ADC($\infty$), where exact results are recovered by definition.

In Dyson GF theory the diagrammatic expansion builds on top of a reference state that is particle-number conserving and that typically respects spherical symmetry. While such an expansion can suitably address doubly closed-shell systems, it becomes inefficient or even breaks down in open-shell systems due to the degeneracy of the reference state with respect to particle-hole excitations. With the wish to retain the simplicity of a single-reference method, a possible solution consists in working, from the outset, with a symmetry-breaking reference state. In particular, breaking U(1) symmetry associated with particle number conservation111In the case of atomic nuclei proton and neutron numbers are conserved individually, and therefore it is always intended U(1)${}_{N}\otimes$U(1)Z where one of the two or both are broken. while maintaining spherical symmetry allows us to efficiently capture pairing correlations, thus gaining access to (singly) open-shell nuclei.

In this spirit, Ref. Gorkov (1958) generalised DSCGF to a U(1) symmetry-breaking scheme based on the use of a Hartree-Fock-Bogolyubov reference state and we refer to this approach as Gorkov self-consistent Green’s function (GSCGF) theory. The resulting four (two normal and two anomalous) Gorkov propagators can be conveniently recast in a $2\times 2$ matrix notation via Nambu formalism Nambu (1960). Hence, all standard GF equations are rewritten in a Nambu-Gorkov matrix form. Moreover, Dyson diagrammatics can be generalised to a Gorkov framework, with minor complications arising from the presence of the four different one-body propagators Somà et al. (2011). The introduction of a chemical potential guarantees that the number of particles is the correct one on average. Eventually, the broken symmetry has to be restored. While symmetry-restored formalism has been developed for other (post-Hartree-Fock) many-body methods Duguet (2014); Duguet and Signoracci (2017); Qiu et al. (2017, 2018, 2019), it remains to be formulated for GSCGF. GSCGF theory has been recently implemented in the context of nuclear physics within the ADC(2) truncation scheme Somà et al. (2011, 2013, 2014b) . Thus, the present paper reports results of DSCGF calculations up to ADC(3) and GSCGF calculations up to ADC(2)222When closed-shell systems are considered, a Gorkov calculation automatically reduces to a Dyson one..

IV.2 Model-space convergence

The present calculations are performed using a spherical HO model space that includes up to the $e_{\text{{\rm max}}}\equiv$ max $(2n+l)$ = 13 shell. Any $k$-body operator is to be represented in the same space and one should truncate the corresponding $k$-body basis consistently according to $e_{k\text{max}}=k\,e_{\text{max}}$. The matrix elements of one- and two-body operators are always included in full. However, this is not feasible for three-body interactions due to the rapid increase in the number of their matrix elements and therefore these are restricted to $e_{3\text{max}}=16<3\,e_{\text{max}}$. The dependence on the basis parameters was tested by computing ground-state observables for different harmonic oscillator frequencies, $\hbar\Omega$, and model space sizes, $e_{\text{max}}$. Results for ground-state energies, $E$, and root-mean-square (rms) radii, $\langle r^{2}_{\text{ch}}\rangle^{1/2}$, computed at the ADC(2) truncation level are displayed in Fig. 5 for two representative nuclei, 36Ca and 68Ni, and for $N\!N\!$+$3N\text{(lnl)}$ and NNLO${}_{\text{sat}}$ Hamiltonians.

Focusing on ground-state energies, both interactions show a typical convergence pattern consisting in curves that gradually become independent of $\hbar\Omega$ and closer to each other as the basis increases. For both nuclei, the change from $e_{\text{max}}=9$ to $e_{\text{max}}=13$ at the $\hbar\Omega$ minimum is larger for NNLO${}_{\text{sat}}$ than for $N\!N\!$+$3N\text{(lnl)}$, consistently with the SRG-evolved character of the latter. In 36Ca, going from $e_{\text{max}}=11$ to $e_{\text{max}}=13$ results in a 1.5 MeV gain for NNLO${}_{\text{sat}}$ and a 300 keV gain for $N\!N\!$+$3N\text{(lnl)}$. In 68Ni gains are 4.9 and 2 MeV respectively. Note that the basis limitations on three-body forces do not affect the lighter systems considered in this work and $e_{3\text{max}}$ = 16 is normally sufficient to converge isotopes around 40,48Ca. However, this truncation can introduce some uncertainties for heavier masses. We shall quantify these errors in Sec. IV.4, when discussing the neutron-rich nuclei in the $pf$-shell.

Charge radii also show their usual convergence pattern, with the $\hbar\Omega$ dependence decreasing as the model space truncation increases. For $e_{\text{max}}=13$ calculations, the $N\!N\!$+$3N\text{(lnl)}$ curves are already rather flat in an interval including both the energy minimum and smaller values of $\hbar\Omega$. Results with NNLO${}_{\text{sat}}$, on the other hand, still present a manifest $\hbar\Omega$ dependence. Hence, a precise determination of the rms charge radius would require the use of extrapolation techniques. As the conclusions of the large-scale systematic analyses presented in this work would not be impacted, such extrapolation is not performed here and is left for future studies.

IV.3 Many-body convergence

Next, let us investigate convergence with respect to the many-body truncation. In Fig. 6, energies per nucleon and rms charge radii computed within ADC(1), ADC(2) and ADC(3) approximations are displayed for the same two representative cases of 36Ca and 68Ni. A clear convergence pattern is visible in all cases. For ground-state energies, ADC(1) results depend strongly on the interaction, with the softer $N\!N\!$+$3N\text{(lnl)}$ more bound than NNLO${}_{\text{sat}}$.

ADC(2) calculations are already sufficient to grasp the bulk of correlation energy. Interestingly, ADC(2) values from the two interactions are similar, which reflects the ability of this self-consistent scheme to resum relevant many-body contributions even in the presence of an SRG-unevolved (though relatively soft) interaction. Going from ADC(2) to ADC(3) results in a further gain in correlation energy, which shows that the ADC(3) truncation level is necessary for precise estimates of total ground-state energies even when evolved or soft interactions are employed. Quantitatively, when going from ADC(2) to ADC(3), one gains 7.6 MeV (9.2 MeV) absolute energy in 36Ca and 16.2 MeV (20.1 MeV) in 68Ni with $N\!N\!$+$3N\text{(lnl)}$ (NNLO${}_{\text{sat}}$)333This corresponds to a 7.9$\%$ (4.5$\%$) increase of correlation energy in 36Ca and a 9.0$\%$ (4.7$\%$) increase in 68Ni with $N\!N\!$+$3N\text{(lnl)}$ (NNLO${}_{\text{sat}}$).. Extrapolating this convergence sequence, one may expect ADC($n$) with $n\geq 4$ to add as little as $\sim 1\%$ correlation energy. As discussed further in Secs. IV.4 and  IV.5, while corrections from ADC(3) are important when confronting total ground-state energies they tend to remain rather constant across whole isotopic chains. This implies that ADC(2) already yields reliable predictions for trends and differential quantities such as two-nucleon separation energies. On the other hand, the ADC(3) is also known to be important to reproduce affinities and ionization spectra in molecules von Niessen et al. (1984) and, correspondingly, it gives significant corrections to the spectra of dominant quasiparticle states discussed in Sec. IV.6.

In the case of calculations performed in the Gorkov framework, an additional source of error comes from the fact that the broken U(1) symmetry is not presently restored. While the number of particles remains the correct one on average, this leads to a dispersion in $N$ and/or $Z$ depending on the open-shell nature of neutron or protons. Here only calculations in semi-magic nuclei are reported, for which the proton variance remains zero. At the ADC(2) level, the maximum variance in neutron number amounts to $\sigma_{N}^{2}\approx$ 1.2, 1.6 and 1.9 for oxygen, calcium and nickel chains respectively (independently of the employed interaction).

For rms radii the convergence pattern results are even more favourable. While ADC(1) already provides a reasonable account of the charge radius, ADC(2) is necessary to reach an essentially converged value, especially in the case of NNLO${}_{\text{sat}}$. Eventually, ADC(3) adds at most 1.2$\%$ to the rms radii that is well converged with respect to the many-body truncation for both $N\!N\!$+$3N\text{(lnl)}$ and NNLO${}_{\text{sat}}$. The above findings on both spectra and radii are consistent with the analysis performed on 34Si and 36S in Ref. Duguet et al. (2017). The same formalism applied in the context of quantum chemistry also showed a similar ADC convergence behaviour von Niessen et al. (1984); Danovich (2011); Degroote et al. (2011); Barbieri et al. (2012).

Table 3 displays a benchmark between the NCSM and SCGF for 16O and the three Hamiltonians. This isotope is still light enough that it can be computed using the importance-truncated NCSM Roth and Navrátil (2007); Roth (2009), although the presence of matrix elements from the $3N$ interaction limits the largest possible basis truncation to $N_{\rm max}=10$.  This issue is more severe for the harder NNLO${}_{\text{sat}}$ interaction. Thus, the NCSM results have been extrapolated using a standard exponential trend with respect to $N_{\rm max}$ and a polynomial dependence on the importance-truncation parameter $\kappa_{min}$ Roth and Navrátil (2007). Table 3 reports the extrapolated values together with the uncertainties estimated from this procedure by repeating the extrapolations with different subsets of data. Taking these uncertainties under consideration, the comparison among the two ab initio approaches is extremely satisfactory and confirms the reliability of SCGF for the computation of medium-mass isotopes presented in the following. From Table 3, one also notices that the $N\!N\!$+$3N(400)$—which was the first Hamiltonian to successfully predict the oxygen ground state energies ab initio—already displays a slight tendency to overbind, even for this nucleus. The NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$ correct this effect, with the latter performing a bit better in comparison to the experiment.

IV.4 Ground-state energies

The following subsections study the performances of the three Hamiltonians –$N\!N\!$+$3N(400)$, $N\!N\!$+$3N\text{(lnl)}$ and NNLO${}_{\text{sat}}$– along three representative medium-mass isotopic chains, namely oxygen, calcium and nickel. Based on the considerations of Secs. IV.2 and IV.3, all the following calculations are performed with an $e_{\text{max}}=13$ model space (14 shells), $e_{3\text{max}}$ = 16, and oscillator frequencies fixed at $\hbar\Omega=20$ MeV for NNLO${}_{\text{sat}}$ and $\hbar\Omega=18$ MeV for $N\!N\!$+$3N\text{(lnl)}$. Similar studies have shown that $N\!N\!$+$3N(400)$ has an optimal minimum at $\hbar\Omega=28$ MeV Somà et al. (2014b), which is used here for this Hamiltonian. For the three isotopic chains, protons maintain a good closed-shell character, i.e. all isotopes are at least semi-magic, which generally ensures that deformation does not play a major role444A possible exception is represented by some nickel isotopes between 56Ni and 68Ni, as discussed later.. Ground-state properties (total binding energies, charge radii, density distributions) of even-even nuclei as well as excitation spectra of odd-even nuclei are investigated to provide a comprehensive benchmark of the three interactions.

Total ground-state energies of oxygen, calcium and nickel isotopes are displayed in Figs. 7(a), 8(a) and 9(a) respectively. ADC(2) results (coloured points and lines), covering all even-even isotopes, are shown together with ADC(3) calculations in doubly closed-shell nuclei (coloured horizontal bars) and compared to available experimental data (black points). Corresponding two-neutron separation energies are shown in panels (b). Following the analysis of Secs. IV.2 and IV.3, model-space convergence errors for NNLO${}_{\text{sat}}$ ($N\!N\!$+$3N\text{(lnl)}$) are estimated to be at most 1% (0.5%) of the total binding energy up to the calcium isotopes and 2% (1%) for the nickels up to 68Ni. Many-body truncation errors are 4% for ADC(2) and below 1% for ADC(3), generally underestimating the binding energy. Uncertainties for $N\!N\!$+$3N(400)$ are the same as for $N\!N\!$+$3N\text{(lnl)}$.

All three interactions yield similar results for ground-state energies of the oxygen isotopes and are generally close to experimental values. While for $N\!N\!$+$3N(400)$ and NNLO${}_{\text{sat}}$ the agreement is excellent through the whole chain, $N\!N\!$+$3N\text{(lnl)}$ shows some mild underbinding for the most neutron-rich systems. Although additional correlations coming in at the ADC(3) level tend to provide additional binding, one notices that this effect is not large in oxygen. For all interactions the dripline at 24O is correctly reproduced, as also visible in Fig. 7(b). For the model space parameters used here, the two N3LO Hamiltonians predict 28O to be less bound than 26O, while the opposite is found for NNLO${}_{\text{sat}}$. However, we find that computed ground-state energies for the unbound 28O depend sensibly on $e_{\text{max}}$ and $\hbar\Omega$ which is consistent with a discretization of the continuum imposed by the HO space.

For heavier systems like calcium and nickel, the $N\!N\!$+$3N(400)$ Hamiltonian is known to produce strong overbinding with respect to experimental data Binder et al. (2014); Somà et al. (2014a). This is confirmed by present calculations as visible in Figs. 8(a) and 9(a). Instead, one notices that the two most recent Hamiltonians, NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$, largely correct for this overbinding. For instance, on the light-mass side, the ADC(2) energy for 36Ca goes from 20.4 MeV (7.2$\%$) overbinding for $N\!N\!$+$3N(400)$ to 11.8 MeV (4.1$\%$) underbinding for NNLO${}_{\text{sat}}$ and 7.0 MeV (2.4$\%$) underbinding for $N\!N\!$+$3N\text{(lnl)}$. Among the heavier isotopes, 68Ni goes from 64.8 MeV (10.9$\%$) overbinding for $N\!N\!$+$3N(400)$ to 45.0 MeV (7.6$\%$) underbinding for NNLO${}_{\text{sat}}$ and 15.9 MeV (2.6$\%$) underbinding for $N\!N\!$+$3N\text{(lnl)}$.

Many-body correlations beyond ADC(2) provide additional binding and ground-state energies of all considered isotopes are lower by 2-3% when switching to ADC(3). While this aggravates the overbinding of $N\!N\!$+$3N(400)$ Somà et al. (2014a), it is expected to reduce the underbinding of the other two potentials. The latter expectation is corroborated by ADC(3) results of closed-shell nuclei along the two chains. Once ADC(3) corrections are included, binding energies computed with both NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$ Hamiltonians are in excellent agreement with experimental data. For the above examples, differences with experiment reduce to 0.9$\%$ and 0.2$\%$ in 36Ca and to 4.2$\%$ and 0.05$\%$ in 68Ni for NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$ respectively and are comparable with the theoretical uncertainties due to the model-space convergence. The extra binding obtained within ADC(3) is therefore crucial if one is after precise comparisons on total ground-state energies. In accordance with the findings of Sec. IV.3, ADC(3) corrections are systematically larger for the SRG-unevolved NNLO${}_{\text{sat}}$ than for the SRG-evolved $N\!N\!$+$3N\text{(lnl)}$. Note that Ref. Hagen et al. (2016b) already reported very poor convergence of the 78Ni isotope with NNLO${}_{\text{sat}}$, mainly due by the truncation of three-body matrix elements. Here, it is found that changing $e_{3\text{max}}$ from 14 to 16 leads to a $\approx$40 MeV variation, which adds to the model space uncertainties discussed in Sec. IV.2. All isotopes beyond 68Ni are likely to be affected in an analogous way.

Further insight can be gained by looking at energy differences. Two-neutron separation energies along the two chains are displayed in Figs. 8(b) and 9(b). For calcium, the effects of overbinding in $N\!N\!$+$3N(400)$ cancels out to a large extent, with the residual mass dependence showing up in the most proton- and neutron-rich isotopes. Calculations with NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$ lead to similar values and closely follow experimental data. Interestingly, all major gaps correctly emerge. While the $N=20$ gap is quantitatively well reproduced by NNLO${}_{\text{sat}}$, it appears somewhat overestimated with $N\!N\!$+$3N\text{(lnl)}$. The opposite holds for the $N=28$ gap, with $N\!N\!$+$3N\text{(lnl)}$ providing a very accurate description.

One compelling question, attracting much attention both experimentally and theoretically, relates to the position of the dripline in calcium isotopes. Very recently, the experimental knowledge was extended with the mass measurements of 55-57Ca Michimasa et al. (2018) and the first evidence of a bound 60Ca Tarasov et al. (2018). While the present theoretical framework is in principle able to be applied beyond 60Ca, convergence problems were encountered with the present settings and the (inconclusive) results beyond 60Ca are not shown here. In particular, the SRG evolution of three-body operators is performed in a three-body HO space in Jacobi coordinates Navrátil et al. (2000) and requires very high total angular momenta (up to $J$=35/2) to resolve all matrix elements up to $e_{3\text{max}}$ = 16 used in the SCGF model space. Calculations for isotopes above 60Ca were found to be affected by model space truncations in two different ways: one is the direct dependence on $e_{3\text{max}}$ in the SCGF basis and the other is the truncations needed to evolve $N\!N\!$+$3N\text{(lnl)}$ through SRG. As an example, in Fig. 10 total energies of 40Ca and 70Ca are displayed as a function of $e_{3\text{max}}$ for two different ways of truncating the Jacobi basis during the SRG evolution. While calculations in 40Ca are converged with respect to both variations, in 70Ca a clear dependence on the latter parameter is visible. This problem arises only beyond $N=40$ for the two SRG-evolved Hamiltonians and is also responsible for a lack of binding in neutron-rich nickel isotopes, which is reflected in an evident kink in the binding energy curve after 68Ni (see Fig. 9). Future technical improvements as well as a better treatment of the continuum, which appears to be crucial beyond 60Ca, are therefore necessary for a correct determination of the calcium dripline.

Two-neutron separation energies in the nickel chain are generally less accurate. The old $N\!N\!$+$3N(400)$ Hamiltonian struggles to catch the experimental trend, with an unrealistic large gap appearing at $N=32$. The other two interactions show a clear improvement. Interestingly, both of them predict a flat trend for proton-rich isotopes, in contrast with the AME data extrapolation. For $N\!N\!$+$3N\text{(lnl)}$, the agreement with experiment between 56Ni and 64Ni is remarkable. As explained above, after 68Ni the results appear to be affected by convergence issues. For NNLO${}_{\text{sat}}$, the description remains reasonable except for most neutron-rich systems. The reproduction of experimental data also deteriorates between $N=28$ and $N=40$, which is likely to be linked with the onset of deformation although this does not appear to be problematic for the soft $N\!N\!$+$3N\text{(lnl)}$ interaction.

IV.5 Charge radii and density distributions

Next let us examine rms charge radii $\langle r^{2}_{\text{ch}}\rangle^{1/2}$. In the present approach rms charge radii are computed from rms point-proton radii by correcting for the finite charge distributions of protons and neutrons, as well as for the Darwin-Foldy term (see Ref. Cipollone et al. (2015) for details). ADC(2) results for absolute rms charge radii of oxygen, calcium and nickel isotopes are displayed in Figs. 11(a), 12(a) and 13(a), all compared to available experimental data. Relative rms radii $\Delta\langle r^{2}_{\text{ch}}\rangle^{1/2}$, i.e. charge radii differences555Note that the relative rms radii used here, differences of $\langle r^{2}_{\text{ch}}\rangle^{1/2}$, differ from the mean square shifts $\delta\langle r^{2}\rangle$ sometimes found in the literature, defined as differences of $\langle r^{2}_{\text{ch}}\rangle$. relative to a reference isotope, are shown in the corresponding panels (b). Following the analysis of Secs. IV.2 and IV.3, conservative errors for NNLO${}_{\text{sat}}$ ($N\!N\!$+$3N\text{(lnl)}$ and $N\!N\!$+$3N(400)$) are estimated to amount up to 1.8% (1,5%) of the charge radius up to the calcium isotopes and of 2.6% (1%) for the nickels, including 78Ni. Note that these errors are dominated by model space convergence in the case of NNLO${}_{\text{sat}}$ and by the many-body truncations for $N\!N\!$+$3N\text{(lnl)}$ and $N\!N\!$+$3N(400)$.

For what concerns absolute radii, a large variation between the different interactions is observed in all cases. For oxygen isotopes, studies with $N\!N\!$+$3N(400)$ and NNLO${}_{\text{sat}}$ exist in the literature Lapoux et al. (2016); Cipollone et al. (2015), where it was shown that already in these light systems $N\!N\!$+$3N(400)$ leads to a strong underestimation of the size of nuclei. From Fig. 11(a) one notices that the new $N\!N\!$+$3N\text{(lnl)}$ interaction improves on $N\!N\!$+$3N(400)$ results producing charge radii that are $\sim$0.1 fm larger, reducing by a factor 2 or better the discrepancy with experiment. A similar picture emerges from the analysis of calcium and nickel chains. In calcium isotopes, see Fig. 12(a), $N\!N\!$+$3N(400)$ strongly underestimates measured radii, with discrepancies of about 12-15$\%$ along the whole chain. $N\!N\!$+$3N\text{(lnl)}$ significantly improves on $N\!N\!$+$3N(400)$ results producing charge radii that are $\sim$0.3 fm larger. Still, experimental data are underestimated by about 5-6$\%$ across all isotopes. NNLO${}_{\text{sat}}$, on the other hand, succeeds in reproducing the bulk values of Ca rms radii, thus maintaining for this observable the good performances already observed for lighter nuclei Duguet et al. (2017). Importantly, the present results with NNLO${}_{\text{sat}}$ are in good agreement with previous coupled-cluster calculations performed on closed-shell and neighbouring isotopes Garcia Ruiz et al. (2016). Similar conclusions can be inferred by inspecting results for nickel isotopes, reported in Fig. 13(a). Absolute rms charge radii obtained with NNLO${}_{\text{sat}}$ maintain their good agreement with data even for this mass region. A kink is visible at 56Ni, in accordance with its good closed-shell character. Beyond this point, the calculation follow the trend of the limited available data, slightly departing from experiment as neutron number increases. Radii obtained with the other two interactions, again, severely underestimate experimental results. Discrepancies are in line with what observed in calcium isotopes, namely $N\!N\!$+$3N(400)$ is more than 15$\%$ off and $N\!N\!$+$3N\text{(lnl)}$ about 8-9$\%$ off. For all chains, ADC(3) calculations for closed-shell systems are also shown. Contrarily to total ground-state energies, radii are essentially converged at the ADC(2) truncation level and additional ADC(3) correlations do not change the overall picture. In particular, it is clear that the discrepancies with experiment cannot be removed by improving the many-body truncation.

When comparing radii obtained with SRG-evolved and bare Hamiltonians, an important caveat relates to the present omission of potentially relevant higher-body radius operators induced by the SRG transformation. In some recent calculations Schuster et al. (2014); Miyagi et al. (2019), such induced operators have been properly included for the $N\!N\!$+$3N(400)$ interaction but have not led to any sizeable improvement. Thus, although it remains to be seen whether the same holds for $N\!N\!$+$3N\text{(lnl)}$, the discrepancies with the experimental data point to intrinsic deficiencies of the Hamiltonian that will have to be addressed in the future.

The systematic flaws in the underestimation of radii appear to be uniform across each isotopic chain. Thus, one may expect that they cancel out to a good extent in differential quantities putting in evidence the isospin dependence for each interaction. In Fig. 11(b) rms radii differences relative to 16O are shown for oxygen. Clearly, the spread of results is appreciably reduced, with small discrepancies showing up for the less stable systems. Radii along the calcium chain, relative to 48Ca, are shown in Fig. 12(b). The lightest isotopes 34-40Ca are the most sensible to the employed interaction, with rather different trends. Interestingly, recent measurements in 36-38Ca Miller et al. (2019) appear to be in better agreement with $N\!N\!$+$3N\text{(lnl)}$ results rather than the ones obtained with NNLO${}_{\text{sat}}$, which predicts a somewhat steeper slope. Both interactions based on the N3LO two-nucleon force of Ref. Entem and Machleidt (2003) predict an inversion of this trend when going down to mass $A$=34 while NNLO${}_{\text{sat}}$ does the opposite. All calculations roughly reproduce the fact that charge radii for 40Ca and 48Ca are basically the same. However, none of them is capable of accounting for the parabolic behaviour between these two isotopes. This is not surprising since this feature has been associated, in the contexts of particle-vibrations coupling and shell model calculations Barranco and Broglia (1985); Caurier et al. (2001), to the presence of highly collective many particle-many hole configurations that are missing in the many-body approach employed here. Energy density functionals are also striving to reproduce the experimental trend, with only recent applications based on Fayans functionals Reinhard and Nazarewicz (2017) able to capture the peculiar behaviour. After 48Ca, NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$ do improve on the poor trend of $N\!N\!$+$3N(400)$ but still fail to reproduce quantitatively the steep slope leading to 52Ca. The charge radius of the latter, recently measured in laser spectroscopy experiments Garcia Ruiz et al. (2016), thus remains a challenge for many-body calculations. Relative radii in nickel isotopes can be examined in Fig. 13(b). Here two distinct regions can be identified. Below 58Ni, we find a similar behaviour to the one of Ca: $N\!N\!$+$3N(400)$ and $N\!N\!$+$3N\text{(lnl)}$ follow the similar trends and actually predict an increase of rms radii with decreasing neutron number, while NNLO${}_{\text{sat}}$ does the opposite and decreases towards 48Ni. More experimental data on both proton rich Ca and Ni would be very useful to pin down this effect. Above 58Ni, $N\!N\!$+$3N\text{(lnl)}$ and NNLO${}_{\text{sat}}$ predict a very similar behaviour, while $N\!N\!$+$3N(400)$ shows a rather steep increase all the way up to 78Ni. The limited amount of available data gives a stronger support to the former trend. Also in this respect, an extension of our experimental knowledge to some of the neutron-rich nickel isotopes would be very valuable.

An even more stringent test for many-body calculations is represented by the description of the charge density distributions, from which rms charge radii have been traditionally extracted. Charge distributions are experimentally accessed via electron scattering measurements and thus currently limited to stable nuclei, although considerable progress is being made towards extending this technique to unstable systems Suda and Simon (2017). In the present framework the nuclear charge density is computed through the folding of the nuclear point-proton density distribution with the charge density distribution of the proton, see Ref. Duguet et al. (2017) for details. One representative isotope for each of the three chains is studied, namely 16O, 40Ca and 58Ni, respectively displayed in Figs. 14, 15 and 16. In all cases, unsurprisingly, NNLO${}_{\text{sat}}$ calculations are the closest to experiment. For 16O and 40Ca the agreement is remarkable. The description slightly deteriorates for 58Ni where, in spite of an excellent reproduction of the experimental charge radius, the theoretical result mildly deviates from the measured distribution. As for charge radii, $N\!N\!$+$3N\text{(lnl)}$ calculations largely improve on the poor performance of $N\!N\!$+$3N(400)$ but fail to reach the accuracy achieved by NNLO${}_{\text{sat}}$. ADC(3) calculations are also reported (dashed lines) in the case of 16O and 40Ca. As for radii, one concludes that density distributions are largely converged at the ADC(2) level, although ADC(3) correlations do provide a refined description, e.g. in the central region of 16O.

IV.6 Excitation spectra

Spectroscopic properties reflect some of the general features of a Hamiltonian (for instance the ability to reproduce magic gaps) but at the same time are sensitive to finer details, e.g. depending on the spin and parity of the excited state. In Green’s function theory, one-nucleon addition and removal (i.e., separation) energies are naturally accessed from the spectral representation of the one-body propagator, see Refs. Dickhoff and Barbieri (2004); Barbieri and Carbone (2017) for details. The generalisation to Gorkov Green’s functions allows for an analogous spectral form that also contains information on separation energy spectra of odd-even neighbours Somà et al. (2011). While the ADC(2) approximation does introduce dynamical correlations that induce a fragmentation of the mean-field spectral function, one might ask whether such correlations are too crude for a quantitative description of (low-lying) excitation spectra. The ADC(3) truncation scheme, by coupling the bare two particle-one hole (2p1h) and 2h1p (or three-quasiparticle in Gorkov theory) configurations introduced in ADC(2), stabilises dominant quasiparticle peaks, usually compresses the spectra and generates further fragmentation Cipollone et al. (2015); von Niessen et al. (1984).

In order to test the two levels of approximation, one-neutron removal (addition) spectra from (to) 48Ca are studied in detail in Fig. 17 using NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$ Hamiltonians. Starting with one-neutron removal, i.e. states in 47Ca, one first notices that, for both interactions, ADC(2) spectra are too spread out, with the first excited states at 5-6 MeV to be compared with about 2 MeV in experiment. Since such states are associated to the removal of a neutron in the $sd$-shell, this is a direct consequence of the overestimation of the $N=20$ gap, see Fig. 8(b). Note that the overestimation is more severe in $N\!N\!$+$3N\text{(lnl)}$ calculations, which is reflected in higher excitation energies and a larger splitting between the $1/2^{+}$ and $3/2^{+}$ states as dictated by its underestimation of radii. Including ADC(3) correlations helps in compressing the spectrum, although the effect of an overestimated $N=20$ gap remains. Interestingly, in NNLO${}_{\text{sat}}$ calculations the correct ordering of $1/2^{+}$ and $3/2^{+}$ states is re-established. In addition, negative-parity states $1/2^{-}$ and $3/2^{-}$ appear. As opposed to the positive-parity levels that are obtained as a simple removal from the $sd$-shell, such states correspond to more complex configurations involving particle-hole excitations across the $N=28$ gap and are not captured by the simpler ADC(2) approximation.

The situation is different for one-neutron addition spectra, displayed in Figs. 17(c) and 17(d). Here low-lying states computed at the ADC(2) level are already in reasonably good agreement with experimental values for both interactions. Again, the quality of the description is correlated with the (excellent) reproduction of the $N=28$ and $N=32$ gaps over which the excitations take place. ADC(3) correlations mainly impact NNLO${}_{\text{sat}}$ results, with the position and splitting of $1/2^{-}$ and $5/2^{-}$ changing by a few hundreds keV and moving closer to data. An exception is represented by the $9/2^{+}$ state, which is high-lying in ADC(2) and gets shifted down by a few MeV in ADC(3).

Following this analysis let us turn to the case of 54Ca and look at all one-nucleon removal and addition spectra, i.e. at its four possible odd-even neighbours. These isotopes are of significant interest, with experiments that either have been recently performed or that are planned for the near future, and should complement the currently scarse data in the region. Results from both ADC(2) and ADC(3) calculations with NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$ are reported in Fig. 18. In 53Ca two excited states have been measured around 2 MeV with tentative spin-parity assignments of $5/2^{-}$ and $3/2^{-}$ Perrot et al. (2006); Steppenbeck et al. (2013). Both interactions yield the two states and support the spin assignments. However, $N\!N\!$+$3N\text{(lnl)}$ does a better job in reproducing both the position and the energy splitting between them. In 55Ca, in addition to one-neutron addition states to 54Ca, one-neutron removal states from 56Ca in the ADC(2) approximation are also shown. The spectra generated by the two interactions display the same low-lying states, although the one from NNLO${}_{\text{sat}}$ is more compressed than the one from $N\!N\!$+$3N\text{(lnl)}$. In both cases the excited state corresponding to the main one-neutron addition quasiparticle, with spin-parity $9/2^{+}$, shows a large correction from ADC(3). In 53K the two Hamiltonians predict a different ground state, with $N\!N\!$+$3N\text{(lnl)}$ agreeing with the tentative experimental assignment. Finally, in 55Sc it is the first excited states to be different, with NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$ predicting, respectively, $3/2^{-}$ and $5/2^{-}$ states on top of the $7/2^{-}$ ground state.

The identification of the ground-state spin in 53K is of particular interest and is related to a series of present-day experimental efforts along potassium isotopes. From 37K up to 45K, the ground-state spins have been known to be $3/2^{+}$, as a naive shell model picture would suggest. Several years back, 47K was shown to have a $1/2^{+}$ ground state via a laser spectroscopy experiment Touchard et al. (1982), with $3/2^{+}$ becoming a low-lying excited state at 360 keV. Recently, high-resolution collinear laser spectroscopy measurements determined that the ground-state spin inversion is maintained in 49K but a re-inversion occurs for 51K Papuga et al. (2013). Available experimental data are summarised in Fig. 19, where the energy difference between $1/2^{+}$ and $3/2^{+}$ states is displayed for even $N$ potassium isotopes. At the time, GSCGF calculations were performed with the $N\!N\!$+$3N(400)$ interaction Papuga et al. (2014), which resulted in the red curve reported in Fig. 19. Although the calculations parallel the experimental trend, the energy gap between the two states is largely overestimated, and the spin inversion in 47K is absent. The same observables have been computed here using the two more recent interactions. NNLO${}_{\text{sat}}$ captures the trend as $N$ increases, but presents a shift compared to data that generates the inversion already at 43K. After that, the ground-state is predicted to have always spin-parity $1/2^{+}$. Instead, $N\!N\!$+$3N\text{(lnl)}$ succeeds in reproducing experimental data, including the inversion and re-inversion of the ground-state spin-parity and the position of the first excited state with remarkable accuracy. Note that Fig. 19 displays results at the ADC(2) level. The ADC(3) corrections for the 53K gap, from Fig. 18, are of at most 0.4 MeV and suggest that missing many-body truncations could shift slightly these curves but are unlikely to alter our conclusions.

V Discussion and conclusions

Figure 20 summarises the performance of the three Hamiltonians on the different observables considered in the present work. Representative ground-state energies, rms charge radii and one-nucleon separation energies are displayed. The older $N\!N\!$+$3N(400)$ interaction served as a workhorse in the early applications of ab initio calculations with chiral $2N+3N$ forces, with empirical success in light isotopes, up to oxygen isotopes and neighbouring elements. In particular, a notable achievement was the correct reproduction of the oxygen dripline Otsuka et al. (2010); Hergert et al. (2013); Cipollone et al. (2013). Typically employed in combination with SRG evolution, it has allowed important benchmarks between many-body methods, both non-perturbative and perturbative, that gave the practitioners confidence in the quality of the different many-body approximations Hebeler et al. (2015). However, its flaws appeared evident early on with strong overbinding being generated as mass number increases Binder et al. (2014) and a severe underestimation of nuclear radii even for oxygen isotopes Lapoux et al. (2016). Furthermore, formal issues were recently raised Schiavilla (2018); Gazit et al. (2019), which question the consistency of the calculations that employed this interaction. The overbinding and underestimation of radii emerging in $N\!N\!$+$3N(400)$ calculations are clearly visible in Fig. 20.

Present results also confirm the overall empirical quality of the NNLO${}_{\text{sat}}$ interaction. Ground-state energies are well reproduced even beyond the light nuclei that were used in the fit of its coupling constants. A mild underbinding is observed for heavier nickel isotopes. However, before drawing definitive conclusions a careful study of model-space convergence (in terms of the one- and three-body truncations $e_{\text{max}}$ and $e_{3\text{max}}$, see also Figs. 5 and 10) has to be performed. Nevertheless, contributions to charge radii are excellently reproduced by this interaction, even for nickel isotopes. The fact that radii, their associated density distributions and trends in the ground-state energies are already well converged at the ADC(2) opens the way to systematic calculations of full isotopic chains within the GSCGF approach. Concerning spectroscopic properties, NNLO${}_{\text{sat}}$ had proven very accurate in the neutron $pf$-shell for 34Si and 36S Duguet et al. (2017). This is to a good extent confirmed here in the neutron addition and removal spectra of 52Ca and 54Ca. The agreement with experiment however deteriorates when looking at the $sd$-shell below $N=20$, with this Hamiltonian struggling to reproduce the observed inversion and re-inversion of the ground and first excited states along potassium isotopes.

So far, the novel $N\!N\!$+$3N\text{(lnl)}$ interaction had been applied only to specific cases Leistenschneider et al. (2018); Gysbers et al. (2019); Chen et al. (2019), but never tested in a systematic way. In the present work its main ground-state properties as well as some selected excitation spectra have been studied extensively in light and medium-mass nuclei. Results in light systems are very encouraging, with NCSM calculations in overall good agreement with experiment even for spectra that are known to be particularly sensitive to nuclear forces. Total energies are well reproduced across the whole light sector of the nuclear chart. In medium-mass nuclei, present calculations focused on three representative isotopic chains. Total binding energies are found to be in remarkable agreement with experimental values all the way up to nickel isotopes once ADC(3) correlations are included, thus correcting for the overbinding generated with $N\!N\!$+$3N(400)$. ADC(2) calculations of differential quantities, where ADC(3) contributions essentially cancel out, are also very satisfactory and are able to capture main trends and magic gaps in two-neutron separation energies along all three chains. As evidenced in Fig. 20, although largely improving on $N\!N\!$+$3N(400)$, rms charge radii obtained with the $N\!N\!$+$3N\text{(lnl)}$ interaction still underestimate experimental results and do not reach the quality of NNLO${}_{\text{sat}}$. On the other hand this interaction yields an excellent spectroscopy, also where NNLO${}_{\text{sat}}$ strives to give even a qualitatively correct account of experimental data. One-nucleon addition and removal spectra in neutron-rich calcium are well reproduced. Impressively, the evolution of the energy differences between the ground and first excited states along potassium isotopes follows closely the experimental measurements.

Further insight can be gained by gauging the importance of $3N$ operators in the two interactions. In Fig. 21 the ratio of $3N$ over $2N$ contributions to the total energy is displayed for a selection of nuclei as a function of mass number $A$ for NNLO${}_{\text{sat}}$ and $N\!N\!$+$3N\text{(lnl)}$. In the former, $3N$ operators are much more relevant, reaching almost 20% of the $2N$ contribution in heavier systems. In contrast, the ratio stays rather low, around 5%, for $N\!N\!$+$3N\text{(lnl)}$. This has first of all practical consequences, as in the majority of many-body calculations the treatment of $3N$ operators is usually not exact, following either a normal-ordered two-body approximation (see, e.g., Ref. Roth et al. (2012)) or some generalisation of it Carbone et al. (2013). Hence a strong $3N$ component is in general not desirable. On top of that, one might worry about the hierarchy of many-body forces from the standpoint of EFT, and possible need to include subleading $3N$ or $4N$ operators that could have a sizeable effect.

Finally, let us compare $N\!N\!$+$3N\text{(lnl)}$ and NNLO${}_{\text{sat}}$ to an interaction that has been extensively employed in nuclear structure studies in the past few years. Usually labelled as $1.8/2.0$ (EM) and first introduced in Ref. Hebeler et al. (2011), it has proven to yield an accurate reproduction of ground-state energies (as well as low-energy excitation spectra) over a wide range of nuclei Simonis et al. (2016, 2017); Gysbers et al. (2019); Holt et al. (2019). Furthermore, it leads to a satisfactory description of infinite nuclear matter properties Hebeler et al. (2011); Drischler et al. (2014, 2019). In Fig. 22 binding energies per particle obtained within in-medium similarity renormalisation group (IM-SRG) calculations with the $1.8/2.0$ (EM) interaction Simonis et al. (2017) are compared, for a set of closed-shell systems, to the ones computed at the ADC(3) level with $N\!N\!$+$3N\text{(lnl)}$ and NNLO${}_{\text{sat}}$. The three sets of calculations achieve an overall excellent reproduction of experimental data. While NNLO${}_{\text{sat}}$ is superior in light nuclei, it tends to slightly underbind some of the heavier systems. One also notices a striking resemblance of the results obtained with the $N\!N\!$+$3N\text{(lnl)}$ and $1.8/2.0$ (EM) interactions (with the only exception of 60Ca, for which no experimental measurement exists) all the way up to 68Ni. The two potentials indeed present several similarities. First, the bare $NN$ part is the same Entem and Machleidt (2003), even though $N\!N\!$+$3N\text{(lnl)}$ and $1.8/2.0$ (EM) are subsequently SRG-evolved to different scales, $\lambda{=}2$ fm*-1* and $\lambda{=}1.8$ fm*-1* respectively. Second, the $3N$ part builds on N2LO operators and, in the case of $1.8/2.0$ (EM), non-local regulators are applied. A difference comes from the fact that for $1.8/2.0$ (EM) $3N$ forces are not SRG-evolved consistently with the $NN$ operators, but rather the LECs of the three-body contact terms are re-fitted a posteriori to the energy of 3H and the charge radius of 4He. In the end, this results in values ($c_{D}=1.264$ and $c_{E}=-0.120$) that are not very different from the ones of $N\!N\!$+$3N\text{(lnl)}$ (see Table 2).

The present systematic analysis shows that the novel $N\!N\!$+$3N\text{(lnl)}$ Hamiltonian represents a promising alternative to existing nuclear interactions. In particular, it has the favourable features of (i) being adjusted solely on $A=2,3,4$ systems, thus complying with the ab initio strategy, (ii) yielding an excellent reproduction of experimental energies all the way from light to medium-heavy nuclei and (iii) behaving well under similarity renormalisation group transformations, with small induced four-nucleon forces, thus allowing calculations up to medium-heavy systems with moderate computational costs. A first large-scale application with SCGF calculations along few isotopic chains around $Z=20$ is already underway and confirms its excellent phenomenological properties Somà et al. (2019). In the short term, having such high-quality interactions at hand allows us to make useful predictions and to test in depth existing and forthcoming many-body methods. In the long term, such efforts aim to contribute to the long-standing goal of performing simulations of atomic nuclei with fully controlled theoretical uncertainties.

Acknowledgements

We thank the authors of Ref. Simonis et al. (2017) for providing the data included in Fig. 22. This work was supported by the NSERC Grant No. SAPIN-2016-00033, by the Espace de Structure et de réactions Nucléaires Théorique (ESNT) at CEA and by the United Kingdom Science and Technology Facilities Council (STFC) under Grants No. ST/P005314/1 and No. ST/L005816/1. Computing support for P.N. came from an INCITE Award on the Titan supercomputer of the Oak Ridge Leadership Computing Facility (OLCF) at ORNL and from Westgrid and Compute Canada. SCGF calculations were performed by using HPC resources from GENCI-TGCC (Contracts No. A003057392 and A005057392) and at the DiRAC the DiRAC DiAL system at the University of Leicester, UK, (BIS National E-infrastructure Capital Grant No. ST/K000373/1 and STFC Grant No. ST/K0003259/1). TRIUMF receives federal funding via a contribution agreement with the National Research Council of Canada.