The single-particle propagator describes the probability amplitude for adding (removing) a particle in state α at one time to the ground state and propagating on top of that state until a later time when it is removed (added) in state β [24]. In addition to the conserved orbital and total angular momentum ( ℓ and j , respectively), the labels α and β in Equation 1 refer to a suitably chosen single-particle basis. We employed a coordinate-space basis in our original ⁴⁸Ca calculation in Ref. [25] but have since updated to using a Lagrange basis [26] in all subsequent calculations (including that of ²⁰⁸Pb from Ref. [27]). It is convenient to work with the Fourier-transformed propagator in the energy domain, G ℓ j α , β ; E = 〈 Ψ 0 A | a α ℓ j 1 E − H ̂ − E 0 A + i η a β ℓ j † | Ψ 0 A 〉 + 〈 Ψ 0 A | a β ℓ j † 1 E − E 0 A − H ̂ − i η a α ℓ j | Ψ 0 A 〉 , (1) with E 0 A representing the energy of the non-degenerate ground state | Ψ 0 A 〉 . Many interactions can occur between the addition and removal of the particle (or vice versa), all of which need to be considered to calculate the propagator. No assumptions about the detailed form of the Hamiltonian H ̂ need be made for the present discussion, but it will be assumed that a meaningful Hamiltonian exists that contains two-body and three-body contributions. The application of the perturbation theory then leads to the Dyson equation [24], which is given by G ℓ j α , β ; E = G ℓ 0 α , β ; E + ∑ γ , δ G ℓ 0 α , γ ; E Σ ℓ j * γ , δ ; E G ℓ j δ , β ; E , (2) where G ℓ ( 0 ) ( α , β ; E ) corresponds to the unperturbed propagator (the propagator derived from the unperturbed Hamiltonian, H 0 , which in the DOM corresponds to the kinetic energy) and Σ ℓ j * ( γ , δ ; E ) is the irreducible self-energy [24]. The hole spectral density for energies below ε F is obtained from S ℓ j h α , β ; E = 1 π I m G ℓ j α , β ; E , (3) where the h superscript signifies it is the hole spectral amplitude. For brevity, we drop this superscript for the rest of this review. The diagonal element of Equation 3 is known as the (hole) spectral function identifying the probability density for the removal of a single-particle state with quantum numbers α ℓ j at energy E . The single-particle density distribution can be calculated from the hole spectral function in the following way, ρ ℓ j p , n r = ∑ ℓ j 2 j + 1 ∫ − ∞ ε F d E S ℓ j p , n r , r ; E , (4) where the ( p , n ) superscript refers to protons or neutrons, and ε F = 1 2 ( E 0 A + 1 − E 0 A − 1 ) is the average Fermi energy which separates the particle and hole domains [24]. The number of protons and neutrons ( Z , N ) is calculated by integrating ρ ℓ j ( p , n ) ( r ) over all space. In addition to the particle number, the total binding energy can be calculated from the hole spectral function using the Migdal–Galitski sum rule [24], E 0 N , Z = 1 2 ∑ α β ∫ 0 ε F d E ⟨ α | T ̂ | β ⟩ S h α , β ; E + δ α β E S h α , α ; E . (5)

This expression assumes that the dominant contribution involves the two-nucleon interaction [28, 29] and the ℓ j labels have been subsumed in α and β .

To visualize the spectral function of Equation 3, it is useful to sum (or integrate) over the basis variables, α , so that only energy dependence remains, S ℓ j ( E ) . The spectral strength S ℓ j ( E ) is the contribution at energy E to the occupation from all orbitals with angular momentum ℓ j . It reveals that the strength for a shell can be fragmented, rather than being isolated at the independent particle model (IPM) energy levels. Figure 1 shows the spectral strength for a representative set of neutron shells in ²⁰⁸Pb that would be considered bound and fully occupied in the IPM. The location of the peaks in Figure 1 corresponds to the energies of discrete bound states with one nucleon removed. For example, the s 1 / 2 spectral function in Figure 1 has four peaks, three below ε F corresponding to the 0s 1 / 2 , 1s 1 / 2 , and 2s 1 / 2 quasihole states, and one above ε F corresponding to the 3s 1 / 2 quasiparticle state. The quasihole wave functions of these bound states can be obtained by transforming the Dyson equation into a nonlocal Schrödinger-like equation by disregarding the imaginary part of Σ * ( α , β ; E ) , ∑ γ 〈 α | T ℓ + R e Σ ℓ j * ε ℓ j n | γ 〉 ψ ℓ j n γ = ε ℓ j n ψ ℓ j n α , (6) where 〈 α | T ℓ | γ 〉 is the kinetic energy matrix element, including the centrifugal term. The wave function, ψ ℓ j n ( α ) , is the overlap between the A and A − 1 systems and the corresponding energy, ε ℓ j n , is the energy required to remove a nucleon with the particular quantum numbers n ℓ j , ψ ℓ j n α = 〈 Ψ n A − 1 | a α ℓ j | Ψ 0 A 〉 , ε ℓ j n = E 0 A − E n A − 1 . (7)

