Self-consistent Green's function theory for atomic nuclei

Nuclear structure theory has recently gone through a major renewal with the development of ab initio techniques that can be applied to a large number of atomic nuclei, well beyond the light sector that had been traditionally targeted in the past. Self-consistent Green's function theory is one among these techniques. The present work aims to give an overview of the self-consistent Green's function approach for atomic nuclei, including examples of recent applications and a discussion on the perspectives for extending the method to nuclear reactions, doubly open-shell systems and heavy nuclei.


INTRODUCTION
The theoretical description of atomic nuclei is particularly challenging, for several reasons. First, different energy scales are at play, which is manifest in the rich set of observables and phenomena one is confronted with 1 . As a consequence, the choice of relevant degrees of freedom might not be universal but depend on the particular properties one is interested in. The standard description in terms of nucleons, i.e. protons and neutrons, leads to a many-body Schrödinger equation for up to a few hundred particles, which are too many to be easily treated exactly but too few to be amenable to a statistical treatment 2 . Furthermore, interactions between nucleons should be derived from quantum chromodynamics (QCD) in its nonperturbative regime, which prevents direct calculations and requires an additional modelling effort 3 . In the past, all these features have hindered full solutions of the many-body Schrödinger equation and favoured the development of a plentitude of models, following different strategies and characterised by different ranges of application. Although successful in reproducing experimental observations, these models are not always comparable to each other and do not provide a coherent and unified description of nuclear systems.
Only relatively recently so-called ab initio calculations of atomic nuclei, i.e. systematically improvable solutions of the A-body Schrödinger equation that start solely from the knowledge of inter-nucleon interactions, have become available for a large number of isotopes. These advances were made possible by the concomitance of different factors. First, new formal developments of many-body techniques were carried out. Second, chiral effective field theory (χ-EFT) was introduced in nuclear physics [3,4,5], providing a systematic and consistent framework in which the nuclear Hamiltonian can be modelled. Third, similarity renormalisation group (SRG) techniques were applied to χ-EFT Hamiltonians [6], which V. Somà nuclei, going from the oxygen region [34,35] to silicon and sulfur [36], calcium [37,38] and nickel [39] regions, up to the very recent computation of tin and xenon isotopes [40]. As discussed in Secs. 2 and 4.2, excited spectra of odd-even nuclei are also easily accessible in SCGF theory and were studied in particular in Refs. [35,36,39]. Excited states of even-even systems were addressed in the form of the electromagnetic dipole response in Ref. [41]. In Ref. [42], ab initio optical potentials were computed and applied in elastic scattering off 16 O and 40 Ca. Ref [43] instead discussed lepton-nucleus scattering, with particular attention to neutrino scattering off 40 Ar. Dedicated applications studied effective charges (typically employed in shell-model calculations) [44] and the scale dependence of effective single-particle energies and other non-observable quantities [45]. The possibility of using nucleon-nucleon interactions derived from lattice QCD calculations was also explored in Ref. [2]. Last but not least, several collaborations with experimental groups have led to testing the method on e.g. state-of-the-art measurements of nuclear masses [46,47], energies and spins of excited states [48,49,50,51], charge and matter radii [52].
The present article discusses a few of these examples, without any pretension of exhaustivity but with the aim of giving the reader a perception of the reach and versatility of SCGF method, as well we the great deal of possible applications to atomic nuclei. A short introduction to the formalism and implementation in finite nuclei is presented beforehand; however, again, it is far from being complete. The reader interested in a comprehensive treatment of GF theory in a modern form is referred e.g. to the book of Dickhoff and Van Neck [53]. Older, yet valuable sources are the books of Nozières [54] and Abrikosov, Gorkov and Dzyaloshinski [55]. An interesting work, although tailored to solid-state physics, is the one of Economou [56]. An extensive review covering SCGF applications to nuclear physics appeared in 2004 [27]. A more recent pedagogical introduction to the basics of SCGF formalism and implementation in both finite nuclei and infinite nuclear matter (including computational details and examples of numerical codes) can be found in Ref. [57]. A numerical code including a second-order evaluation of the self-energy in the Dyson framework is publicly available [58].
This manuscript is organised as follows. In Sec. 2 the most important concepts and equations of GF theory are introduced. Section 3 describes the actual implementation of the methods in calculations of finite nuclei, briefly reviewing the most commonly used self-energy truncations, the working equations, the choice of basis and interaction. In Sec. 4 some representative applications to ground-state properties, excitation spectra and lepton-nucleus collisions are discussed. Finally, section 5 addresses perspectives and present as well as future challenges, focusing on three directions of research: the consistent treatment of nuclear structure and reactions, the generalisation to doubly open-shell systems and the extension to heavy nuclei.

BASIC CONCEPTS AND EQUATIONS
Many-body Green's function 4 theory comprises a set of techniques that originated in quantum field theory (QFT) and have been subsequently imported in the (non-relativistic) quantum many-body problem. The late 1950's and the 1960's marked the beginning of the field, with flow of QFT ideas and development of formalism. Since the 1970's technical developments were realised and the approach was applied throughout several disciplines and types of many-body problems, ranging from many-electron systems such as crystals, molecules and atoms to many-nucleon systems such as nuclei and nuclear matter. Starting from the 1990's such techniques were implemented as an ab initio method in nuclear physics.

V. Somà
As for other many-body methods, the purpose of such techniques is to achieve an approximate (yet systematically improvable) solution of the A-body Schrödinger equation. In standard (or Dyson) many-body Green's function theory, this is realised by rewriting the Schrödinger equation in terms of one-, two-, ..., A-body objects g I (= g), g II , ..., g A named propagators or, indeed, Green's functions (GFs). Each of these objects is then expanded in a perturbation series, which in practical applications is truncated to include a subset of all possible contributions. In self-consistent schemes such series are themselves expressed in terms of the exact GFs, which requires an iterative solution and makes the method intrinsically non-perturbative, effectively resumming an infinite subsets of perturbative terms. The x-body GF g x allows one to compute all x-body observables in the A-body ground state 5 . Therefore, for most applications one is mainly interested in the one-body GF.
Formally g is defined as the expectation value of a time-ordered product of annihilation and creation operators in the A-body ground state with Greek indices labelling basis states of the one-body Hilbert space H 1 . One usually works with the Fourier transform of (1) in the energy domain and recasts its perturbation series into the Dyson equation where g (0) represents some initial ansatz for g, e.g. stemming from the solution of Hartree-Fock (HF) equations. The (irreducible) self-energy Σ encodes all terms of the expansion, which is truncated in actual calculations.
Once the one-body GF is obtained as the solution of Eq. (2), the expectation value of any one-body operatorÔ 1B ≡ αβ O 1B αβ a † α a β can be computed as where C↑ denotes an integral closed on the upper imaginary plane and the one-body density matrix has been introduced.
Additionally, for a Hamiltonian with one-and two-body operators, the one-body propagator gives access to the total energy by means of the Galitski-Migdal-Koltun (GMK) sum-rule [59,60]

V. Somà
where T αβ denote the matrix elements of the one-body operator. Nowadays, realistic nuclear structure calculations require the inclusion of at least a three-body interaction in the starting Hamiltonian. In this case, the GMK sum rule needs to be generalised to [32] where the ground-state expectation value of the three-nucleon operatorŴ has to be evaluated. Such a term requires in principle the knowledge of the three-body propagator g III . This is however currently out of reach and in most of practical applications the last term in Eq. (6) is computed as i.e. by approximating the full three-body density matrix with the antisymmetrised product of the one-body one. In Ref. [34] this was shown to introduce errors smaller than 250 keV for the binding energy of oxygen isotopes.
In addition to giving access to the ground-state properties of the A-body system, the one-body GF contains information on neighbouring (A ± 1) nuclei. It becomes evident when rewriting the propagator (1) in the Lehmann representation where |Ψ A±1 i represent eigenstates of (A ± 1)-body systems while ) are one-nucleon addition and removal energies, respectively.
In Dyson GF theory the expansion starts from a particle-number conserving (e.g. a Hartree-Fock) reference state providing g (0) . On top of this, spherical symmetry is typically imposed. While such an expansion can suitably address closed-shell systems, it becomes inefficient or even breaks down as soon as pairing and/or quadrupole correlations become important. If one wishes to stick with 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 conservation 6 while maintaining spherical symmetry allows one to efficiently capture pairing correlations, thus gaining access to singly open-shell nuclei.
Dyson GFs were thus generalised to a U(1) symmetry-breaking (typically Hartree-Fock-Bogolyubov) reference state originally by Gorkov [61]. The formalism was then adapted and implemented for applications to finite nuclei in Ref. [30]. Technically, the extension is achieved by working with an A-body ground state V. Somà that is a linear combination of states with different particle numbers This leads to the definition of four one-body propagators two of which (g 11 and g 22 ) involve normal combinations of a and a † and are associated to the standard density matrix (3). The remaining two propagators (g 12 and g 21 ) invoke so-called anomalous contributions of a and a † (interpreted as the annihilation or the creation of a nucleon pair) and lead to the definition of an anomalous (or pairing) density matrix In Eqs. (10) creation operators {ā † α } define a one-body basis dual to {a † α } and are obtained viā which correspond to exchanging the state a with its time-reversal partnerã up to the phase η a [30]. The four Gorkov propagators (10) can be conveniently recast in a 2 × 2 matrix notation via Nambu's formalism [62] g αβ (t, t ) ≡ All quantities (operators, self-energy, ...) can be generalised in an analogous fashion such that one ends up with the Gorkov equation Similarly, all standard GF equations including Eqs. (3)-(8) are rewritten in a matrix form. Last but not least, a chemical potential λ needs to be introduced to guarantee that the number of particles is the correct one on average. This amounts to replacing the HamiltonianĤ with the grand potential

V. Somà
As a consequence of the symmetry breaking, observables might be contaminated by components associated to different particle numbers. Even if in practice the variance is expected to remain small 7 , the broken symmetry has to be eventually restored. While symmetry-restored formalism has been developed for other (post-Hartree-Fock-Bogolyubov) many-body methods [63,64], it remains to be formulated for Gorkov GFs.

Choice of approximation scheme
Self-consistent GF approximation schemes are defined by the content of the irreducible self-energy, which is expressed as a function of the exact GFs and encodes its perturbative expansion. There exist several ways of approximating the self-energy. The most basic one simply amounts to truncating the perturbative expansion at a certain order. More refined techniques resort to including infinite subsets of perturbationtheory terms via the definition of implicit equations. This is the case, e.g., of the so called ladder or in-medium T-matrix approximation that resums all multiple particle-particle scattering contributions 8 . This approximation scheme gained considerable attention in early ab initio applications because of its ability of tackling nucleon-nucleon interactions with strong short-range components [66,67,68,69,70,71]. These truncations and resummations are typically conveniently expressed in terms of Feynman diagrams, which facilitate the manipulation of the various terms and give an insight in their physical content.
An alternative, although not orthogonal, route was proposed in the context of quantum chemistry [72,73] and instead exploits the analytical structure of the self-energy. Similarly to the one-body GF, the exact (dynamical, i.e. energy-dependent, part of the) self-energy displays a Lehmann representation where the matrices M, N couple the single-particle motion of the nucleons (i.e. the one-body propagator in which the self-energy is inserted) to intermediate multiparticle-multihole configurations, whose energies ("bare" and resulting from the interference between them) are encoded in the matrices E > , E < and C, D (respectively). The algebraic diagrammatic construction at order n [ADC(n)] is built by demanding that, in addition to including all perturbation-theory contributions up to a given order n, the approximated self-energy has the same analytical structure as the exact one, i.e., in particular, is the same function of the energy. The latter condition requires the self-energy to contain additional sets of contributions, e.g. infinite resummations that would necessitate ad-hoc procedures are in this way automatically included in the ADC formalism. The first order, ADC(1), is simply the standard Hartree-Fock (or Hartree-Fock-Bogolyubov, in the case of Gorkov GFs) approximation. ADC(2) introduces lowest-order dynamical correlations in terms of two particle-one hole and two hole-one particle contributions. ADC(3) builds couplings between such configurations and, as a result, includes infinite-order resummations of both particle-particle/hole-hole and particle-hole ladders. Higher orders build on higher-rank particle-hole excitations in a similar but not identical fashion as in other popular many-body methods like coupled-cluster (CC) or in-medium SRG (see e.g. [74] for a connection between GF and CC formalisms). Importantly, by preserving at each order 7 Away from closed-shell systems, contributions from components with A = A are assumed to cancel out to some extent. The largest contamination is expected in differential observables across closed shells, where one of the two systems does not spontaneously break particle-number symmetry. 8 Other types of resummation are employed in other domains. For instance, in solid-state physics the resummation of particle-hole (i.e., ring) diagrams (typically in the so-called GW approximation) allows resolving the long-range features of the Coulomb force [65].

V. Somà
the analytical properties of the exact self-energy, the ADC expansion ensure that causality is not violated. Moreover, this form allows the derivation of an energy-independent auxiliary eigenvalue problem that significantly simplifies the numerical solution of Dyson or Gorkov equation, as discussed in Sec. 3.2.
The ADC scheme has been developed in the context of finite nuclei in the past few years. At present, ADC(1), ADC(2) and ADC(3) self-energies are implemented in the Dyson framework [34,36], while ADC(1) and ADC (2) are available in the Gorkov case [30,31].  is under construction within a project whose goal is to automatise the generation of the associated self-energy diagrams [75]. In Dyson theory, the ADC formalism has been recently generalised to the presence of three-body interactions [33].
Generally speaking, ADC(n) defines a truncation scheme that is systematically improvable in the sense that going to higher orders should provide results closer and closer to the exact solution, recovered in the limit of ADC(∞). Nevertheless, in the case of perturbation theory (and, by consequence, resummation methods building on MBPT), there does not exist a well-defined expansion parameter informing on the accuracy associated to a certain truncation level. Hence, there are two ways of assessing the accuracy of a given approximation: i) by comparing successive orders in the expansion and ii) via benchmarks with exact methods. Concerning the first possibility, a typical convergence behaviour for ADC(1-3) is shown in Fig. 1 for ground-state energies and root-mean-square (rms) charge radii of two representative mediummass nuclei. Two different interactions are used, the "soft" NN+3N (lnl) and the "harder" NNLO sat , see Sec. 3.4 for details. One observes a clear pattern going from ADC(1) to ADC(3) for all quantities. For total energies, while ADC(2) already yields a qualitatively good account, additional correlations introduced at  V. Somà the ADC(3) level are deemed necessary for an accurate description. In contrast, charge radii are already largely converged at the ADC(2) level. After analysing the ADC convergence behaviour, one expects ADC(4) corrections to be small for total energies and negligible for charge radii. The accurateness of the ADC(3) scheme is also confirmed by a direct benchmark against no-core shell model in 16 O, with total energies computed in the two approaches differing by less than 1% for any tested interaction [39].