When solutions to Equation 6 are found near the Fermi energy where there is naturally no imaginary part of the self-energy, the normalization of the quasihole is well defined as the spectroscopic factor, Z ℓ j n = 1 − ∂ Σ ℓ j * α q h , α q h ; E ∂ E ε ℓ j n − 1 , (8) where α q h corresponds to the quasihole state in Equation 7. The quasihole peaks in Figure 1 get narrower as the levels approach ε F , which is a consequence of the imaginary part of the irreducible self-energy decreasing when approaching ε F . In fact, the last mostly occupied neutron level in Figure 1 (2p 1 / 2 ) has a spectral function that is essentially a delta function peaked at its energy level, where the imaginary part of the self-energy vanishes. For these orbitals, the strength of the spectral function at the peak corresponds to the spectroscopic factor in Equation 8. The spectroscopic factor can be probed using the exclusive ( e , e ′ p ) reaction which will be discussed in Section 2.4 (see also Refs. [31, 32]).

The Dyson equation, Equation 2, simplifies the complicated task of calculating G ( α , β ; E ) from Equation 1 to finding and inverting a suitable Σ * ( α , β ; E ) (suppressing the ℓ j labels). It was recognized long ago that Σ * ( α , β ; E ) represents the potential that describes elastic-scattering observables [33]. The link with the potential at negative energy is then provided by the Green’s function framework as was realized by Mahaux and Sartor who introduced the DOM as reviewed in Ref. [21]. The analytic structure of the nucleon self-energy allows one to apply the dispersion relation, which relates the real part of the self-energy at a given energy to a dispersion integral of its imaginary part over all energies. The energy-independent correlated Hartree–Fock (HF) contribution [24] is removed by employing a subtracted dispersion relation with the Fermi energy used as the subtraction point [21]. The subtracted form has the further advantage that the emphasis is placed on energies closer to the Fermi energy for which more experimental data are available. The real part of the self-energy at the Fermi energy is then still referred to as the HF term and is sufficiently attractive to bind the relevant levels at about the correct energies. In practice, the imaginary part is assumed to extend to the Fermi energy on both sides while being very small in its vicinity. The subtracted form of the dispersion relation employed in this work is given by R e Σ * α , β ; E = R e Σ * α , β ; ε F − P ∫ ε F ∞ d E ′ π I m Σ * α , β ; E ′ 1 E − E ′ − 1 ε F − E ′ + P ∫ − ∞ ε F d E ′ π I m Σ * α , β ; E ′ 1 E − E ′ − 1 ε F − E ′ , (9) where P is the principal value. The static term, R e Σ * ( α , β ; ε F ) , is denoted by Σ H F from here on. Equation 9 constrains the real part of Σ * ( α , β ; E ) by empirical information of its HF and imaginary parts which are closely tied to experimental data. Initially, standard functional forms for these terms were introduced by Mahaux and Sartor who also cast the DOM potential in a local form by a standard transformation which turns a nonlocal static HF potential into an energy-dependent local potential [34]. Such an analysis was extended in Refs. [35, 36] to a sequence of Ca isotopes and in Ref. [37] to semi-closed-shell nuclei heavier than Ca. The transformation to the exclusive use of local potentials precludes a proper calculation of the nucleon particle number and expectation values of the one-body operators, like the charge density in the ground state (see Equation 4). This obstacle was eliminated in Ref. [38], but it was shown that the introduction of nonlocality in the imaginary part was still necessary in order to accurately account for particle number and the charge density [22]. Theoretical work provided further support for this introduction of a nonlocal representation of the imaginary part of the self-energy [39, 40]. A review detailing these developments was published in Ref. [23].