Working equations
In practice, solving Dyson or Gorkov equations in their forms (2) and (14) is problematic, specially for finite systems where a solution has to be achieved for numerous (discrete) values of the energy ω. A method to overcome the problem was proposed once again in quantum chemistry [76], and consists of two steps. In the first one, after exploiting the analytical energy dependence of the propagator, Eq. (2) is transformed into an eigenvalue equation where the (energy-dependent) self-energy acts as an effective one-body potential. The second steps makes use of the analytical energy dependence of the self-energy in the form (16) and rewrites the problem as an energy-independent eigenvalue equation. The latter constitutes the working equation to be solved iteratively, and whose solutions directly provide amplitudes and energies entering Eq. (8) for the one-body propagator (see Ref. [57] for a more detailed discussion). The derivation of this energy-independent eigenvalue problem has been generalised to Gorkov theory in Ref. [30].
One disadvantage of the energy-independent formulation is that the number of energy poles, i.e. the eigenvalues and the dimensionality of the energy-independent Dyson or Gorkov matrix, increases at each self-consistent iteration. In practice, this growth is reduced via the application of Krylov projection techniques [76,77], typically implemented by means of a Lanczos algorithm (see Ref. [31] for a detailed discussion and a numerical study).

Choice of basis
Equations presented in Sec. 2 are general, i.e. are valid in any basis of choice. In an actual calculation, one needs to specify a basis in which operators, together with all relevant quantities, are expanded on. In the case of atomic nuclei one typically employs a one-body spherical harmonic oscillator (HO) basis, whose eigenfunctions are well suited to the description of a confined system 9 . A HO basis is characterised by two parameters: the oscillator inverse length Ω and the number of considered HO wave functions. The latter is usually determined by the parameter e max ≡ max(2n + l), which sets the energy threshold of a basis eigenfunction. Many-body bases are subsequently built as direct products of one-body bases. While naturally for a k-body operator one would set e kmax = k·e max , the storage of three-body matrix elements for realistic values of e max presently constitutes an issue and obliges one to work with e 3max 3 · e max . In current state-of-the-art implementations, typical values of e max = 12 − 15, e 2max = 2 · e max and e 3max = 14 − 18 are used.
The first step of a calculation consists in studying the convergence properties (of the observables of interest) with respect to the basis size and oscillator parameter. An example is showed in Fig. 2, where the basis dependence of ground-state energy and charge radius is investigated for two different nuclei and two different Hamiltonians with ADC(2) calculations. One sees that in 36 Ca an e max = 13 model space guarantees sufficiently converged results. In 68 Ni and for the higher-cutoff Hamiltonian NNLO sat (see next section for details) the convergence is not yet optimal. While for radii the crossing point of different e max 9 On the other hand, the asymptotic behaviour of HO wave functions do not correctly account for the fall-off of the nuclear wave function. As a consequence, a HO basis is not well suited to describe states near the particle continuum, where the long-range part of the wave function is particularly important. In this case, a possibility consists in complementing the HO with a basis specifically designed to account for resonances and non-resonant continuum, e.g. the Berggren basis [78,79].  Figure 2. Ground-state energies (top panels) and rms charge radii (bottom panels) of 36 Ca and 68 Ni computed with the NN+3N (lnl) and NNLO sat Hamiltonians as a function of the harmonic oscillator parameter Ω and for increasing size e max of the one-body basis. In all cases e 3max = 16 was used. The figure is taken from Ref. [39]. curves can be taken as a good estimate of the converged value, for total energies extrapolation techniques might have to be considered [80,81].

Choice of input interaction
Like most of ab initio techniques, self-consistent GFs can take in principle any nucleon-nucleon plus threenucleon (NN+3N) Hamiltonian as input. In the early 2000s, nearly all applications were performed using semi-phenomenological 10 potentials like CD-Bonn [82] or Nijmegen [83], possibly complemented with three-nucleon forces [71,84]. Although these interactions had proven successful in the description of light nuclei, their "hard" character, i.e. the associated large momentum cutoff, required the use of sophisticated resummation schemes, thus hindering applications in medium-mass nuclei. Starting from 2010, the use of Hamiltonians derived in the context of chiral effective field theory (χ-EFT) [5] began to spread. Compared to the previous phenomenological models, χ-EFT interactions present several advantages: 1. By explicitly taking into account only low-energy degrees of freedom, they have a much smaller associated cutoff; V. Somà 3. In principle, a theoretical error can be associated to a given χ-EFT Hamiltonian (which relates to the employed truncation in the EFT expansion).
Nowadays, χ-EFT interactions constitute the standard for ab initio nuclear structure calculations. The first feature, i.e. their "softness", is often further amplified by the use of similarity renormalisation group (SRG) techniques [6], i.e. unitary transformations of the Hamiltonian that further decouple low-and high-momentum modes, leading to much improved convergence properties. The third point, in contrast, has been only marginally exploited up to now (see e.g Refs. [85,86,87,88]) and will require important formal and technical developments in the near future.
Mainly three different χ-EFT Hamiltonians have been employed in recent GF calculations, all of which are discussed here. The first one, labelled NN+3N (400), is based on the next-to-next-to-next-to-leading order (N 3 LO) nucleon-nucleon potential from Entem and Machleidt [89,4] combined with the N 2 LO 3N interaction with a local regulator [90]. The 2N interaction of Ref. [89] was built with a cutoff of 500 MeV/c, however, a 400 MeV/c regulator was used for the 3N sector [91]. This Hamiltonian has been systematically applied to p-and sd-shell nuclei and yields a good reproduction of oxygen, nitrogen and fluorine binding enrgies [92,34,35]. Nevertheless, it leads to overbinding in medium-mass nuclei starting in the calcium chain and underpredicts nuclear radii even for O isotopes [93,38,52].
With the main objective of improving on the description of radii, a chiral Hamiltonian with terms up to N 2 LO was developed in Ref. [94]. 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 14 C and 16,22,24,25 O as well as charge radii of 14 C and 16 O. The resulting interaction, named NNLO sat , successfully describes the saturation of infinite nuclear matter [94] as well as various observables in mid-mass nuclei, including charge radii [52,41,42,43]. Unlike the NN+3N (400) interaction, NNLO sat employs a non-local regulator.
Motivated by the success of NNLO sat , a novel interaction named NN+3N (lnl) was presented recently [39]. The goal was to amend the original NN+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 N 3 LO potential, which provides a better description of nucleon-nucleon data compared to the lower-order NNLO sat , it yields superior features in light systems, e.g. a better reproduction of spectroscopy of natural parity states in p-and light sd-shell nuclei. NN+3N (lnl) has been also shown to provide a very good description of energy observables (ground-state energies and energy spectra) in medium-mass nuclei up to mass A ∼ 60 [39].

Ground-state properties
The total ground-state energy, or binding energy, of a nucleus constitutes the most basic nuclear structure observable. In Green's function theory, total energies are preferably computed via the generalised GMK sum-rule (6). While earlier applications made use of a 2N -only Hamiltonian, possibly complemented by a phenomenological correction to compensate for missing 3N [28,29,37], starting from 2013 calculations with realistic 2N + 3N interactions could be routinely performed. A representative example concerns the oxygen chain [34] and is shown in Fig. 3 (left). In the bottom panel, ADC(3) ground-state energies are displayed for closed-shell oxygen isotopes, computed with the NN+3N (400) interaction respectively V. Somà and dependence on SRG found in Ref. [28], we expect an accuracy of at least 5% on binding energies. All calculated values agree with the experiment within this limit. Note that the interactions employed were only constrained by 2N and 3 H and 4 He data.
the unbound d 3=2 orbit starts being filled. Other bound quasihole states are lowered resulting in additional overall binding. As a result, the inclusion of NNLO 3NFs consistently brings calculations close to the experiment and reproduces the observed dripline at 24 O [41][42][43]. Our calculations predict 25 O to be particle unbound by 1.54 MeV, larger than the experimental value of 770 keV [44] but within the estimated errors. The ground state resonance for 28 O is suggested to be unbound by 5.2 MeV with respect to 24 O. However, this estimate is likely to be affected by the presence of the continuum which is important for this nucleus but neglected in the present work.
The same mechanism affects neighboring isotopic chains. This is demonstrated in Fig. 4 for the semimagic odd-even isotopes of nitrogen and fluorine. Induced 3NF forces consistently under bind these isotopes and even predict a 27 N close in energy to 23 N. This is fully corrected by full 3NFs that strongly bind 23 N with respect to 27   excluding and including original 3N operators 11 . One notices that the addition of 3N forces is crucial for a quantitative reproduction of experimental data. In particular, when only a 2N interaction is considered, the neutron dripline is wrongly located at N = 20, while it is correctly reproduced at 24 O in the presence of 3N forces.
In such a context it can be instructive to inspect one-neutron addition and removal energies associated to dominant quasiparticle peaks, see top-left panel of Fig. 3. One sees that the 3/2 + fragment becomes bound in neutron-rich isotopes when the 2N -only Hamiltonian is employed. When 3N interactions are switched on, it is instead pushed up and remains unbound all the way to 28 O, thus explaining the position of the dripline. This observation confirmed the repulsive character of the Fujita-Miyazawa 3N interaction, as previously discussed in Ref. [95].
This result was one of the first successful applications of ab initio techniques beyond light nuclei. The oxygen chain also constituted, for a few years, a testbed where calculations from various approaches could be benchmarked, demonstrating the reliability of the different many-body truncations [92,96,23,97]. More recently, the availability of new Hamiltonians prompted calculations of heavier systems, from V. Somà calcium up to the nickel chain. An example is constituted by Gorkov GF calculations exploring the performance of three interactions along few semi-magic chains, as reported in Ref. [39]. A typical outcome is the one displayed in Fig. 3 (right), where ground-state energies of calcium isotopes as well as their differences, two-neutron separation energies, are shown as a function of neutron number. The overbinding generated by the older NN+3N (400) clearly stands out. In contrast, the newer Hamiltonians NNLO sat and NN+3N (lnl) yield an excellent reproduction of total as well as differential ground-state energies, specially once ADC(3) corrections are taken into account. A good performance is found also in the nickel chain, up to the point where current computational limitations hinder a complete model-space convergence of the calculations [39]. In addition to semi-magic chains, the theory was tested in relation to novel experimental measurements in potassium [46], titanium [47] and argon [98] chains.
The performance of these three Hamiltonians was further investigated on nuclear radii in Ref. [39]. It was found that NNLO sat provides a good account of rms charge radii all the way up to nickel. Specifically, the bulk contributions are well described already at the ADC(2) level, while finer details (e.g. the parabolic behaviour observed between 40 Ca and 48 Ca) need further improvement of the many-body truncation and/or the interaction. Density distributions provide even further insight into the way nucleons arrange themselves in the correlated nuclear medium. Starting from the point-proton distribution ρ p (r), one obtains the charge density distribution after folding with the finite proton size 12 where the sets (θ i , r i ) come from having parameterised the charge density distribution of the proton as a linear superposition of three Gaussians and have been adjusted to reproduce the proton charge form factor from electron scattering data [99]. As for radii, densities are computed directly from the one-body GF and can therefore be routinely evaluated for medium-mass systems.
An interesting example relates to the possible presence of a depletion in the central part of the charge density profile, usually referred to as bubble. One of the most likely candidates has been identified in the nucleus 34 Si [100,101,102,103,104,105,106]. In Ref. [36] this system, together with its Z + 2 partner 36 S, has been thoroughly investigated by means of both Dyson and Gorkov GF calculations. The resulting point-proton and charge density distributions are shown in Fig. 4 (left). The effect of folding with the finite size of the proton is specially visible in the centre, with an attenuation of oscillations in the charge profiles. The agreement between the computed and the measured charge distribution for the stable nucleus 36 S is excellent, which gives confidence in the prediction for the unstable 34 Si, yet unknown experimentally. For the latter, a depletion in the region below 2 fm is indeed found. Its magnitude is comparable or even larger than what previously obtained in EDF [104,105] or shell model calculations [103].
A measurement of the charge distribution of 34 Si would require electron scattering on radioactive ions. These experiments are becoming feasible only now [108], with first results on the heavier 132 Xe obtained by the SCRIT collaboration [109]. In the coming years, the case of 34 Si will thus constitute an interesting objective for electron-nucleus scattering facilities. Such a measurement would extract the electromagnetic 12 In addition, centre of mass and relativistic Darwin-Foldy corrections are taken into account by employing an effective position variable r eff i , see Ref. [36] for details.  charge form factor, related to the charge profile via

V. Somà
where q is the transferred momentum, itself related to the incident momentum p and the scattering angle θ via q = 2p sin θ/2. The calculated charge form factors for 300 MeV electron scattering on 34 Si and 36 Si are displayed in Fig. 4 (right). Clear differences appear in the angular dependence for the two systems, with a higher magnitude and a displaced position of the second minimum for 34 Si. This analysis gives indications on what range of transferred momenta, and consequently which luminosities, are necessary for identifying possible depleted density profiles in this mass region.