To use the DOM self-energy for predictions, the parameters of the self-energy are constrained through weighted χ 2 minimization (using the Powell method [44]) by measurements of elastic differential cross sections ( d σ d Ω ) , analyzing powers ( A θ ) , reaction cross sections ( σ react ) , total cross sections ( σ tot ) , charge density ( ρ ch ) , energy levels ( ε n ℓ j ) , particle number, and the root-mean-square charge radius ( R c h ) . The angular dependence of Σ ( r , r ′ ; E ) is represented in a partial-wave basis, and the radial component is represented in a Lagrange basis using Legendre and Laguerre polynomials for scattering and bound states, respectively. The bound states are found by diagonalizing the Hamiltonian in Equation 6 and the propagator is found by inverting the Dyson equation, Equation 2, whereas all scattering calculations are performed in the framework of R -matrix theory [26]. Whereas it has been suggested in Refs. [45–47] that charge-exchange reactions to isobaric analog states could further constrain the isovector potential, charge-exchange data were not included in the fits reviewed in this article. Reasonable cross sections are obtained with our DOM potential, suggesting that these data, although important, are not sufficient to alter the conclusions of our work significantly. This may be due to the use of nonlocal potentials as opposed to the local ones used in Refs. [45, 46] based on Ref. [48].

When constraining the ⁴⁸Ca self-energy, the isoscalar part is largely determined by the nearby N = Z ⁴⁸Ca nucleus. Therefore, using our ⁴⁰Ca parametrization from Ref. [31] as a starting point, we only needed to fit the asymmetric parameters of the ⁴⁸Ca potential [25, 30]. This resulted in a ⁴⁸Ca self-energy that closely reproduced all training data [25]. In the case of ²⁰⁸Pb, there is not a nearby nucleus with N = Z ; therefore, we started from the ⁴⁸Ca parameters of Ref. [32] and varied both the isoscalar and isovector parameters to reproduce experimental data. To illustrate how well this method works, we show the result of the ²⁰⁸Pb fit below.

Proton reaction cross sections together with the DOM result are displayed in panel (a) of Figure 2. The neutron total cross section is shown in panel (b) of Figure 2. Both aggregate cross sections play an important role in determining volume integrals of the imaginary part of the self-energy, thereby providing strong constraints on the depletion of IPM orbits. The elastic differential cross sections of proton and neutrons up to 200 MeV are shown in panel (a) of Figure 3. Panel (b) contains the analyzing powers for neutrons and protons which strongly constrain the spin-orbit components of the self-energy.

The charge density of ²⁰⁸Pb is shown in panel (a) of Figure 4. The experimental band is extracted from elastic electron scattering differential cross sections [49]. This dataset is well reproduced after using the DOM charge density from Figure 4 as the ingredient in a relativistic elastic electron scattering code [52]. The corresponding elastic electron scattering cross section is shown in panel (b) of Figure 4 and compared to experiment with all available data transformed to electron energy of 502 MeV in the center-of-mass frame [51].

In Figure 5, single-particle levels calculated using Equation 6 are compared to the experimental values for protons and neutrons in panels (a) and (b), respectively. The middle column consists of levels calculated using the full DOM and the right column contains the experimental levels. The first column of the figures represents a calculation using only the static part of the self-energy, corresponding to the Hartree–Fock (mean-field) contribution. It is clear from these level diagrams that the mean-field overestimates the particle-hole gap (see also Ref. [53]). The inclusion of the dynamic part of the self-energy is necessary to reduce this gap and properly describe the energy levels [21]. Furthermore, the effect of including the dynamic part of the self-energy on the proton levels is stronger than the effect on the neutron levels. This suggests that protons deviate more from the IPM than neutrons in ²⁰⁸Pb.

The number of neutrons and protons in the DOM fit of ²⁰⁸Pb, calculated by integrating Equation 4 using shells up to ℓ ≤ 20 , is shown in Table 1. As there are 82 protons and 126 neutrons in ²⁰⁸Pb, the reported values are accurate to within a fraction of a percent. The binding energy of ²⁰⁸Pb was fit to the experimental value using Equation 5. As there is no way at present to assess the contribution of three-body interactions to the ground-state energy, we employ the present approximation which applies when only two-body interactions occur in the Hamiltonian, to ensure that enough spectral strength occurs at negative energy which has implications for the presence of high-momentum components. Also shown in Table 1 is R ch 208 calculated as the RMS radius of the charge density displayed in Figure 4.

The reproduction of all available experimental data indicates that we have realistic self-energies of ²⁰⁸Pb and similarly for ⁴⁸Ca [25, 32] capable of describing both bound-state and scattering processes. With these self-energies, we can therefore make predictions of observables such as the neutron skin. Additionally, a parallel DOM analysis of these and other nuclei was conducted using Markov Chain Monte Carlo (MCMC) to optimize the potential parameters employing the same experimental data and a very similar functional form but with a reduced number of parameters. All observables from this MCMC fit fell within one standard deviation of those presented above [55, 56].