Excited-state properties
An asset of GF theory resides in the rich content of the one-body propagator, which does not solely provide information on the targeted (even-even) system with mass number A but also on the four neighbouring (oddeven) A ± 1 nuclei. This information is explicit in the Lehmann representation of the GF, Eq. (8). The poles E ± i of the function (8) correspond to one-nucleon addition and removal energies, as schematically depicted in Fig. 5 (top left). In addition, the associated amplitudes in the numerator represent the probabilities to reach a specific eigenstate |Ψ A+1 n (|Ψ A-1 k ) of the A + 1 (A − 1) system by adding (removing) a nucleon in a single-particle state to (from) the ground state |Ψ A 0 of the even-even system. Those amplitudes can be expanded in a single-particle basis {a † α } according to Next, spectroscopic probability matrices for the nucleon addition and removal can be built 13 , S + n ≡ U n U † n and S − k ≡ V * k V T k , respectively. Their elements read as Taking the trace over the one-body Hilbert space H 1 leads to spectroscopic factors which are the norms of the spectroscopic amplitudes. A spectroscopic factor thus sums the probabilities that an eigenstate of the A + 1 (A − 1) system can be described as a nucleon added to (removed from) a single-particle state on top of the ground state of the A-nucleon system.
The complete spectroscopic information associated with one-nucleon addition and removal processes can be collected into the spectral function S(ω), defined as the energy-dependent matrix on H 1 where the first (second) sum runs over eigenstates of H in the Hilbert space H A+1 (H A−1 ) associated with the A + 1 (A − 1) system. Taking the trace of S(ω) gives the spectral strength distribution (SDD) which is a basis-independent function of the energy. Equations (19)- (23) can be generalised to the Gorkov formalism [30].
An example of SDD computed at three different levels of approximation, ADC(1-3), is shown in Fig. 5 for the nucleus 34 Si. ADC(1), i.e. Hartree-Fock, is a mean-field (or independent-particle) approximation, which translates into a series of quasiparticle peaks with unity spectroscopic factors, i.e. there is a one-to-one correspondence between the single-particle basis {a † α } and the many-body states |Ψ A±1 i accessible via the process of adding or removing a nucleon. One can nevertheless identify the main qualitative features of the SDD, e.g. associate the ground state of 33 Si ( 35 Si) with the first peak on the left (right) of the Fermi energy characterised by spin and parity 3/2 + (7/2 − ). ADC(2) introduces the lowest-order dynamical correlations that lead to a first fragmentation of the spectral distribution. A number of fragments with small spectroscopic factors appear in the vicinity of the ADC(1) peak, which is now shifted in energy and reduced in strength, i.e. it has SF ± i < 1. ADC(3) correlations further fragment the quasiparticle strength, giving rise to a large number of small peaks and a further reduction of the main-peak spectroscopic factor. One notices that, around the Fermi energy, fragments with a good quasiparticle character, i.e. large spectroscopic 13 Here bold symbols denote matrices in the one-body Hilbert space H 1 .

V. Somà
ree Gorkovwing set of from which follows that 2 The poles of the propagators 3 are given by ω k ≡ Ω k −Ω 0 . The relation of such poles to separation energies between 2 Similarly to Eq. 5, we may equivalently write Eq. 41 asŪ k a = +U k a andV k a = −V k a . 3 As discussed later on, eigensolutions of Gorkov's equations come in pairs (ω k , −ω k ) such that one should only sum on positive solutions in Eq. 39.
) and e cent p are tic behaviour of or e cent p < 0 the s asymptotically ,  factors, survive, in accordance to Landau's Fermi liquid theory. In contrast, away from the the Fermi level the strength is spread over a wide energy interval and one can hardly identify single-particle-like excitations.
From the SDD depicted in Fig. 5 one can extract an excitation spectrum of the A − 1 (A + 1) system, by looking at increasingly negative (increasingly positive) energies on the left (right) of the Fermi surface (see top-left panel for a schematic representation). By comparing ADC(2) and ADC(3) truncations, one concludes that the latter is necessary for an accurate description of the main energy peaks [36,39], in accordance with analogous calculations in quantum chemistry [110,111,112,113]. Conservatively, one can associate uncertainty bands with ADC(2)-ADC(3) differences, as done in Fig. 6 where spectra resulting from one-nucleon addition and removal from and to 54 Ca are displayed. These four systems sit at the limits of what can be presently studied in radioactive ion beam facilities. Yet, not many experimental data are available, such that these calculations mostly represent predictions that can be presumably tested in the near future. In 53 Ca, where the ground state and two excited states are identified, NN+3N (lnl) GF calculations succeed in reproducing the measured spectrum in with good accuracy. NNLO sat instead   (3) where available. In all panels, states with E x < 5 MeV and SF < 10% are shown. The figure is adapted from Ref. [39], with the addition of new experimental data from Refs. [50,51].
mildly overestimates the splitting between the two excited states. In 53 K, the difference between the two interactions is even more striking. While NN+3N (lnl) reproduces the measured first excited state with good precision, NNLO sat predicts a wrong ordering of the 3/2 + and 1/2 + states, resulting in a ground state with incorrect spin. The potassium chain indeed constitutes an interesting case because of an unusual inversion (at N = 28) and re-inversion (at N = 32) of the spins of the ground-and first excited states [114,49]. As discussed in Refs. [51,39], Gorkov GF calculations with the NN+3N (lnl) Hamiltonian yield ax excellent reproduction of the whole trend, from 37 K to 53 K. Such studies, complementary to those focusing on ground-state observables, are not only relevant in relation to current experimental programs, but also provide a unique testing ground for the development of nuclear interactions.
The spectral representation of the one-body GF thus gives access to spectroscopic studies of odd-even nuclei. If one is interested in the excitations of an even-even system, the two-body GF (in one of its possible time orderings) has to be considered instead. A typical example is the polarisation propagator, which provides the response of the nuclear system to an external operator. In analogy to Eq. (8), its Lehmann representation reads as where n π label the excited states of the A-body system. In the numerators, the residues  14,16,22,24 O computed with the NNLO sat interaction and the SCGF many-body sing an ADC(3) self-energy. The curves are obtained by folding the discrete spectra or 16 O in (b) are from Ahrens et al. [47] (red squares) and from Ishkhanov et al. [49] eistenschneider et al. [48]. 68 Ni was rizability data are ed SCGF d around gmy and Table IV lities α D of method in tal spectra Expt.

7(3) fm 3 7(22) fm 3
for a comparison with the closest peaks in the computed discrete RPA spectrum, which is also displayed in Fig. 9. In particular, the computed strength at low energy is fragmented in two principal peaks at 10.68 MeV and 10.92 MeV, located at higher energy than the experimental PDR. For the GDR, Table IV reports the centroid calculated from the DRPA response around the main peak after the Lorentzian folding.
The α D computed by integrating the DRPA spectrum is in agreement with the experiment, also reported in Table IV. The 3.88(31) fm 3 value is obtained by including corrections from a theoretical extrapolation of the low-energy and highenergy parts of the spectrum [6], which were not accessible in the experiment of Rossi et al. [52]. Both the discrete peaks and the convoluted response in Fig. 9 confirm that the computed spectrum is somehow shifted towards higher energy as compared to the experimental excitation energies. The strength of the PDR is also underestimated.
The lack of strength in the low-energy part of the spectrum could point to insufficient constraints on the isospin-violating  36,40,48,52,54,70 Ca. Experimental data for 40 Ca in (b in (c) are from Birkhan et al. [50] (green circles) and from Ahrens et al. [47] rescaled contact terms of the NNLO sat interaction, giving a too soft symmetry energy. This is also found in calculations of the infinite nucleonic matter within the microscopic Brueckner-Hartree-Fock [53] and SCGF [54] approaches. However, the limitations scheme p interaction try energy 054327-9  [115,116,117]. The figure is taken from Ref. [41].
represent particle-hole matrix elements between excited states of the A-nucleon system. The poles appearing in the denominator instead constitute energy differences between excited states of the A-nucleon system and its ground-state. The polarisation propagator (24) is obtained as a solution of the Bethe-Salpeter equation, where Π (0) (ω) is the free polarisation propagator, and the particle-hole irreducible interaction K (ph) plays a similar role as that of the self-energy in Eq. (2) for the single-particle propagator. While the corresponding formalism has been developed and implemented for Dyson GFs, the generalisation to the Gorkov framework remains to be carried out.
In Refs. [118,41] the polarisation propagator was calculated with the aim of accessing the nuclear isovector electric dipole (E1) response. A dressed random phase approximation (DRPA) was adopted for the computation of Π(ω). The scheme makes use of a correlated, e.g. ADC(3), one-body propagator as the starting point for the RPA equations, thus going beyond a simple HF-based particle-hole resummation. Figure 7 shows results for the integrated isovector E1 photoabsorption cross section, which is directly obtained from the polarisation propagator [41]. Two representative examples are shown here, 16 O and 48 Ca. Calculations make use of the NNLO sat interaction, ensuring that nuclei have the correct size, which is a crucial property for this application. In 16 O the peak associated to the giant dipole resonance (GDR) is well reproduced. At higher excitation energies, the calculation underestimates the experimental spectrum, presumably due to missing correlations beyond the simple RPA. A similar picture emerges for 48 Ca, with the GDR peak in good agreement with recent experimental measurements [117] and the high-energy tail missing some strength. Several other closed-shell nuclei, for some of which experimental data are not yet V. Somà available, are discussed in Ref. [41]. A generalisation of this formalism to the Gorkov framework would allow us to extend these promising results to a large number of open-shell nuclei.

Lepton-nucleus scattering
The spectral function introduced in Eqs. (22)- (23) carries information about the energy-momentum distribution of the correlated nucleons. Certain scattering processes, where an external probe scatters off the nucleus, under certain kinematical conditions (e.g. characterised by a sufficiently large momentum transfer), can be described as an incoherent sum of scattering amplitudes on bound nucleons. Then, the cross section can be computed in the so-called impulse approximation and expressed as a sum of one-body terms containing a convolution with the nucleon spectral function. This is the case, e.g., of electron and neutrino scattering in the region of the quasielastic peak. Here the double differential cross section for inclusive lepton-nucleus scattering can be written as [119] where L µν is the leptonic tensor and k = (E k , k) and k = (E k , k ) are the laboratory four-momenta of the incoming and outgoing leptons, respectively. The factor C = α/(k − k ) 4 for electrons and C = G/8π 2 for neutrinos, where G = G F for neutral current (NC) and G = G F cos θ c for charged current (CC) processes. The electroweak coupling constants are α 1/137, G F = 1.1803 × 10 −5 GeV −2 [120] and cos θ c = 0.97425 [121]. The hadron tensor W µν incorporates the transition matrix elements from the target ground state |Ψ A 0 to the final states |Ψ A f due to the hadronic currents, including additional axial terms for neutrino scattering. The impulse approximation consists in factorising |Ψ A f → |p ⊗ |Ψ A−1 n , i.e. it allows to work with the outgoing nucleon of momentum p and the residual nucleus, left in a state |Ψ A−1 n . This leads to [122,123] W µν 1b (q, ω) = where ω represents the energy transfer, M N is the nucleon mass, e(p) the energy of a nucleon with momentum p. The one-body current operators j µ s depend on the spin-isospin degrees of freedom s and S h s (p, E) is the one-body spectral function normalised to the total number of nucleons. For two-body currents and hadron production, Eq. (29) extends non trivially in terms of one-and two-body spectral functions [124,125,126,123].
Following the above formalism, electron and neutrino scattering off argon and titanium isotopes was investigated in Ref. [43]. The interest in studying lepton scattering off these nuclei resides in the fact that future-generation neutrino experiments (e.g. DUNE [129]) will use liquid-argon time-projection chambers, which rely on scattering of neutrinos off 40 Ar. The nuclear component of the cross section has to be well determined for a meaningful interpretation of the measured events, in particular to reconstruct the neutrino energy with sufficient accuracy. In this respect, a tailored electron scattering experiment 14 was designed and recently performed at JLab [130,127,128]. In Ref. [43] nuclear spectral functions of stable argon, calcium and titanium isotopes were computed at the ADC(2) level with the NNLO sat interactions. A good reproduction of available charge radii and density distributions was found. Starting from these spectral functions, inclusive electron-48 Ti and electron-40 Ar cross sections were computed. They are shown in Fig. 8 (left) as a function of the energy transfer and compared to the recent experimental data from the JLab E12-14-012 collaboration [127,128]. Calculations closely follow the quasi elastic peak, thus validating the theoretical approach and the impulse approximation in particular. Next, the quasielastic neutral and charged current cross sections were studied for 1 GeV neutrino scattering. Results are displayed in the right panel of Fig. 8, also compared to scattering off 12 C. The quasielastic peak is found at a similar energy, with an increase in the magnitude of the cross section consistent with super scaling properties of inclusive reactions. The use of a Ti proton spectral function instead of an Ar neutron spectral function was also tested and found to be an excellent approximation. Further studies will be needed to thoroughly assess theoretical uncertainties.

Towards the description of nuclear reactions
While the good knowledge of electro(weak) interactions facilitate the modelling of lepton-nucleus scattering, nuclear reactions, because of the complexity of strong interactions, require a more involved V. Somà theoretical description. As a consequence, very few ab initio methods are presently capable of going beyond structure properties and directly computing reaction observables. Nevertheless, nucleon-nucleus and nucleus-nucleus reactions constitute nowadays the tools of choice to study the properties of atomic nuclei, such that progress on their ab initio description would be highly valuable.
A relatively clean process that can be used as a training field for ab initio reaction theory is quasifree nucleon knockout [131]. Quasifree scattering represents a process in which an incident nucleon (typically a proton 15 ) with an energy of few hundred MeV knocks out a bound nucleon from the isotope of interest (see Fig. 9, left panel, for a schematic illustration). Kinematical conditions are chosen such that the process can be preferentially described by a single, localised interaction between the incident and the struck nucleons, thus minimising multiple collisions for the incoming nucleon. This sudden removal mechanism suggests that the remaining A − 1 nucleons can be treated as spectators, which translates into an impulse approximation analogous to the one discussed in Sec. 4.3. As a result, the total cross section can be separated into a structure and a reaction part. The latter involves nucleons that are usually described in terms of distorted waves to account fo their propagation under the influence of the nuclear medium. Under these assumptions the differential cross section, labelled distorted-wave impulse approximation (DWIA), can be schematically written as The first factor on the right-hand side is the spectroscopic factor of the struck nucleon, encoding the structure properties of the nucleus. The second term represents an in-medium proton-nucleon cross section, determining the probability of the collision between the projectile and one of the bound nucleons. The third term contains the scattering matrices for the (effective) degrees of freedom at play, i.e. respectively i) proton and initial A-body system, ii) proton and final (A − 1)-body system, iii) struck nucleon and final (A − 1)-body system, plus the scattering wave function of the outgoing nucleon.
An ab initio calculations for some of the above quantities was recently performed within the GF formalism, in connection with a quasifree neutron knockout study on 54 Ca [50]. Specifically, spectroscopic factors and one-body overlap functions for the struck neutron were computed. The latter enter the evaluation of the scattering matrices, where they are convoluted with nucleon-nucleus optical potentials. Overlaps functions (OFs) for three different neutron states, corresponding to the first three (ground and excited) states in 53 Ca, are shown in Fig. 9 (right). GF calculations were performed at the ADC(3) level with two different interactions, NNLO sat and NN+3N (lnl). For comparison, OFs from an EDF calculation with the SLy4 Skyrme parameterisation are displayed. An overall good agreement is found between the three sets of calculations. Looking more in detail, one notices that NNLO sat better reproduces the maximum of the overlap at around 2-5 fm, i.e. the region of the nuclear surface. This clearly relates to the ability of this interaction to well describe nuclear sizes (cf. discussion in Sec. 4.1). In contrast, NN+3N (lnl) is in better agreement in the tail of the OF. This is in line with the more accurate description of the low-lying spectrum of 53 Ca, cf. Fig. 6.
In order to compute the cross section (30) for the 54 Ca(p, pn) 53 Ca reaction, in Ref. [50] the GF input was complemented by phenomenological optical potentials and in-medium nucleon-nucleon cross sections, yielding results in good agreement with shell model calculations. An analogous study had been performed 15 Since many current experiments concern unstable nuclei, they are performed in inverse kinematics, whence the use of a proton target. In the following, the case of an incident proton will be thus considered.  [132]. OFs corresponding to the main quasiparticle fragments for three angular momentum and parity channels J Π = 1/2 + , 3/2 + , 5/2 + are shown. For comparison, EDF calculations [133,134] performed with the SLy4 parameterisation are displayed. These OFs were employed in the calculation of one-proton knockout cross sections in Ref. [51].
previously for the 14 O(d, t) 13 O and 14 O(d, 3 He) 13 N reactions, leading to similar conclusions [48]. These applications might be seen as a first step towards a consistent approach to structure and reaction, and show that ab initio ingredients can be as efficient as phenomenological ones. In fact, the in-medium nucleonnucleon cross section could be already extracted from nuclear matter calculations, see e.g. Ref. [135]. Moreover, a nucleon-nucleus potential can be directly computed from the one-body self-energy, see Ref. [42] for the first applications to oxygen and calcium isotopes. Thus, in the future the full ab initio calculation of the cross section (30) can be envisaged.