Spectroscopic factors come directly from the self-energy through Equation 8, making the DOM ideal for predicting ( e , e ′ p ) cross sections (see the ⁴⁰Ca ( e , e ′ p ) 39 K analysis in Ref. [31]). When we tried to calculate ⁴⁸Ca ( e , e ′ p ) 47 K using the fit from Ref. [25], we found that the spectroscopic factors were too large to describe the data. Unlike in ⁴⁰Ca, there is a lack of high-energy ( E > 100 MeV) proton reaction cross-sectional data in ⁴⁸Ca. This allowed the fit of Ref. [25] to predict proton reaction cross sections which fell off for higher energies. Consequentially, the ⁴⁸Ca proton spectroscopic factors were too large to describe ⁴⁸Ca ( e , e ′ p ) 47 K data to the same degree of accuracy achieved for ⁴⁰Ca [31]. Observing for E c.m. > 150 MeV that σ react ( E ) is close to constant, we used the ratio of σ react ( E ) measurements of ⁴⁰Ca and ⁴⁸Ca at 700 MeV [57] to scale the ⁴⁰Ca σ react ( E ) data such that it could be used as a constraint for ⁴⁸Ca. Thanks to the dispersion relation, Equation 9, the increased I m Σ ( r , r ′ ; E ) to accommodate higher reaction cross sections at positive energies pulls strength from below ε F . This reduced the spectroscopic factors which then allowed for accurate descriptions of ⁴⁸Ca ( e , e ′ p ) 47 K cross sections [32]. This only altered the proton parameters; thus, the neutron skin remained unchanged at R skin = 0.25 fm. This demonstrates that once a sufficiently complete set of data is used, the DOM is capable of making accurate predictions.

The valence spectroscopic factors in ²⁰⁸Pb are consistent with the observations of Ref. [58] and the interpretation of Ref. [59]. The past extraction of spectroscopic factors using the ( e , e ′ p ) reaction yielded a value around 0.65 for the valence 2 s 1 / 2 orbit [60] based on the results of Refs. [61, 62]. Although the use of nonlocal optical potentials may slightly increase this value as shown in Ref. [31], it may be concluded that the value of 0.69 obtained from the DOM analysis is consistent with the past result. Nikhef data obtained in a large missing energy and momentum domain [63] can now be consistently analyzed employing the complete DOM spectral functions.

Correlations can also be studied through the momentum distribution, n ( k ) , which represents the diagonal of the double Fourier transform of the single-particle density matrix. The calculated DOM momentum distributions of ⁴⁸Ca and ²⁰⁸Pb are shown in Figure 6. The high-momentum tail of n ( k ) arises from short-range correlations (SRCs), which is another manifestation of many-body correlations beyond the IPM description of the nucleus [64]. This high-momentum content can be quantified by integrating the momentum distribution above the Fermi momentum. Using k F = 270 MeV/c, 13.4% of protons and 10.7% of neutrons have momenta greater than k F in ²⁰⁸Pb, whereas ⁴⁸Ca has 14.6% high- k protons and 12.6% high- k neutrons. These numbers are in qualitative agreement with what is observed in the high-momentum knockout experiments conducted by the CLAS collaboration at Jefferson Lab [65]. Furthermore, the fraction of high-momentum protons is larger than the fraction of high-momentum neutrons. These features were predicted by ab initio calculations of asymmetric nuclear matter reported in Refs. [66–68] which demonstrated unambiguously that the inclusion of the nucleon–nucleon tensor force, when constrained by nucleon–nucleon scattering data, is responsible for making protons more correlated with increasing nucleon asymmetry at normal density. This supports the n p -dominance picture in which the dominant contribution to SRC pairs comes from n p SRC pairs which arise from the tensor force in the nucleon–nucleon interaction [69, 70]. Due to the neutron excess in ²⁰⁸Pb and ⁴⁸Ca, there are more neutrons available to make n p SRC pairs which lead to an increase in the fraction of high-momentum protons.

In the DOM, this high-momentum content is determined by how much strength exists in the hole spectral function at large, negative energies. The hole spectral function is constrained in the fit by the particle number, binding energy, and charge density. Whereas the particle number and charge density can only constrain the total strength of the hole spectral function, the binding energy constrains how the strength of the spectral function is distributed in energy. This arises from the energy-weighted integral in Equation 5, which will push strength of the spectral function to more negative energies in order to achieve more binding. This, in turn, alters the momentum distribution, thus partially constraining the high-momentum content. It should be noted that the DOM does not exhibit the characteristic energy dependence of high-momentum strength distributions [71] as reported in Ref. [22]. Such a dependence is more difficult to be implemented as it requires abandoning the factorization of spatial and energy dependence of the DOM self-energy (see Equation 14).