Towards doubly open-shell nuclei
The development of Gorkov GF theory [30], and its subsequent implementation to finite nuclei [37,31] proved that symmetry breaking can be a powerful tool in the context of ab initio calculations. The generalisation of other many-body techniques to a symmetry-breaking scheme [92,136,13] further confirmed the validity of this strategy. Such advances have allowed to extend the reach of ab initio calculations in mid-mass systems from a few closed-shell nuclei to a large number of open-shell isotopes, e.g. to complete semi-magic isotopic or isotonic chains. In their current implementation, however, these methods do not break rotational symmetry. This results in an inefficient account of quadrupole correlations, such that the description of (doubly) open-shell systems displaying significant deformation can be problematic.
As opposed to pairing, where the strong static correlations at the Fermi surface cause the breakdown of the particle-hole expansion, in the presence of deformation one can usually produce converged calculations, i.e. compute few orders in the expansion. Nevertheless, one expects the accuracy to deteriorate with the strength of the deformation. This has been indeed observed in Gorkov GF calculations around the calcium chain, in particular for titanium and chromium isotopes characterised by mid-shell protons. For instance, by V. Somà studying neutron gaps 16 at the neutron traditional magic number N = 28 one finds an excellent agreement with experiment up to Z = 22, after which symmetry-restricted GF calculations clearly depart from data [98]. Furthermore, one can identify [137] a correlation between the deviation to experimental data and the amount of deformation (e.g. estimated by EDF calculations [138]). This situation suggests that the additional breaking of the SU(2) symmetry associated to rotational invariance could be beneficial in the extension of correlation-expansion methods to doubly open-shell nuclei. Most of existing approaches are indeed being generalised along these lines [139,140]. In the case of Gorkov GF such an extension will presumably require the use of importance truncation [141] and/or tensor factorisation techniques [142].

Towards heavy nuclei
Over the past years, GF ab initio calculations have extended their reach across the Segrè chart, going from the first application to the oxygen chain (A ∼ 20) [34] to recent computations of nickel isotopes (A ∼ 70) [39]. Such calculations rely on sophisticated numerical codes that make extensive use of available high-performance computing resources. Although the management of the computing time 17 could be problematic, the bottleneck that currently prevents (converged) calculations beyond A ∼ 100 is related to the storage of the matrix elements of 3N operators. Indeed, presently employed truncations on the three-body basis of e 3max = 14 − 18 allow keeping the size of 3N matrix elements below 100 GB, which is roughly the order of magnitude of the available memory on a single node of a state-of-the-art supercomputer. Going considerably beyond this size thus constitutes an issue. At the same time, such values of e 3max are enough to achieve reasonably converged results in the (A ∼ 60 − 70), while they become insufficient for larger isotopes [143,39].
Different strategies are being explored to overcome this issue. One possibility would be to discard beforehand, for a given e 3max , a subset of the initial 3N matrix elements. While performing a selection 18 directly on the original set might be problematic, a promising technique based on tensor factorisation algorithms has been put forward recently [141]. Since in the majority of applications 3N forces are included in the normal-ordered two-body approximation 19 , another option could consist in performing the normal ordering procedure in a different (smaller) basis than the HO one, e.g. in momentum space.
Even with current limitations, however, some meaningful results can be produced for nuclei above mass A = 100. Indeed, as discussed in Sec. 3.3, not all observables show the same convergence pattern. In particular, while ground-state energy curves get lower and lower as the basis size is increased, radii tilt around a fixed point that can be assumed to correspond to the infinite-basis result. This allow to provide an estimate of the radius (and, correspondingly, of the density distribution) in bases for which the energy is far from being converged. Based on this observation, charge densities of closed-and open-shell tin and xenon isotopes have been recently calculated within GF theory [40]. Examples are reported in Fig. 10, where the NNLO sat interaction has been employed. In the left panel, the charge distribution of several nuclei is shown, exemplifying the typical range of system that can be presently accessed. In the right panel, the charge density of 132 Xe is displayed and compared to a two-point Fermi distribution fitted on the recent experimental data from the SCRIT collaboration [109]. The two are in very good agreement at the 16 Neutron gaps, defined as differences of two-neutron separation energies, are one of possible observables that meaningfully estimate the "magic" character of a neutron number. 17 In this respect, while the numerical cost of Dyson GF calculations grows with the mass number, the one of Gorkov GFs solely depend on the basis dimension [31]. 18 I.e., using a different truncation than the e 3max truncation. 19 It consists of two steps: i) a normal ordering of the (3N ) Hamiltonian with respect to a reference state and ii) the disregard of operators of rank higher than two (see Ref. [144] for a pedagogical description and the generalisation to the case of symmetry-breaking reference states). terval of e ective momentum transfers between 0.8 fm ≠1 and 1.1 fm ≠1 being slightly o the error bars. To discard the density oscillations within the nucleus as the source of the discrepancy, we fitted a two-point Fermi density to the radius and surface predicted by the theory. Calculations using this Fermi distribution gave results within the band obtained from the genuine SCGF density. This confirms the inability of the experiment to give insights on the internal structure of the nucleus without going past the second minimum in the cross-section. As a comparison, the results obtained with the recently proposed NN + 3N (lnl) interaction [34] are displayed as well, scaled upwards for clarity. Contrary to NNLOsat, it fails at reproducing the experimental values, as expected with an underestimated charge radius. This demonstrates the unique capacity of NNLOsat to reproduce radii and density distributions, and sets an important precedent in the use of SCGF with the NNLOsat interaction for pre-or post-diction of experimental results from electron scattering o exotic nuclei. In particular, this motivates experimental measurements at higher momentum transfer to properly gauge the internal structure of nuclei.
Having proved the capacity of SCGF and NNLOsat to give meaningful insights on the charge radius and density distributions of 132 Xe, charge densities have been calculated for 100 Sn, 132 Sn, 136 Xe and 138 Xe as well and are displayed in Fig. 4. The behaviour of the charge distributions is qualitatively similar for all of them, with oscillations in the density within the nucleus and the possibility of a light depletion at its center.
The charge radii extracted from our calculations are displayed for di erent Sn and Xe isotopes in Tab. I and compared with experimental results [56]. Our results show overall a good reproduction of the experimental data and are a proof of the capacity of NNLOsat to produce accurate results in the heavy nuclei regime, even despite the inability to obtain converged values for the groundstate energy. In the future, more accurate calculations with smaller errors may uncover slight di erences between NNLOsat and the experimental values. Among the nuclei studied, 100 Sn stands out as a particularly interesting case. Sitting close to the proton dripline [57], at the end of super-allowed --decay chains [58,59] and with the largest strength known in alloweddecay [60], and being expected to be the heaviest doubly-magic nucleus with N = Z [61], experimental data in its area are scarce [62].
In particular, neither its spectrum nor its radius have been measured yet. While its spectrum has recently been predicted from first principles [30], Tab. I displays the first ab initio prediction of its charge radius. Conclusion. Our calculations demonstrated the capacity of SCGF and the NNLOsat interaction to give a meaningful estimation of the charge radius and charge density distribution of heavy nuclei up to mass A = 138 which had never been studied before. We computed successfully the charge radius and density distribution of 132 Sn, 132 Xe, 136 Xe, and 138 Xe, mostly agreeing with known experimental values, and gave the first ab initio prediction for the charge radius and density distribution of 100 Sn. In particular, we reproduced the experimental cross-section of the SCRIT electron scattering experiment for 132 Xe, demonstrating the capacity of ab initio methods with well-designed chiral interactions to be used for the internal structure study of heavy exotic nuclei, alongside new experimental facilities. Our errors bars, though conservative, are small enough to shed light on discrepancies with experimental values, informing theory and putting constraints on experiments. In particular, our results are a motivation for measurements at higher momentum transfer to probe the internal structure of the nuclei.
The authors are grateful to Petr Navrátil for providing matrix elements of the NNLOsat interaction and to him and Carlotta Giusti for several useful discussions. This work is supported by the UK Science and Technology  at the ADC(2) level. As such, we do not discuss di erences between ADC(2) and ADC (3) results any further in this Letter. In the following, we will hence represent our results as a band obtained for frequencies from 10 to 14 MeV at N max = 13 and 12 to 14 MeV at N max = 11, for E 3max = 16. From this procedure, the charge radius of 132 Xe is estimated to be 4.824 ± 0.124 fm, which agrees with the value extracted from the SCRIT experiment recently, namely Èr 2 Í 1/2 = 4.79 +0.11 ≠0.08 fm [10]. For comparison, the calculations have been reproduced using the newlyproposed NN + 3N (lnl) interaction [34], which is known to have good convergence properties with respect to the model space size and to give results similar to the very succesful 1.8/2.0(EM) interaction [31]. In contrast with NNLOsat, the charge radius obtained for 132 Xe is 4.070 ± 0.045 fm, largely underestimating the experimental value consistently with studies on lighter nuclei [34]. Despite this failure at reproducing the experimental value, one notices that NN + 3N (lnl) yield better-converged values than NNLOsat as expected.
Additionally to the sole charge radius, another quantity that can be computed from SCGF calculations is the charge density distribution. In the case of 132 Xe, the SCRIT group extracted the parameters c and t for a two-parameter Fermi charge distribution fl(r) = fl 0 / {1 + exp[4 ln 3(r ≠ c)/t]}. Fig. 2 displays this twopoint Fermi distribution as a dotted line with a gray band representing the error bars, while the green band represents our SCGF calculations. It can be observed that while the SCGF calculations agree with the 2-point Fermi distribution at the surface of the nucleus, though slightly over-predicting the charge radius, we obtain an oscillating behaviour for the density inside the nucleus  Figure 10. Left: Charge density distributions of selected Sn and Xe isotopes computed with the NNLO sat interaction in Gorkov GF theory at the ADC(2) level. Each curve is shifted upwards with respect to the one below by 0.025 fm −3 for clarity purposes. Coloured band account for uncertainties stemming from model-space convergence. Right: Charge density distributions of 132 Xe computed with the NNLO sat interaction in Gorkov GF theory at the ADC(2) level. The dotted line with grey band corresponds to the two-point Fermi distribution with parameter and error bars extracted from Ref. [109]. Figures are taken from Ref. [40].
surface and in the tail of the distribution. In the interior, the two-point Fermi distribution and the lack of high-momentum transfer data lead to flat behaviour for the experimental distribution, whereas the computed density shows a well defined oscillation pattern. This example shows that even in present implementations of GF calculations it is possible to provide relevant predictions in the region A = 100 − 150.