Toward an accurate strongly-coupled many-body theory within the equation of motion framework

TL;DR

This paper develops an advanced non-perturbative many-body nuclear response theory within the equation of motion framework, incorporating complex configurations like two quasiparticles coupled to two phonons to improve spectral predictions.

Contribution

It introduces a novel approach that explicitly includes beyond-mean-field correlations up to six fermions, enhancing the accuracy of nuclear response calculations.

Findings

Higher configurations improve spectral detail accuracy.

The method successfully describes both high-energy and low-energy nuclear responses.

Results align well with experimental data for selected nuclei.

Abstract

Non-perturbative aspects of the quantum many-body problem are revisited, discussed and advanced in the equation of motion framework. We compare the approach to the two-fermion response function truncated on the two-body level by the cluster expansion of the dynamical interaction kernel to the approach known as time blocking approximation. Such a comparison leads to an extended many-body theory with non-perturbative treatment of high-order configurations. The present implementation of the advanced theory introduces a new class of solutions for the response functions, which include explicitly beyond-mean-field correlations between up to six fermions. The novel approach, which includes configurations with two quasiparticles coupled to two phonons (2q$\otimes$2phonon), is discussed in detail for the particle-hole nuclear response and applied to medium-mass nuclei. The proposed developments are implemented numerically on the basis of the relativistic effective meson-nucleon Lagrangian and compared to the models confined by two-fermion and four-fermion configurations, which are considered as state-of-the-art for the response theory in nuclear structure calculations. The results obtained for the dipole response of $^{42,48}$Ca and $^{68}$Ni nuclei in comparison to available experimental data show that the higher configurations are necessary for a successful description of both gross and fine details of the spectra in both high-energy and low-energy sectors.

Tables

Table 1. Table 1: Table of symbols

Symbol	Name	Equations and Figures
$G_{{11}^{\prime }}\left(t-t^{\prime }\right)$	One-fermion in-medium propagator	Eqs. (1,8)
$G_{12,1^{\prime }{2}^{\prime }}\left(t-t^{\prime }\right),G_{12,1^{\prime }{2}^{\prime }}^{pp}\left(\omega \right)$	Two-time particle-particle in-medium propagator	Eqs. (11,15)
$H$	Many-body Hamiltonian	Eq. (4)
$v_{1234},{\overline{v}}_{1234}$	Non-symmetrized and antisymmetrized two-fermion bare interaction	Eq.(6)
$R_{12,1^{\prime }{2}^{\prime }}\left(t-t^{\prime }\right),R_{12,1^{\prime }{2}^{\prime }}\left(\omega \right)$	Two-time particle-hole in-medium propagator (response function)	Eqs. (12,13)
$T_{{11}^{\prime }}\left(t-t^{\prime }\right)$	One-fermion T-matrix	Eq. (26)
${\Sigma }_{{11}^{\prime }}\left(t-t^{\prime }\right),{\Sigma }_{{11}^{\prime }}\left(\omega \right)$	Irreducible one-fermion self-energy	Eq. (30)
${\Sigma }_{{11}^{\prime }}^{\left(0\right)}$	Static part of the one-fermion self-energy	Eq. (35)
${\Sigma }_{{11}^{\prime }}^{\left(r\right)}\left(t-t^{\prime }\right),{\Sigma }_{{11}^{\prime }}^{\left(r\right)}\left(\omega \right)$	Dynamical part of the one-fermion self-energy	Eqs. (46,50), Figs. 1,3
${𝒩}_{{121}^{\prime }{2}^{\prime }}$	Two-fermion norm kernel	Eq. (65)
$R_{12,1^{\prime }{2}^{\prime }}^{\left(0\right)}\left(\omega \right)$	Uncorrelated particle-hole propagator	Eq. (68)
$T_{12,1^{\prime }{2}^{\prime }}\left(t-t^{\prime }\right)$	Two-time particle-hole T-matrix	Eq. (LABEL:F1)
$K_{12,1^{\prime }{2}^{\prime }}\left(t-t^{\prime }\right)$	Irreducible two-time particle-hole interaction kernel	Eqs. (72,73)
$K_{12,1^{\prime }{2}^{\prime }}^{\left(0\right)}$	Static part of the particle-hole interaction kernel	Eqs. (74,80), Fig. 4
$K_{12,1^{\prime }{2}^{\prime }}^{\left(r\right)}\left(t-t^{\prime }\right)$	Dynamical part of the particle-hole interaction kernel	Eqs. (74,LABEL:Fr-83), Figs. 5-8
${\tilde{V}}_{1234}$	Static effective two-fermion interaction	Eq. (96)
$V^{\left(e\right)}$(12,34)	Energy-dependent effective two-fermion interaction	Eq. (96)
${\tilde{\Sigma }}_{{11}^{\prime }}$	Static part of the one-fermion self-energy in PVC approaches	Eq. (97)
${\Sigma }_{{11}^{\prime }}^{\left(e\right)}\left(\omega \right)$	Dynamical part of the one-fermion self-energy in PVC approaches	Eq. (97)
${\Phi }_{12,1^{\prime }{2}^{\prime }}\left(\omega \right)$	Dynamical particle-hole interaction kernel in PVC-TBA	Eq. (119), Fig. 10
${\Phi }_{12,1^{\prime }{2}^{\prime }}^{\left(n\right)}\left(\omega \right)$	Dynamical particle-hole interaction kernel in EOM/R(Q)TBA(n)	Eqs. (124,129), Fig. 13

Equations

$$G(1,1^{\prime })\equiv G_{11^{\prime }}(t-t^{\prime })=-i⟨T\psi (1){\psi }^{†}(1^{\prime })⟩,$$

$$\psi (1)=e^{iH{t}_1}{\psi }_1{e}^{-iH{t}_1},{\psi }^{†}(1)=e^{iH{t}_1}{{\psi }^{†}}_1{e}^{-iH{t}_1},$$

$$[{\psi }_1,{{\psi }^{†}}_{1^{\prime }}{]}_{+}\equiv {\psi }_1{{\psi }^{†}}_{1^{\prime }}+{{\psi }^{†}}_{1^{\prime }}{\psi }_1={\delta }_{11^{\prime }},$$

$${\left[{\psi }_1,{\psi }_{1^{\prime }}\right]}_{+}\equiv {\psi }_1{\psi }_{1^{\prime }}+{\psi }_{1^{\prime }}{\psi }_1=0,$$

$${\left[{{\psi }^{†}}_1,{{\psi }^{†}}_{1^{\prime }}\right]}_{+}\equiv {{\psi }^{†}}_1{{\psi }^{†}}_{1^{\prime }}+{{\psi }^{†}}_{1^{\prime }}{{\psi }^{†}}_1=0.$$

$$H=H^{(1)}+V^{(2)}+W^{(3)}+...$$

$$H^{(1)}=\sum _{12}{t}_{12}{\psi }_1^{†}{\psi }_2+\sum _{12}{v}_{12}^{(MF)}{\psi }_1^{†}{\psi }_2\equiv \sum _{12}{h}_{12}{\psi }_1^{†}{\psi }_2,$$

$$V^{(2)}=\frac{1}{4}\sum _{1234}{\bar{v}}_{1234}{{\psi }^{†}}_1{{\psi }^{†}}_2{\psi }_4{\psi }_3,$$

$$W^{(3)}=\frac{1}{36}\sum _{123456}{\bar{w}}_{123456}{{\psi }^{†}}_1{{\psi }^{†}}_2{{\psi }^{†}}_3{\psi }_6{\psi }_5{\psi }_4$$

$$G_{11^{\prime }}(\epsilon )=\sum _n\frac{{\eta }_1^n{\eta }_{1^{\prime }}^{n\ast }}{\epsilon -(E_n^{(N+1)}-E_0^{(N)})+i\delta }+$$

$$+\sum _m\frac{{\chi }_1^m{\chi }_{1^{\prime }}^{m\ast }}{\epsilon +(E_m^{(N-1)}-E_0^{(N)})-i\delta }.$$

$${\eta }_1^n=⟨0^{(N)}|{\psi }_1|n^{(N+1)}⟩,{\chi }_1^m=⟨m^{(N-1)}|{\psi }_1|0^{(N)}⟩.$$

$$G$$

$$=$$

$$R(12,1^{\prime }{2}^{\prime })$$

$$R_{12,1^{\prime}2^{\prime}}(\omega)=\sum\limits_{\nu>0}\Bigl{[}\frac{\rho^{\nu}_{21}\rho^{\nu\ast}_{2^{\prime}1^{\prime}}}{\omega-\omega_{\nu}+i\delta}-\frac{\rho^{\nu\ast}_{12}\rho^{\nu}_{1^{\prime}2^{\prime}}}{\omega+\omega_{\nu}-i\delta}\Bigr{]}$$

$${\rho }_{12}^{\nu }=⟨0|{\psi }_2^{†}{\psi }_1|\nu ⟩$$

$$i{G}_{12,1^{\prime }{2}^{\prime }}^{pp}(\omega )=\sum _{\mu }\frac{{\alpha }_{21}^{\mu }{\alpha }_{2^{\prime }{1}^{\prime }}^{\mu \ast }}{\omega -{\omega }_{\mu }^{(++)}+i\delta }-\sum _{\varkappa }\frac{{\beta }_{12}^{\varkappa \ast }{\beta }_{1^{\prime }{2}^{\prime }}^{\varkappa }}{\omega +{\omega }_{\varkappa }^{(--)}-i\delta },$$

$${\alpha }_{12}^{\mu }=⟨0^{(N)}|{\psi }_2{\psi }_1|{\mu }^{(N+2)}⟩,{\beta }_{12}^{\varkappa }=⟨0^{(N)}|{\psi }_2^{†}{\psi }_1^{†}|{\mu }^{(N-2)}⟩$$

$${\partial }_t{G}_{11^{\prime }}(t-t^{\prime })=-i\delta (t-t^{\prime })⟨[{\psi }_1(t),{{\psi }^{†}}_{1^{\prime }}(t^{\prime }){]}_{+}⟩+$$

$$+⟨T[H,{\psi }_1](t){{\psi }^{†}}_{1^{\prime }}(t^{\prime })⟩,$$

$$[H,{\psi }_1](t)=e^{iHt}[H,{\psi }_1]e^{-iHt}.$$

$$[H,{\psi }_1]=-{\epsilon }_1{\psi }_1+[V,{\psi }_1],$$

$$(i{\partial }_t-{\epsilon }_1)G_{11^{\prime }}(t-t^{\prime })={\delta }_{11^{\prime }}\delta (t-t^{\prime })+i⟨T[V,{\psi }_1(t)]{{\psi }^{†}}_{1^{\prime }}(t^{\prime })⟩.$$

$$R_{11^{\prime }}(t-t^{\prime })=i⟨T[V,{\psi }_1](t){{\psi }^{†}}_{1^{\prime }}(t^{\prime })⟩.$$

$$R_{11^{\prime }}(t-t^{\prime })\overleftarrow{{\partial }_{t^{\prime }}}$$

$$[H,{{\psi }^{†}}_{1^{\prime }}]={\epsilon }_{1^{\prime }}{{\psi }^{†}}_{1^{\prime }}+[V,{{\psi }^{†}}_{1^{\prime }}],$$

$$R_{11^{\prime }}(t-t^{\prime })(-i\overleftarrow{{\partial }_{t^{\prime }}}-{\epsilon }_{1^{\prime }})$$

$$G_{11^{\prime }}(\omega )=G_{11^{\prime }}^{(0)}(\omega )+\sum _{22^{\prime }}{G}_{12}^{(0)}(\omega )T_{22^{\prime }}(\omega )G_{2^{\prime }{1}^{\prime }}^{(0)}(\omega ),$$

$$T_{11^{\prime }}(t-t^{\prime })$$

$$T_{11^{\prime }}^{(0)}(t-t^{\prime })$$

$$T_{11^{\prime }}^{(r)}(t-t^{\prime })$$

$$G(\omega )=G^{(0)}(\omega )+G^{(0)}(\omega )T(\omega )G^{(0)}(\omega ),$$

$$G(\omega )=G^{(0)}(\omega )+G^{(0)}(\omega )\Sigma (\omega )G(\omega )$$

Full text

I Introduction

Response of many-body quantum systems to external perturbations encompasses a large class of problems actively studied in various areas of quantum physics. The observed characteristics of response are associated with quantum correlation functions, which encode the complete information about strongly coupled complex media. The general notion of correlation functions forms the common underlying background connected across the fields of quantum chromodynamics, dynamics of hadrons, astrophysics, condensed matter, solid state, atomic, molecular, nuclear structure physics, and formal aspects of quantum field theory (QFT). Thus, advancements toward an exact theory of correlations in strongly coupled quantum systems have potentially a broad impact on many fields of research.

Response of atomic nuclei to various external perturbations represents a very rich playground to study correlation functions due to the availability of numerous experimental probes, which was further extended with the advent of rare beam facilities Tanihata (1998); Glasmacher et al. (2017). It has a long history of theoretical studies based on QFT. It was recognized rather early that the random phase approximation (RPA) Bohm and Pines (1951) appears as a good approach to the gross features of nuclear spectra, such as the positions and total strengths of collective excitations, however, the response theory which, in principle, should provide the full spectral composition has to be extended by configurations beyond one-particle one-hole ($1p-1h$) ones considered in RPA. The idea of coupling between single-particle and emergent collective degrees of freedom in nuclei Bohr and Mottelson (1969, 1975); Broglia and Bortignon (1976); Bortignon et al. (1977); Bertsch et al. (1983); Soloviev (1992), which explained successfully many of the observed phenomena, was, in fact, linked to the non-perturbative versions of QFT-based and, in principle, exact, equations of motion for correlation functions in nuclear medium Ring and Schuck (1980). Indeed, the equation of motion method developed in Ref. Rowe (1968) and further elaborated, e.g., in Refs. Schuck (1976); Adachi and Schuck (1989); Danielewicz and Schuck (1994); Dukelsky et al. (1998); Schuck and Tohyama (2016a), was shown to produce a hierarchy of approximations to the dynamical kernels of the equations for one-fermion and two-time two-fermion propagators. In particular, the non-perturbative versions of those kernels, which include full resummations in the particle-hole and particle-particle channels, can be mapped to the kernels of the phenomenological nuclear field theories (NFT), where such kernels are commonly referred to as particle-vibration coupling (PVC) or quasiparticle-phonon models (QPM) Bohr and Mottelson (1969, 1975); Broglia and Bortignon (1976); Bortignon et al. (1977); Bertsch et al. (1983); Soloviev (1992). This mapping, although it has to be corrected for the accurate lowest-order limit, provides an understanding of the emergent collective phenomena, explaining the mechanism of their formation from the underlying strongly interacting degrees of freedom. Moreover, the EOM method gives clear insights into the relationship between the bare nucleon-nucleon interaction and its modification in the strongly coupled medium. It reveals, in particular, that the latter is not reducible to static ’potentials’, but splits into a static part calculable from the bare interactions beyond the Hartree-Fock approximation Schuck and Tohyama (2016a) and a well-defined dynamical component Adachi and Schuck (1989); Dukelsky et al. (1998).

Although these aspects were widely ignored in the majority of semi-phenomenological PVC models based on effective in-medium interactions over the years Bortignon et al. (1978); Bortignon and Broglia (1981); Bertsch et al. (1983); Mahaux et al. (1985); Bortignon et al. (1986); Bortignon and Dasso (1997); Colò and Bortignon (2001); Tselyaev (1989); Kamerdzhiev et al. (1997); Ponomarev et al. (2001, 1999); Lo Iudice et al. (2012), these and other approaches provided invaluable knowledge about the importance of coupling between single-nucleon and collective degrees of freedom in nuclear structure. Lately, this type of approaches was linked to the contemporary density functional theories Litvinova and Ring (2006); Litvinova et al. (2008, 2010); Tselyaev (2013); Afanasjev and Litvinova (2015); Tselyaev et al. (2016); Niu et al. (2015, 2018), advancing the PVC models to self-consistent frameworks, and applied to experimental data analyses Endres et al. (2010); Poltoratska et al. (2014); Massarczyk et al. (2012); Lanza et al. (2014); Özel-Tashenov et al. (2014).

However, very little progress has been made on the conceptual advancements of the many-body aspects of the non-perturbative NFT’s. Some specific rare topics, such as PVC with charge-exchange phonons Litvinova (2016); Robin and Litvinova (2018) or PVC-induced ground state correlations Kamerdzhiev et al. (1991); Tselyaev (1989); Kamerdzhiev et al. (1997); Robin and Litvinova , have been addressed in the time blocking approximation (TBA), however, relatively little effort has been made on developments and numerical implementations of NFT’s beyond the two-particle two-hole $2p-2h$ level Tselyaev et al. (2018); Litvinova (2015), although the phenomenological multiphonon approach Soloviev (1992); Ponomarev et al. (2001); Savran et al. (2006); Andreozzi et al. (2008) indicates the possibility to meet the shell-model standards in large model spaces. Also, the problem of consistent linking those approaches to the underlying bare interactions remains unsolved.

The goal of the present work is to bring together the past and recent developments of the EOM method and to compare them to the PVC approach in the time blocking approximation. As the EOM method starts from the bare interaction between fermions, while the existing versions of PVC models for atomic nuclei are based on the effective static interactions, it only makes sense to compare the dynamical kernels of the equations for the two-fermion propagators. We will show that such a comparison allows one to recognize and justify the use of cluster expansions truncated on the two-body level in the PVC models and to substantiate their extensions to higher-order configurations. We will also discuss another important feature of the EOM method, namely its capability of deriving the emergence of the long-range correlations and the corresponding collective degrees of freedom from the bare short-range correlations in a parameter-free way. Eventually, by confronting the EOM method with the class of approaches of the PVC type, we develop a non-perturbative, consistent and systematically improvable theory for one- and two-fermionic correlation functions, which is directly based on the bare two-fermion interactions.

We want to keep the presentation sufficiently general to be useful also in other than just the nuclear physics field. The dynamical aspects of the in-medium fermionic correlations, in particular, the PVC are surely also of prime importance in electronic condensed matter systems where, e.g., we know that particle-phonon coupling can even reverse the sign of the repulsive Coulomb interaction to give rise to superconductivity Cooper (2010); Broglia and Zelevinsky (2013). Also, even though in atoms and molecules collective features are not well born out, particle-hole correlations are important in screening the Coulomb force Tiago et al. (2008); Martinez et al. (2010); Sangalli et al. (2011). Similar techniques as the one we outline here have indeed been presented recently in the field of chemical physics Olevano et al. (2018).

In Sections II - IV the EOM’s are reviewed following the formalism of Refs. Schuck (1976); Adachi and Schuck (1989); Danielewicz and Schuck (1994); Dukelsky et al. (1998); Schuck and Tohyama (2016a), Section V discusses the PVC model in the time blocking approximation (PVC-TBA) in terms of Refs. Kamerdzhiev et al. (1997); Litvinova et al. (2007a) and its mapping to the EOM, Section VI is devoted to a non-perturbative advancement of the theory for the particle-hole response beyond $2p-2h$ configurations, and Section VII presents the first numerical implementations of the extended theory. Finally, Section VIII provides conclusions and outlook. In Table 1 we collect the most important symbols used in the article with the references to the corresponding equations and figures, where they are introduced and employed, to help navigate throughout the text.

II One-fermion and two-fermion propagators in a correlated medium

A convenient way of describing a strongly correlated many-body system is calculating various correlation functions, or propagators, which are directly linked to observables. For instance, the one-fermion in-medium propagator, or Green function, is defined as follows:

[TABLE]

where $T$ is the chronological ordering operator, $\psi(1),{\psi^{\dagger}}(1)$ are one-fermion (for instance, one-nucleon) fields:

[TABLE]

and the subscript ’1’ denotes the full set of the one-fermion quantum numbers in an arbitrary representation. The fermionic fields obey the usual anticommutation relations:

[TABLE]

The averaging in Eq. (1) $\langle...\rangle$ is performed over the formally exact ground state of the many-body system of $N$ particles described by the Hamiltonian $H$:

[TABLE]

Here, the operator $H^{(1)}$ is the one-body contribution to the Hamiltonian:

[TABLE]

whose matrix elements $h_{12}$, in general, combine the kinetic energy $t$ and the mean-field $v^{(MF)}$ part of the interaction. The operator $V^{(2)}$ describes the two-body sector associated with the two-fermion interaction

[TABLE]

the operator $V^{(3)}$ generates the three-body forces

[TABLE]

with the antisymmetrized matrix elements ${\bar{v}}_{1234}$ and ${\bar{w}}_{123456}$, respectively, and so on. We will make an explicit derivation of the equations of motion assuming that the Hamiltonian is confined by the two-body interaction, however, the theory can be naturally extended to multiparticle forces. The inclusion of a three body force, as may be necessary in nuclear physics for an ab initio approach, can formally be handled straightforwardly, but at the cost of inflating the amount of formulas. Therefore, we will refrain from doing this here, but may consider an extension to the case of the presence of the three-body forces in a future work.

It is often convenient to work in the basis, which diagonalizes the one-body (also named single-particle) part of the Hamiltonian (5), so that $h_{12}=\delta_{12}\varepsilon_{1}$. We will use this basis from the beginning, while on the way to the final equations of motion it will be redefined as soon as the mean-field part of the Hamiltonian will absorb additional contributions from the two-body sector. Furthermore, for the case of self-bound finite Fermi systems the common practice is to work with intrinsic Hamiltonians which are obtained by isolating the kinematics of the center of mass from the relative motion. In this case the two-body part of the Hamiltonian absorbs the kinetic energy of the relative motion, see for instance, Ref. Hergert et al. (2016) and references therein for details.

The single-particle propagator (1) depends explicitly on a single time difference $\tau=t-t^{\prime}$, and the Fourier transform with respect to $\tau$ to the energy domain leads to the spectral (Lehmann) expansion:

[TABLE]

This expansion is a sum of simple poles with the residues composed of matrix elements of the field operators between the ground state $|0^{(N)}\rangle$ of the $N$-particle system and states $|n^{(N+1)}\rangle$ and $|m^{(N-1)}\rangle$ of the neighboring systems with $N+1$ and $N-1$ particles, respectively:

[TABLE]

By definition, these matrix elements give the weights of the given single-particle (single-hole) configuration on top of the ground state $|0^{(N)}\rangle$ in the $n$-th ($m$-th) state of the $(N+1)$-particle ($(N-1)$-particle) system, respectively. The residues are associated with the observable occupancies of the corresponding states, often called spectroscopic factors. The poles are located at the energies $E^{(N+1)}_{n}-E^{(N)}_{0}$ and $-(E^{(N-1)}_{m}-E^{(N)}_{0})$ of those systems, respectively, i.e., related to the ground state of the $N$-particle system.

The two-fermion, three-fermion and, in general, $n$-fermion propagators, or Green functions, are defined in analogy to Eq. (1):

[TABLE]

The response of a many-body system to external perturbations, which can be associated with one-body operators, is expressed via the two-time two-fermion particle-hole propagator (response function):

[TABLE]

where we imply that $t_{1}=t_{2}=t,t_{1^{\prime}}=t_{2^{\prime}}=t^{\prime}$. The Fourier transformation of Eq. (12) to the energy (frequency) domain leads to the spectral expansion

[TABLE]

which, similarly to the one for the one-fermion propagator (8), satisfies the general quantum field theory requirements of locality and unitarity. The residues of this expansion are products of the matrix elements

[TABLE]

called transition densities, which are properly normalized and represent the weights of the pure particle-hole configurations on top of the ground state $|0\rangle$ in the model (ideally, exact) excited states $|\nu\rangle$ of the same N-particle system. The poles are the excitation energies $\omega_{\nu}=E_{\nu}-E_{0}$ relative to the ground state.

The response to the probes with pair transfer is associated with the two-time two-fermion Green function (II) whose spectral expansion reads

[TABLE]

where the residues are composed of the matrix elements

[TABLE]

and the poles $\omega_{\mu}^{(++)}$ and $\omega_{\varkappa}^{(--)}$ are formally exact states of the $(N+2)$- and $(N-2)$-particle systems, respectively. The spectral expansions of Eqs. (8,13,15) are model independent and valid for any physical approximations to the many-body states $|n\rangle,|m\rangle$, $|\nu\rangle$, $|\mu\rangle$, and $|\varkappa\rangle$. The sums in Eqs. (8,13,15) are formally complete, i.e., run over the discrete spectra and imply continuous integrals over the continuum states.

Thus, one can see that, indeed, due to their direct links to the observables, both one-fermion and two-fermion Green functions as well as the particle-hole response function are of high importance for the characterization of strongly correlated quantum many-body systems, in particular, of atomic nuclei. In the next two sections we consider equations of motion for the one-fermion Green function $G_{11^{\prime}}(t-t^{\prime})$ and for the particle-hole response $R_{12,1^{\prime}2^{\prime}}(t-t^{\prime})$ and outline possible strategies on the way to their accurate descriptions. As the correlation functions provide, in principle, complete characteristics of quantum many-body systems, essentially the same theory is applicable to atomic physics, quantum chemistry, condensed matter and other areas of physics dealing with fermionic systems.

III Equation of motion for one-fermion propagator

III.1 General formalism

Let us consider the time evolution of the one-fermion propagator (1). Taking the derivative with respect to $t$, one obtains

[TABLE]

where we defined

[TABLE]

With the help of the commutator

[TABLE]

assuming that $V\equiv V^{(2)}$ and discarding the three-body interaction, Eq. (17) leads to

[TABLE]

Let us introduce a function $R_{11^{\prime}}(t-t^{\prime})$ corresponding to the last term on the right hand side of Eq. (20):

[TABLE]

As we will see in the following, it is useful to determine the equation of motion for this function with respect to $t^{\prime}$:

[TABLE]

Using the commutator

[TABLE]

one arrives at the EOM for $R_{11^{\prime}}(t-t^{\prime})$:

[TABLE]

Combining it with the first EOM (20) and performing the Fourier transformation to the energy (frequency) domain with respect to the time difference $t-t^{\prime}$, we obtain

[TABLE]

where we introduced the free (uncorrelated) one-fermion propagator $G^{(0)}_{11^{\prime}}(\omega)=\delta_{11^{\prime}}/(\omega-\varepsilon_{1})$ and the interaction kernel (one-body T-matrix, not to be confused with the time ordering operator):

[TABLE]

Here and in the following we use the superscript ”(0)” to denote the static parts of the interaction kernels and ”(r)” for their dynamical time-dependent parts, which are associated with retardation effects in our approach. The EOM (25) which, in the operator form, is

[TABLE]

can be transformed to the Dyson equation:

[TABLE]

with the interaction kernel $\Sigma(\omega)$, such as

[TABLE]

from which it follows that the operator $\Sigma$ represents the one-fermion self-energy (also called mass operator) as the irreducible (with respect to one-fermion line) part of the kernel $T$: $\Sigma=T^{irr}$. Analogously to Eq. (26), the self-energy is decomposed into the instantaneous mean-field part $\Sigma^{(0)}$ and the energy-dependent dynamical part $\Sigma^{(r)}(\omega)$:

[TABLE]

Notice here that the decomposition of the kernels (26,30) into the static and time- (energy-) dependent, or dynamical, parts is a generic feature and the direct consequence of the time-independence of the bare interaction $V$ of Eq. (6).

The first static (instantaneous) terms of both kernels coincide and read:

[TABLE]

Here we need first to evaluate the commutator $[V,\psi_{1}]$ which, with help of the anticommutation relations (3), can be obtained as

[TABLE]

where the Latin indices have the same meaning as the number indices and the definition of the antisymmetrized interaction matrix elements ${\bar{v}}_{1234}=v_{1234}-v_{1243}$ was taken into account. Evaluating the anticommutator

[TABLE]

one gets:

[TABLE]

Thus, the first (instantaneous) part $\Sigma^{(0)}$ of the mass operator (30) is associated with the mean field contribution:

[TABLE]

where $\rho_{li}=\langle{\psi^{\dagger}}_{i}\psi_{l}\rangle$ is the ground-state one-body density and we have applied the (anti)symmetry properties of the antisymmetrized interaction matrix elements: ${\bar{v}}_{1234}=-{\bar{v}}_{1243}=-{\bar{v}}_{2134}={\bar{v}}_{2143}$. The second (dynamical) part $\Sigma^{(r)}(\omega)$ of the mass operator comprises all retardation effects induced by the nuclear medium.

In order to understand the dynamical part $\Sigma^{(r)}(\omega)$ of the self-energy $\Sigma(\omega)$, let us first evaluate its reducible counterpart $T^{(r)}_{11^{\prime}}(t-t^{\prime})$. Here we can use the result of Eq. (32) for the commutator $[V,\psi_{1}]$, and the following result for the second commutator:

[TABLE]

so that

[TABLE]

or, returning to the number indices,

[TABLE]

where we assume that $t_{2}=t_{3}=t_{4}=t$ and $t_{2^{\prime}}=t_{3^{\prime}}=t_{4^{\prime}}=t^{\prime}$, as dictated by the instantaneous interaction. With the help of of the three-fermion Green function Eq. (LABEL:Tr) can be rewritten as

[TABLE]

Here we realize that, although the EOM for one-fermion propagator $G(\omega)$ (28) is formally a closed equation, however, its interaction kernel $\Sigma(\omega)=T^{irr}(\omega)$ is defined by Eqs. (26), (35) and (39), where its time-dependent part contains the three-fermion Green function $G(432^{\prime},23^{\prime}4^{\prime})$. Up to here the theory is exact, but, in order to calculate the three-body propagator, one would need to generate a series of equations of motion for higher-rank propagators. However, with a very good accuracy the problem can be truncated at the two-body level. The three-fermion Green function, according to Refs. Martin and Schwinger (1959); Vinh Mau ; Mau and Bouyssy (1976); Ring and Schuck (1980), can be (approximately) decomposed as follows 111Accuracy of this approximation will be examined elsewhere.:

[TABLE]

where

[TABLE]

contains all uncorrelated three-body contributions. This kind of truncation of the many-body EOM hierarchies is referred to as cluster-expansion approach in condensed matter physics Kira and Koch (2006, 2011). Both Eqs. (LABEL:G3) and Eq. (41) include terms with one-fermion propagators connecting the same time moments, either $t$ or $t^{\prime}$. As these propagators form closed loops, obviously such contributions are connected by one-fermion lines, i.e., are reducible. Thus, the irreducible part of the three-fermion propagator of Eq. (LABEL:G3) reads:

[TABLE]

The two-body propagators in Eq. (42), such as $G(32^{\prime},23^{\prime})$, contain both uncorrelated and correlated terms. It is convenient to introduce its connection to the response functions $R(12^{\prime},21^{\prime})$ and to the fully correlated parts of the two-body propagators ${\bar{R}}(12^{\prime},21^{\prime})$:

[TABLE]

thus related as follows:

[TABLE]

Collecting the terms associated with $\tilde{R}$ components, one can separate fully correlated and fully uncorrelated $G^{(0)}$ parts of $G^{irr}(432^{\prime},23^{\prime}4^{\prime})$ according to:

[TABLE]

The corresponding dynamical part of the self-energy shown diagrammatically in Fig. 1 is

[TABLE]

Using the symmetry properties of the interaction $\bar{v}$ and renaming the indices $3\leftrightarrow 4,3^{\prime}\leftrightarrow 4^{\prime}$ in the third term, $4^{\prime}\leftrightarrow 3^{\prime}$ in the fourth term, $4\leftrightarrow 3$ in the fifth term and $4\leftrightarrow 3$ in the last term, one obtains

[TABLE]

This form of the self-energy allows for a direct reduction to the second-order (with respect to the bare interaction $\bar{v}$) approach for the dynamical self-energy. Indeed, by dropping the fully correlated terms with $\tilde{R}$ the lowest-order approach takes the form

[TABLE]

which corresponds to the last two terms in Fig. 1. The response theory built on this self-energy leads to an approach called second random phase approximation (SRPA), which has been implemented within the frameworks employing both microscopic Papakonstantinou and Roth (2009, 2010) and effective Gambacurta et al. (2010, 2011); Gambacurta and Grasso (2016) interactions and demonstrates a good ability of generating fragmentation effects on nuclear excitation spectra. Closely related are the methods based on a stochastic one-body transport theory, that distinguish and include in-medium nucleon-nucleon collisions of both coherent and incoherent nature Lacroix et al. (2001).

III.2 Emergent phonons and particle-vibration coupling

The prerequisites for calculating the self-energy beyond the perturbation theory are the correlation functions ${\bar{R}}$ in the particle-hole and particle-particle channels. In practice, however, it is more convenient to deal with the particle-particle Green function and particle-hole response function which contain both correlated and uncorrelated contributions (44). By isolating the particle-hole response function, the self-energy can be reorganized as follows:

[TABLE]

where we separate the terms with the particle-particle (pp) and particle-hole (ph) Green functions coupled to one-hole and one-particle ones, respectively, from the term with the uncorrelated two-particle-one-hole propagator. This equation can serve as foundation for microscopic approaches to the single-particle self-energy, which refer to the phenomenon of particle-vibration coupling, or PVC. The correlation functions $G^{pp}$ and $R^{ph}$ represent the emergent degrees of freedom, phonons, which are the quasibound states of two fermions embedded in the strongly-correlated medium. In nuclear physics the phonons are associated with nuclear vibrations and, in relatively large systems, such as medium-mass and heavy nuclei, can acquire a collective character. They are often called ”vibrations” and attributed to the vibrational motion of the nuclear surface, although, as we will see in this subsection, in the formally exact theory they are more general objects.

The connection to the PVC can be seen more explicitly if one identifies the correlation functions contracted with the interaction matrix elements with the phonon propagators and coupling vertices, as displayed diagrammatically in Fig. 2. For this purpose it is convenient to work with the Fourier image of $\Sigma^{(r)}_{11^{\prime}}(t-t^{\prime})$ in the energy domain

[TABLE]

The first term of Eq. (49) transforms as follows:

[TABLE]

where we denoted the single-particle energies as $\varepsilon_{n}^{(+)}=E^{(N+1)}_{n}-E^{(N)}_{0}$ and $\varepsilon_{m}^{(-)}=E^{(N-1)}_{m}-E^{(N)}_{0}$, and defined the pairing phonon vertex functions according to:

[TABLE]

Then, introducing the amplitude $\Gamma^{pp}_{12,1^{\prime}2^{\prime}}(\omega)$,

[TABLE]

with the pairing phonon propagator

[TABLE]

Eq. (51) can be alternatively obtained by the convolution

[TABLE]

Similarly, the Fourier image of the second term of Eq. (49)

[TABLE]

can be obtained either directly or by introducing the following mapping:

[TABLE]

with the phonon vertices $g^{\nu}$ and propagators $D_{\nu}(\omega)$:

[TABLE]

Then it can be shown that

[TABLE]

In Eqs. (51-59) the index $\sigma=\pm 1$ stands for the forward (”particle”) and backward (”hole”) components of the phonon propagators and vertices, and the spectral representations (8,15) along with the definitions (1,II,12) were applied. Finally, the last term of the self-energy (49) with only uncorrelated single-particle Green functions transforms as follows:

[TABLE]

The complete dynamical part of the one-fermion self-energy (49) is shown in Fig. 3 in the diagrammatic form in terms of the particle-vibration coupling. The first two terms on the right hand side are formed by the topologically similar one-loop diagrams which are analogous to the electron self-energy corrections in quantum electrodynamics, where electron emits and reabsorbs a photon, or to the nucleonic self-energy of quantum hadrodynamics where a single nucleon emits and reabsorbs a meson. Here, they represent the effects of a strongly correlated medium, where a single fermion emits and reabsorbs a phonon of the particle-particle (first term) and the particle-hole (second term) nature, in addition to the emission and reabsorption of an uncorrelated two-particle-one-hole configuration (third term). The mappings established by Eqs. (53,57) and illustrated diagrammatically in Fig. 2 explain the underlying mechanism of the induced in-medium interaction, where the emergent composite bosons are formed by quasibound fermionic pairs.

One may notice that the dynamical self-energy $\Sigma^{(r)}$ recast in the form of Eq. (47) with the separation of fully correlated and non-correlated parts helps to relate the approach to the lowest-order perturbation theory and to assess the role of correlations. It is clear that in the case of weak coupling the uncorrelated term(s) play the leading role and the phonon-exchange interaction can be neglected, however, in the strong-coupling regime the phonon coupling dominates and the lowest-order uncorrelated term does not give the leading contribution. Indeed, in the major applications to nuclear systems only the first two terms are taken into account and, furthermore, coupling to the pairing phonons was found much less important than coupling to the normal particle-hole phonons. These approximations, however, were shown to be justified only within the methods based on the effective nucleon-nucleon interactions. Such interactions are typically obtained by fitting the bulk nuclear properties, such as their masses and radii, on the mean-field level, i.e., assuming that one-fermion self energy contains solely the static part (35) with only the one-body density, which is implicitly coupled to correlations in the dynamical part of the self-energy (49). This coupling is essentially important as

[TABLE]

which means that the one-fermion density matrix entering Eq. (35) is the equal-times limit of the full solution of Eq. (28). This fact is often expressed in terms of the density dependencies of the effective interactions while these dependencies are typically disconnected from a detailed analysis of Eq. (28) with the complete kernel $\Sigma(\omega)$. Inevitably, existing versions of the PVC model which add the dynamical part on top of the effective interactions imply an additional subtraction procedure to remove the double counting of PVC which is contained in the static approximation in the parameters of the phenomenological mean field Tselyaev (2013). Such a subtraction turned out to be a very elegant way of avoiding the double counting, instead of the complicated refitting of the parameters of the mean field, and it is widely applied in calculations of two-body Green functions, in particular, the particle-hole response discussed below. However, an analogous procedure has not been formulated for the case of the one-body propagator.

Regardless what kind of the two-body interaction is used for calculations of the one-fermion Green function from Eq. (28), for an accurate solution beyond the static approximation to the interaction kernel the knowledge about the particle-hole response and particle-particle Green function is needed, as follows from Eqs. (28,30,35,49). As we discuss in the next section following Refs. Dukelsky et al. (1998); Olevano et al. (2018), it is possible to formulate the equation of motion for these two-fermion Green functions in a similar manner as for the one-fermion propagator and to obtain non-perturbative approximations to its closed form. Below we discuss the EOM approach to the particle-hole response function, while the EOM for the particle-particle response can be generated in complete analogy.

IV Equation of motion for the particle-hole response

The equation of motion for the particle-hole response of Eq. (12) can be generated by taking its time derivative with respect to $t$:

[TABLE]

After the evaluation of the first commutator and the commutator with the one-body part of the Hamiltonian, the EOM (63) takes the form

[TABLE]

with the norm kernel ${\cal N}_{121^{\prime}2^{\prime}}$,

[TABLE]

In the case when the one-body density matrix is diagonal, it takes the form ${\cal N}_{121^{\prime}2^{\prime}}=\delta_{11^{\prime}}\delta_{22^{\prime}}(n_{1}-n_{2})\equiv\delta_{11^{\prime}}\delta_{22^{\prime}}{\cal N}_{12}$, where $n_{1}=\langle{\psi^{\dagger}}_{1}\psi_{1}\rangle$ is the occupation number of the fermionic state $1$. In Eq. (64) and in the following $\varepsilon_{12}=\varepsilon_{1}-\varepsilon_{2}$ is the difference between the matrix elements of the one-body Hamiltonian, or the single-particle energies, in the same representation.

The second EOM is generated by differentiating the last term on the right hand side of Eq. (64) with respect to $t^{\prime}$:

[TABLE]

Combining it with the first EOM (64) and performing the Fourier transformation to the energy domain, one obtains

[TABLE]

where the uncorrelated particle-hole response $R^{(0)}(\omega)$ is defined as

[TABLE]

The integral part is determined by the Fourier transform $T_{12,1^{\prime}2^{\prime}}(\omega)$ of the interaction kernel $T_{12,1^{\prime}2^{\prime}}(t-t^{\prime})$ (not to be confused with the time ordering operator)

[TABLE]

which splits naturally into the instantaneous $T^{(0)}$ and the time-dependent $T^{(r)}$ parts

[TABLE]

with the shorthand notation for the product of the diagonal norm kernels ${\tilde{\cal N}}_{121^{\prime}2^{\prime}}={\cal N}_{12}{\cal N}_{1^{\prime}2^{\prime}}$. In the operator form Eq. (LABEL:2bgfb) reads

[TABLE]

The latter equation can be further transformed to a formally closed equation for $R(\omega)$, similar to the Dyson equation for one-fermion Green function. The kernel of this new equation should be irreducible with respect to uncorrelated particle-hole response $R^{(0)}$, which means that

[TABLE]

where

[TABLE]

or $K(\omega)=T^{irr}(\omega)$, i.e. that $K(\omega)$ absorbs the irreducible contributions of $T(\omega)$.

Obviously, the kernel $K(\omega)$ can be also decomposed into the instantaneous (static) and time-dependent (frequency-dependent) terms:

[TABLE]

In a complete analogy to the case of one-fermion EOM, the decomposition of the interaction kernel (LABEL:Ft,74) into the static and time(energy)-dependent, or dynamical, parts is a generic feature of the in-medium interaction in the particle-hole channel and the direct consequence of the time-independence of the bare interaction $V$ of Eq. (6).

Evaluation of the static part of Eq. (LABEL:F1) with the help of the following commutator

[TABLE]

gives for the internal single commutator of $T^{(0)}$:

[TABLE]

and for the double commutator, using again the result of Eq. (75):

[TABLE]

Thus, the static part of Eq. (LABEL:Ft) reads:

[TABLE]

where we introduced the two-fermion density $\rho_{ij,kl}$,

[TABLE]

so that $\sigma^{(2)}_{ij,kl}$ represents its fully correlated part. After that, the static kernel takes the form:

[TABLE]

where the first term isolates the contribution from the bare interaction and the norm factor will be compensated by its inverse, according to Eq. (LABEL:Ft). In transforming Eq. (78) to Eq. (80), the remaining terms with the single-particle mean field can be absorbed in the single-particle energies by replacing $\varepsilon_{1}\to{\tilde{\varepsilon}}_{1}=\varepsilon_{1}+\Sigma^{(0)}_{11}$ in the uncorrelated response function of Eq. (68) and redefining the working basis. From Eq. (80) it becomes clear that, in the absence of correlations contained in the quantities $\sigma^{(2)}$ and $T^{(r)}$, the EOM (72) takes the form of the well-known random phase approximation. The complete static part of the effective two-fermion interaction is shown diagrammatically in Fig. 4 in terms of the full two-body densities $\rho$ represented by the rectangular blocks.

The time-dependent part $T^{(r)}$ of the kernel $T$ can be evaluated by making use of the commutator (76):

[TABLE]

where we have decomposed the kernel $T^{(r)}_{12,1^{\prime}2^{\prime}}(t-t^{\prime})$ into the four terms with different general structure:

[TABLE]

or, returning to the number indices,

[TABLE]

where the two-time two-particle-two-hole (four-fermion) Green function $G(543^{\prime}1^{\prime},5^{\prime}4^{\prime}31)$ was introduced according to:

[TABLE]

The components of the irreducible dynamical kernel $K^{(r)}_{12,1^{\prime}2^{\prime}}(t-t^{\prime})=T^{(r)irr}_{12,1^{\prime}2^{\prime}}(t-t^{\prime})$ are shown diagrammatically in Fig. 5, where we have omitted the factors $\pm i/4$ in front of each diagram.

The general result of Eqs. (LABEL:Fr-83) for the dynamical kernel is known in the literature, see, for instance, Refs. Schuck and Ethofer (1973); Dukelsky et al. (1998). One can see that the dynamical kernel of Eq. (72) requires the knowledge about the two-time two-particle-two-hole (four-body) propagator $G$ of Eq. (LABEL:G4). An exact treatment of the time-dependent kernel would require generating EOM’s for the three-fermion or four-fermion propagators, thus building a hierarchy of EOM’s which is equivalent to that known as Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy. In this work we consider a truncation of this hierarchy on the level of two-fermion correlation functions, as proposed, in particular, in Refs. Schuck and Tohyama (2016a); Olevano et al. (2018). This approach avoids generating EOM’s for the three-fermion and four-fermion propagators and, as we will show below, leads to a closed system of equations for the two-fermion correlation functions.

As all four components of $T^{(r)}(t-t^{\prime})$ (83) contain the same two-particle-two-hole propagator, it is sufficient to analyze in detail one of them, then the other three can be reconstructed straightforwardly. Let us consider $K^{(r;11)}(t-t^{\prime})$ as the reference component keeping only the irreducible (with respect to the particle-hole uncorrelated propagator) contributions to $T^{(r)}(t-t^{\prime})$. The lowest-order approximation is determined by the uncorrelated irreducible part $G^{(0)irr}$ of the four-fermion propagator:

[TABLE]

where

[TABLE]

Thus, the lowest-order approximation requires only the one-fermion (mean-field) Green functions:

[TABLE]

and their backward-going counterparts. The uncorrelated part $K^{(r;11)0}(t-t^{\prime})$ of the irreducible kernel $K^{(r;11)}(t-t^{\prime})$ is given diagrammatically in Fig. 6. Using the uncorrelated one-fermion Green function:

[TABLE]

with $t_{11^{\prime}}=t_{1}-t_{1^{\prime}}$ and $\sigma_{1}=\pm 1$ for particle/hole states, the Fourier transform of Eq. (LABEL:G0irr) can be calculated, so that

[TABLE]

where $\delta\to+0$, the explicit term corresponds to the sum of a) and c) contributions, and the antisymmetrized term ’$\cal AS$’ to the sum of b) and d) contributions shown in Fig. 6.

The next step in building the dynamical part of the kernel is to include nonperturbative resummations between each pair of fermionic lines. In this approximation, the correlated contributions to Eq. (LABEL:G4) are represented by products of the two-body correlated and uncorrelated propagators

[TABLE]

where the term $\cal AS$ absorbs all possible antisymmetrizations and the index ’0’ marks an uncorrelated ground state. In turn, resummations within the remaining uncorrelated fermionic pairs lead to the terms with two two-body correlators:

[TABLE]

The latter approach to the dynamical kernel of Eq. (72) is, thus, the most complete one within the concept of the cluster-expansion truncation on the two-body level. The approximation of Eq. (91) only neglects the additional uncorrelated contribution, which should be not very important in the strong-coupling regimes, as in the case of the one-fermionic kernel (49), and the contributions with explicit correlation functions of three and more fermions. One can see that the nonperturbative approximations (90,91) to the four-fermion propagator (LABEL:G4) contain the particle-hole response functions as well as the two-particle and two-hole Green functions. The singly correlated contributions from $G^{(c)irr}$ to the kernel $K^{(r;11)}$ are shown in Fig. 7, and the doubly correlated ones in Fig. 8. At this stage it becomes clear that for the determination of the dynamical part of the interaction kernel one needs the knowledge about the two-body correlation functions of the particle-hole and particle-particle (hole-hole) types, which enter Eqs. (90,91). To formalize this, we can complement the particle-hole propagator considered above, by the particle-particle and hole-hole ones. Thus, we can define

[TABLE]

Then, the EOM for ${\hat{R}}(\omega)$ can be written in a compact form:

[TABLE]

Remarkably enough, all frequency (time) dependence of the kernel originates from the internal two-fermion propagators which are themselves the main variables. Such a closed approach has been discussed, in particular, in the EOM language for the two-fermion Green functions Adachi and Schuck (1989); Danielewicz and Schuck (1994); Dukelsky et al. (1998) and time-dependent density matrices in Refs. Schuck and Tohyama (2016a, b). A similar approach to the nuclear response has been developed over the years as the method of chronological decoupling of diagrams or the time blocking approximation Tselyaev (1989); Kamerdzhiev et al. (1997); Litvinova and Tselyaev (2007); Tselyaev (2007); Litvinova et al. (2007a, 2008) which has become self-consistent in its later implementations based on meson-nucleon Lagrangians Litvinova et al. (2007a, 2008, 2013); Litvinova (2015, 2016); Robin and Litvinova (2016, 2018). This approach starts from the general Bethe-Salpeter equation for a four-time two-fermion Green function, but after applying a certain time projection technique reduces to the two-time or single-frequency equation of motion of the Dyson type with the kernel, which is topologically equivalent to selected components of the EOM singly-correlated kernel. This model is discussed in detail in the next Section. Analogous approaches of the nuclear field theory, although derived differently, also lead to single-frequency EOM’s with frequency-dependent kernels containing couplings between single-fermion and phonon propagators Broglia and Bortignon (1976); Bortignon et al. (1977); Bertsch et al. (1983); Colò and Bortignon (2001); Niu et al. (2012). All these methods are not related explicitly to the bare nucleon-nucleon interactions and based on phenomenological descriptions of the mean field, instead of the static part of the kernel (78,80), and effective in-medium nucleon-nucleon interactions. The EOM method for two-fermion response functions, taking into account both static and dynamical parts of the interaction kernel, is now being more applied in quantum chemistry Olevano et al. (2018).

It is not likely that Eq. (93) with the dynamical interaction kernels of Eqs. (90,91) can be solved analytically except for some toy models Schuck and Tohyama (2016a, b). In practice, its self-consistent solution can be found iteratively. To initialize an iterative algorithm, one would need a starting approximation to the response functions contained in the kernel $K[{\hat{R}}(\omega)]$. In the cases of phenomenological models discussed above and in the next Section, the starting-point response functions are approximated by the solution of Eq. (93) keeping only the static term in the kernel, that corresponds to the random phase approximation with effective interactions.

V Particle-vibration coupling (PVC) model in the time blocking approximation (TBA)

In this Section we revisit and investigate the physical content of the approach to the particle-hole response function of Eq. (72), known as particle-vibration coupling model, in the context of the EOM method discussed above. In the Green function language, the PVC model was formulated in Refs. Tselyaev (1989); Kamerdzhiev et al. (1997) within a method of chronological decoupling of diagrams. It emerged as an extension of the Migdal’s theory of finite Fermi systems Migdal (1967) beyond the quasiparticle random phase approximation (QRPA). The extended theory was generalized for the case of superfluid pairing in Refs. Tselyaev (2007); Litvinova and Tselyaev (2007) and received a relativistic formulation and a fully self-consistent (parameter-free) implementation in Refs. Litvinova et al. (2007a, 2008); Marketin et al. (2012); Robin and Litvinova (2016). In the latter versions and numerous later developments and applications the method included the same idea of chronological decoupling of diagrams and was renamed to time blocking approximation (TBA).

The final expressions of the PVC-TBA approach for the dynamical kernel of Eq. (72) turned out to be equivalent to those of the NFT Bertsch et al. (1983); Colò and Bortignon (2001) obtained by making use of the perturbation theory for the coupling between the quasiparticles and collective doorway states. Both approaches were based on the assumption about the dynamical part of the nucleonic self-energy in the EOM form of Eq. (49), where the bare interaction is replaced by the effective interaction, and typically kept only the second term of it, which occurred to be dominant. The intermediate fermionic line in this term was approximated by the mean-field one-fermion Green function. The conventional PVC-TBA and NFT response theories were confined by the $ph\otimes phonon$ configurations following the idea of a small parameter hidden in the effective PVC vertices. Although PVC-TBA showed a considerable improvement of the description of nuclear excited states compared to QRPA, its unclear foundation and uncontrollable approximations remained the drawbacks preventing this approach from further development. Later, based on formal similarities with the EOM kernels Schuck (1976), the PVC-TBA method was generalized to the case of $2phonon$ configurations Tselyaev (2007) implemented in the relativistic framework in Refs. Litvinova et al. (2010, 2013). Further extensions were discussed in Refs. Litvinova (2015); Tselyaev et al. (2018), however, without a comparison to the EOM method and without a detailed analysis of the assumptions intuitively based on the EOM derivation. Thus, the goals of this and subsequent Sections are (i) to explicitly compare the EOM with the dynamical kernel discussed in Section IV with the PVC-TBA and, with the insights of this comparison, (ii) develop a systematic approach to the response functions with a non-perturbative treatment of higher configurations in a strongly correlated medium.

Below we will follow the formalism and conventions of Refs. Kamerdzhiev et al. (1997); Litvinova et al. (2007a), so that the index mapping will require the replacement $R(12,34)\to R(21,43)$ in order to compare the equations for the particle-hole response function with those of the previous section. Note that in the PVC-TBA model, one starts from the Bethe-Salpeter equation (BSE) for a more general four-time response Tselyaev (1989); Kamerdzhiev et al. (1997):

[TABLE]

where the summation over the number indices $1$, $2,\dots$ implies integration over the respective time variables, and $V$ is the interaction amplitude irreducible in the $ph$-channel. This amplitude is, in general, a variational derivative of the one-fermion self-energy $\Sigma$ with respect to the exact single-particle Green function:

[TABLE]

and a four-time analog of the kernel of Eq. (74). Similarly, it is a sum of the static effective interaction $\tilde{V}$ and the time-dependent (energy-dependent) one $V^{(e)}$:

[TABLE]

where we implied that $t_{12}=t_{1}-t_{2}$. Despite the fact that in this approach one starts from the formally exact BSE (LABEL:bse) and employs the decomposition (96) of the interaction kernel derived by the EOM method, it does not imply a connection with the underlying bare interaction. Instead, the static part of the self-energy is adjusted to experimental data, for instance, the data on nuclear binding energies and radii assuming that they depend only on one-body density or fitted to reproduce the lowest single-particle excitations obtained in knock-out or transfer reactions. Thus, in this theory,

[TABLE]

where the dynamical part of the self energy is borrowed from that of the EOM (49), neglecting the uncorrelated part, performing the mapping (57,58) with the replacement $\bar{v}\to{\tilde{V}}$ and Fourier transformation to the energy domain. In Eq. (97) ${\tilde{\varepsilon}}_{{1}}$ are the eigenvalues of the phenomenological one-body part of the Hamiltonian. Usually the PVC models do not include coupling to pairing phonons, but in principle it can be added as a term similar to the last line of Eq. (97). Notice that, although both the self-energy (97) and the interaction kernel (96) are again decomposed into the static and dynamical parts, in the effective theory these parts are not equivalent to those in the EOM starting from the bare interaction. However, we can call them topologically equivalent as they have the same internal propagator structure.

In terms of the free response, which is, in the time domain, a product of two one-fermion propagators $R^{0}(12,34)=G(1,3)G(4,2)$, the BSE (LABEL:bse) can be written in the operator form as

[TABLE]

Since the static mean-field part of the interaction kernel is fixed by fitting the global characteristics of the many-body system to data, it is convenient for further analysis to eliminate the exact Green function $G$ and rewrite it in terms of the mean field Green function $\tilde{G}$, such as

[TABLE]

with

[TABLE]

where $\delta\to+0$. Using the connection between $\tilde{G}$ and $G$ in the Nambu form

[TABLE]

one can rewrite Eq. (98) as follows:

[TABLE]

with the uncorrelated mean-field particle-hole response function ${\tilde{R}}^{0}(12,34)={\tilde{G}}(1,3){\tilde{G}}(4,2)$ and $W$ as a new interaction kernel of the form

[TABLE]

where

[TABLE]

The BSE in this form is more convenient to consider since the mean-field Green function ${\tilde{G}}$ is well defined. However, the interaction kernel $W$ in Eq. (104) becomes more complicated. The graphical representation of the Eq. (103) with the interaction kernel defined by Eqs. (104,105) is shown in Fig. 9.

Thus, one can see that the dynamical part of the kernel takes the form, which is similar to the one obtained in the EOM, but still depends on four times. The last term $i\Sigma^{e}(1,3)\Sigma^{e}(4,2)$ of Eq. (105) can be still related to one of the doubly correlated contributions of the EOM, but in the present context it looks very different from the others. As it is discussed in Ref. Kamerdzhiev et al. (1997), this term is of a higher order and compensates multiple counting of the particle-phonon coupling arising from the two previous terms, if the entire dynamical kernel is still unrestricted by only particle-hole pairs of indices. Indeed, terms with backward-going correlation functions in the dynamical kernel are also possible, but they require a special consideration. They have been analyzed, in particular, within a phenomenological approach and included in calculations of nuclear neutral excitations in Refs. Tselyaev (1989); Kamerdzhiev et al. (1991, 1997) as well as more recently within a self-consistent relativistic framework applied to charge-exchange excitations in Ref. Robin and Litvinova .

In the leading resonant time blocking approximation, which is discussed below, the possibility of having particle-particle and hole-hole states as well as the connection between particle-hole and hole-particle states in the dynamical kernel is neglected. This corresponds to the absence of ground state correlations more complex than the particle-hole ones. In this case multiple counting does not take place and the term $i\Sigma^{e}(1,3)\Sigma^{e}(4,2)$ can be omitted.

Let us consider the Fourier transformation of the Eq. (103) to the energy domain. The response function formally depends on four time variables, but, in fact, on three time differences, because of the time translational invariance. Thus, a triple Fourier transform is needed to translate the BSE (103) into the equation with respect to energy variables. In order to obtain the spectral representation of the response, two of them have to be integrated out. These operations lead to the following equation:

[TABLE]

with the subsequent integration over $\varepsilon$:

[TABLE]

One can notice, however, that both the solution of Eq. (LABEL:bsee) $R$ and its kernel $W$ are singular with respect to the energy variables. This is related to the fact that Eq. (103) contains integrations over all time points of the intermediate states. This means that many configurations which are actually more complex than $1p1h\otimes$phonon are contained in the exact response function. In Ref. Tselyaev (1989) a special time-projection technique was introduced to block the $ph$-propagation through these complex intermediate states. It has been shown that for this type of response it is possible to reduce the integral equation (LABEL:bsee) to a relatively simple algebraic equation. Below we will see that this approximation corresponds to retaining certain part of the terms specified in Eq. (90) in the internal Green functions propagators of the dynamical two-time kernel in the EOM method.

The Bethe-Salpeter equation (103) can be solved in two steps. First one can calculate the correlated particle-hole response $R^{(e)}$, which describes the particle-hole propagation under the influence of the interaction $W^{(e)}$

[TABLE]

and contains all the effects of particle-phonon coupling and all the singularities of the integral part of the main BSE (103). Second, the remaining equation for the full response function $R$

[TABLE]

which contains only the static effective interaction $\tilde{V}$, can be easily solved when $R^{(e)}$ is known. Thus, the main problem to address singularities is to calculate the correlated particle-hole response $R^{(e)}$. The latter can be represented as an infinite series of graphs with uncorrelated $ph$-propagators alternated with single interaction events:

[TABLE]

where ${\Gamma}^{(e)}$ is, thus, a reducible analog of $W^{(e)}$ containing correlated two-particle-two-hole blocks connected by the uncorrelated $ph$-propagators. Then, to avoid higher-complexity configurations, a time-projection operator

[TABLE]

with the Heaviside type functions $\theta(t)$ is introduced into the integral part of Eq. (111) Tselyaev (1989) according to:

[TABLE]

so that

[TABLE]

After the Fourier transformation in time we restrict ourselves to the two-time response function $R_{12,34}(\omega)$, because it has to be subsequently contracted with equal-times external field operators:

[TABLE]

which depends only on one energy variable $\omega$. As a result, we obtain an algebraic equation for the spectral representation of the particle-hole response:

[TABLE]

where

[TABLE]

${\tilde{\varepsilon}}_{12}={\tilde{\varepsilon}}_{1}-{\tilde{\varepsilon}}_{2}$ and $\Phi(\omega)$ is the particle-phonon coupling amplitude with the following components: $ph-ph$ matrix element has the form

[TABLE]

where $\{1,1^{\prime},1^{\prime\prime}\}$ are the particle states with ${\tilde{\varepsilon}}_{1},{\tilde{\varepsilon}}_{1^{\prime}},{\tilde{\varepsilon}}_{1^{\prime\prime}}>\varepsilon_{F}$, $\{2,2^{\prime},2^{\prime\prime}\}$ are the hole states with ${\tilde{\varepsilon}}_{2},{\tilde{\varepsilon}}_{2^{\prime}},{\tilde{\varepsilon}}_{2^{\prime\prime}}\leq\varepsilon_{F}$ and $\varepsilon_{F}$ is the Fermi energy. The hp-hp matrix elements $\Phi_{12,1^{\prime}2^{\prime}}^{(hp,hp)}(\omega)$ are obtained by Hermitian conjugation and time-reversal transformation $\omega\to-\omega$. The diagrammatic representation of $\Phi(\omega)$ is shown in Fig. 10.

Thus, we have obtained the expression for the dynamical interaction kernel which can be compared to the kernel of the EOM method. Indeed, a complete matching can be revealed by looking at the parts of the latter kernel associated with the singly-correlated terms and performing the mapping defined by Eq. (58). The matching is illustrated in Figs. 11,12 and can be additionally verified by taking the Fourier transform of the dynamical kernel with the internal Green functions of the type (90) Olevano et al. (2018). One should be, however, aware of the differences between ${\bar{v}}$ of the ab-initio theory and ${\tilde{V}}$ used in effective theories as well as of the double-counting removal correction needed in the latter case.

In the framework of an effective theory, indeed, an important correction has to be done to the Eq. (116). Being adjusted to experimental data, the static interaction $\tilde{V}$ contains, in principle, all beyond-mean-field correlations, in particular, those which are explicitly included into the dynamical kernel ${\tilde{V}}^{(e)}$, in the static approximation. Therefore, a direct addition of the dynamical interaction leads to a double counting of the static correlations contained in $\Phi(\omega)$. In order to avoid this, the dynamical kernel has to be corrected as follows Tselyaev (2013):

[TABLE]

The physical meaning of this subtraction is clear: the average value of the particle-vibration coupling amplitude $\Phi$ in the ground state is supposed to be contained already in the residual effective interaction $\tilde{V}$, therefore, we should take into account only the additional energy dependence, i.e., $\delta\Phi(\omega)=\Phi(\omega)-\Phi(0)$, on top of this effective interaction. Instead of Eq. (116), we finally solve the following response equation:

[TABLE]

In early applications of PVC-TBA Tselyaev (1989); Kamerdzhiev et al. (1997) the phonon coupling vertices $g^{\nu}$ were calculated based on experimental information about deformation parameters for the lowest collective excitations. Their experimental energies were taken as the $\Omega_{\nu}$ values. The common practice was to include only very few collective phonons, to justify neglecting the uncorrelated term in the dynamical kernel. In later applications, such as in Refs. Litvinova and Tselyaev (2007); Litvinova et al. (2007a, 2008); Lyutorovich et al. (2008); Litvinova et al. (2010); Niu et al. (2012); Litvinova et al. (2013); Niu et al. (2014); Lyutorovich et al. (2015); Tselyaev et al. (2016); Robin and Litvinova (2016), the phonon spectra were extracted from the RPA solutions, i.e., the solutions of Eq. (116) without the dynamical kernel (more specifically, keeping only the first term $\tilde{V}$ in the interaction kernel of Eq. (116)). These works adopted larger phonon spaces and the subtraction of Eq. (120) elaborated in detail in Ref. Tselyaev (2013), thus, representing an important step toward a closed calculation scheme for the particle-hole response function containing dynamical medium effects in nuclei.

The first fully self-consistent version of the PVC-TBA was implemented in the relativistic framework based on the effective meson-nucleon Lagrangian of quantum hadrodynamics Litvinova et al. (2007a). In this approach, the most important particle-hole phonons again correspond to the solutions of the same equation (116) without the dynamical kernel. In a more complete version also particle-particle and hole-hole (pairing) phonons should be included as well as the phonons appeared in the form of proton-neutron (charge-exchange) correlated pairs. Contribution of the charge-exchange phonons was investigated within the relativistic formalism in Refs. Litvinova (2016); Robin and Litvinova (2018), where their role was found somewhat smaller than that of the regular like-particle ones, but sizable enough to modify the single-particle structure Litvinova (2016) and spin-isospin strength distributions Robin and Litvinova (2018). The contribution of pairing vibrations to the dynamical kernel of Eq. (72) is commonly neglected, however, a systematic study of those effects is highly desirable to properly assess their role.

But even without the inclusion of the pairing channels, accurate solutions of Eqs. (72,116) require an iterative scheme, as discussed in Section IV. This was not realized in the conventional PVC-TBA, however, the idea of iterative solutions is straightforward and can be, in principle, implemented numerically. In the calculations based on the effective interaction derived from a reasonably good density functional theory the most important phonons are described relatively well already in RPA. Their characteristics do not change significantly (except for acquiring fragmentation and larger widths by the high-frequency and soft modes, which is expected to give second-order effects to the dynamical kernel) in the calculations beyond RPA. This has been verified in Ref. Tselyaev et al. (2018) by direct calculations, where the phonons obtained as full solutions of Eq. (116) were recycled and reused in the dynamical PVC-TBA kernel (116) to compute the dipole response of some medium-mass and heavy nuclei, which remained almost unaffected by these non-linear effects. At the same time, it follows from numerous studies performed with this type of the dynamical kernel that there is a clear need of extensions beyond its configuration complexity. Indeed, as we have shown in Section II, the response function should describe, in principle, the complete excitation spectrum of a quantum many-body system, while in numerical applications of the conventional PVC-TBA method with $1p1h\otimes phonon$ configurations the deficiency of the level density is quite obvious after the comparison to data Endres et al. (2010); Özel-Tashenov et al. (2014).

In this context it becomes interesting to look at the terms which are missing in the conventional PVC-TBA, but present in the EOM of Section IV. One can see after performing the matching shown in Fig. 11 that the singly-correlated terms of $K^{(r;11)}$, such as b), e) and f) of Fig. 7 and their counterparts from $K^{(r;12)}$, $K^{(r;21)}$, and $K^{(r;22)}$ as well as all doubly correlated terms of Fig. 8 and their respective counterparts are not included in the original PVC-TBA. A way to include the terms of such structure to some extent was formulated in Tselyaev (2007) as a two-phonon version of the approach based on the EOM Schuck (1976), where both phonon correlation functions in the dynamical kernel of the type (91) were taken in the random phase approximation. Such a possibility has been also discussed in the context of a formal comparison between the models where the internal Green functions in the dynamical kernel are confined by only the uncorrelated terms of the type (LABEL:G0irr), the singly-correlated terms (90) and those of Eq. (91) Olevano et al. (2018). In the PVC-TBA framework, the two-phonon approach was implemented numerically for the dipole nuclear response and analyzed in Refs. Litvinova et al. (2010, 2013). In the latter work, a quantitative comparison between the results obtained within the approximations of Eq. (119) (the kernel of Fig. 11) and Eq. (91) with two coupled RPA phonons has been made, which showed some improvements of the description of nuclear strength functions, when the doubly-correlated dynamical kernel is used. However, the major deficiencies, which occur because of the too small overall number of excited states in the resulting response function (13) remained unsolved, because this approach does not go beyond the correlated $2p-2h$ configurations.

VI Doubly correlated dynamical kernel beyond $2p-2h$ configurations

The direct way to overcome the latter problem is to calculate the dynamical kernel beyond $2p-2h$ configurations. Indeed, the Green function of Eq. (91) entering the dynamical kernel consists of coupled two-fermion correlation functions, however, it is not limited by any approximation to these correlation functions. Ideally, they should provide a convergent solution of Eq. (93) and contain, in principle, the entire excitation spectrum. We have already mentioned two attempts to go beyond the conventional PVC-TBA kernel of Eq. (119) and Fig. 10, which were investigated in Refs. Litvinova et al. (2010) and Tselyaev et al. (2018). The latter uses the phonons computed beyond RPA as solutions of the full Eq. (116) and the former employs the kernel of Eq. (91) with two RPA response functions. As mentioned above, both models still left some room for improvement.

A different strategy has been outlined in Ref. Litvinova (2015) as a generalized PVC-TBA, which proposes an iterative algorithm for the solution of Eq. (116). After determining the RPA phonons, like in the conventional PVC-TBA method, we can solve the equation of motion (116) for the particle-hole response function $R(\omega)$, which contains configurations of the $ph\otimes phonon$ or two-quasiparticles coupled to phonon ($2q\otimes phonon$) type, for various multipolarities, and after that reiterate the dynamical kernel as follows:

[TABLE]

where the quantities $\zeta$ are the phonon vertex matrices:

[TABLE]

so that the resulting four terms

[TABLE]

correspond to the four diagrams in Fig. 13 with $n=2$. The index $\eta=\pm 1$ denotes upper and lower components in the quasiparticle space, see Refs. Litvinova (2015); Tselyaev (2007) for more details.

Thus, the amplitude $\Phi^{(3)}$ contains the contributions of the graphs shown in Fig. 14 in all orders with respect to the internal propagators. However, the proposed procedure allows calculations of their contribution without explicit evaluation of the diagrams of Fig. 14. It is straightforward to see that these terms contain $2q\otimes 2phonon$ configurations and thereby represent the next, three-particle-three-hole ($3p3h$), level of the configuration complexity, as compared to the original PVC-TBA, which in this implementation for superfluid systems in the relativistic framework was named relativistic quasiparticle time blocking approximation (RQTBA). Thus, we adopt EOM/R(Q)TBA3 as a working name for the approach of Eq. (124).

The amplitude $\Phi^{(3)}$ forms the dynamical interaction kernel for the correlated particle-hole propagator $R^{(e)(3)}$ taking into account 3p3h correlations (to be compared to $R^{(e)(2)}\equiv R^{(e)}$ which includes 2p2h ones):

[TABLE]

Analogously to the conventional PVC-TBA and RQTBA, the remaining equation for the full response function is formulated in terms of the correlated particle-hole propagator $R^{(e)(3)}(\omega)$ as a free term and the static effective interaction as a kernel:

[TABLE]

where the superscript ’$(3)$’ indicates that this response function takes into account $3p-3h$ configurations. Analogously to the $2q\otimes phonon$ RQTBA, the subtraction of the amplitude $\Phi^{(3)}$ at zero frequency from the effective interaction in Eq. (126) eliminates the double counting of the static contribution of the phonon coupling effects. It also improves convergence of the sums in Eq. (124) with respect to the phonon and quasiparticle energies of the intermediate complex configurations Tselyaev (2013). Finally, the observed strength function is a (model-independent) double convolution of the response function with the external field operator:

[TABLE]

The imaginary part $\Delta$ of the energy variable in Eq. (127) is often used in the numerical calculations as a smearing parameter. It provides a smooth envelope of the strength distribution and averages over complex configurations which are not taken into account explicitly. This parameter can also mimic the experimental resolution of the data invoked for the comparison.

Calculations presented below were performed within the approach of Eqs. (116,122-127), however, in principle, the iteration procedure can be continued until convergence is achieved. The initial steps are characterized as follows:

[TABLE]

Then, the chain of operator equations for the correlated particle-hole propagator $R^{(e)(n)}$, the phonon coupling amplitude $\Phi^{(n)}$ and the response function $R^{(n)}$ reads

[TABLE]

where the index ’$(n)$’ with $n\geq 2$ indicates the iteration step. The phonon vertex matrices $\zeta^{\nu}$ can be, in principle, recalculated on each step using the spectral representation of the response function (13) while, as discussed above, the study of Ref. Tselyaev et al. (2018) showed that the effect of such corrections may be small. In the last equation (129) $\Phi^{(n)}(0)$ is subtracted from the static interaction kernel to avoid double counting effects for the case of calculations based on effective interactions. In ab-initio calculations this subtraction is absent. The proposed iterative method allows one to obtain the contributions of the three, four , and higher PVC-loop diagrams without calculating them explicitly. The generalized dynamical kernel is illustrated diagrammatically in Fig. 13.

There are various factors which favor the convergence of the iterative procedure with respect to the internal response functions (129). In an analogy to Eliashberg theory for electron-phonon interaction in solids Eliashberg (1960), in finite nuclei a certain smallness is contained in the phonon vertices which enter the dynamical kernel, that serves as a foundation for the phenomenological PVC approaches Bohr and Mottelson (1975). Indeed, in spherical nuclei the PVC vertices enter the dynamical kernel in combinations ${\bar{g}}^{\nu}_{12}\sim\langle 1||g^{\nu}||2\rangle/(\Omega_{\nu}\sqrt{2j_{1}+1})$ Bohr and Mottelson (1975); Tselyaev (1990); Kamerdzhiev et al. (1997). Thus, on each iteration, a factor of $\sim|{\bar{g}}^{\nu}|^{2}\ll 1$ suppresses higher-order contributions. The equation for the response function is commonly solved for certain angular momentum and parity, which correspond to those quantum numbers of the external field. For this purpose, Eqs. (122-129) have been formulated in the coupled form in Ref. Litvinova (2015), where it is shown that the kernel $\Phi^{(n)}(\omega)$ contains a geometrical factor which enters $\Phi^{(n)}(\omega)$ as a multiplier on each iteration. This factor shows up as a product of two 6j-symbols that contains a smallness originated from the coupling of the angular momentum of the external channel with the angular momenta of the two internal correlation functions. As it is shown below, in the numerical calculations the fragmentation effects generated by the $2q\otimes phonon$ configurations on the two-quasiparticle R(Q)RPA states are considerably stronger than those caused by the subsequent addition of the $2q\otimes 2phonon$ configurations, that can serve as a clear indication of convergence.

VII Numerical details and results

The EOM/RQTBA3 model, which was originally proposed in Ref. Litvinova (2015) and briefly revisited in the previous section, has been implemented numerically and tested in calculations of nuclear dipole response. This type of response is known to dominate nuclear spectra and associated with the largest corpus of available experimental data. The original as well as the evaluated data on the giant dipole resonance (GDR), the high-energy part of the dipole response, can be found in Ref. nnd . Typically the GDR above the particle emission threshold and its low-energy counterpart are measured with different methods, although the newer techniques, such as inelastic proton scattering Tamii et al. (2011); Poltoratska et al. (2014), allow for unified high-quality measurements. The knowledge about dipole strength distributions in nuclei is crucial for many applications in nuclear sciences and astrophysics, see more details in Refs. Paar et al. (2007); Savran et al. (2013). Overall, testing both the high- and low-energy dipole strength distributions is the best benchmark for newly developed many-body models.

As in the base RQTBA model Litvinova et al. (2008), we implement a multi-step parameter-free calculation scheme, but with a few more steps now: (i) the closed set of the relativistic mean field (RMF) Hartree-Bogoliubov equations Serot and Walecka (1986); Ring (1996); Vretenar et al. (2005) are solved using the NL3 parametrization of Refs. Boguta and Bodmer (1977); Lalazissis et al. (1997) for the non-linear sigma-meson model and monopole pairing forces adjusted to reproduce empirical pairing gaps. The obtained single-particle Dirac spinors and the corresponding single-nucleon energies, being the eigenstates and eigenvalues of the relativistic mean-field Hamiltonian, formed the working basis for subsequent calculations of the response, where the same effective NL3 meson-exchange interaction $\tilde{V}$ is also adopted; (ii) the relativistic quasiparticle random phase approximation (RQRPA) equation Paar et al. (2003), which is equivalent to Eq. (116) without the dynamical kernel, is solved to obtain the phonon vertices $g^{\nu}$ and their frequencies $\Omega_{\nu}$. The set of phonons with the $J_{\nu}^{\pi_{\nu}}$ = 2+, 3-, 4+, 5-, 6+ and frequencies $\Omega_{\nu}\leq$ 15 MeV, together with the RMF single-particle basis, forms the $2q\otimes phonon$ configurations for the particle-phonon coupling amplitude ${\Phi}(\omega)$, while the phonon space was additionally truncated according to the values of the reduced transition probabilities of the corresponding electromagnetic transitions: the modes with the values of the reduced transition probabilities $B(EL)$ less than 5% of the maximal one (for each $J_{\nu}^{\pi_{\nu}}$) were neglected. These are the common truncation criteria for the PVC models based on the RMF, see, for instance Litvinova (2016), where a convergence study was presented; (iii) Eq. (116) for the response function is solved in the truncated configuration space, which includes excitations within the energy window of interest 0-25 MeV, for spins and parities $J^{\pi}$ = 0*±* - 6*±* (in this first application we, however, neglected the static part of the kernel of Eq. (116), just to study the power of the correlated internal particle-hole propagator to induce additional fragmentation - in further applications the full kernel will be taken into account); (iv) the obtained response functions, together with the previously obtained phonon characteristics, are used to compute the dynamical kernel of Eq. (124); (v) the correlated particle-hole propagator with $J^{\pi}=1^{-}$ is obtained according to Eq. (125) in the same truncated configuration space; (vi) the full dipole response function is computed by solving Eq. (126) in the momentum space, as described in Litvinova et al. (2007a, 2008); (vii) finally, the strength function is found according to Eq. (127) with the external field operator of the electromagnetic dipole (EME1) character:

[TABLE]

where $Z$ and $N$ are the numbers of protons and neutrons, respectively, $A=N+Z$, and $e$ is the proton charge. The sums in Eq. (LABEL:opE1) are performed over the corresponding nucleonic degrees of freedom. While the response to this operator is often attributed to electromagnetic probes, hadronic probes are associated predominantly Lanza et al. (2014) with the response to the isoscalar dipole (ISE1) operator

[TABLE]

where $\eta=5\langle r^{2}\rangle/3$ and the second term in the brackets eliminates the spurious translational mode Garg and Colò (2018).

A sufficiently large quasiparticle basis in both Fermi (particle) and Dirac (antiparticle) sectors should be used in solving Eqs. (116, 126). Although the dynamical kernels $\Phi(\omega)$, which induce fragmentation of two-quasiparticle configurations, may be cut off outside a window confined by the energy of interest, the static kernel $\tilde{V}$ has to be included in the complete or nearly complete two-quasiparticle space Litvinova et al. (2008). The latter kernel is responsible for the correct location of the simple (R(Q)RPA) modes of the strength distribution and associated mainly with the medium-range correlations, while the former kernels introduce the long-range effects causing the redistribution of the strength. Here the completeness means that the two-quasiparticle basis, in which Eqs. (116, 126) are solved, should include all the single-quasiparticle states which participate in the RMF self-consistent procedure. In our case the basis spans the single-quasiparticle states with the angular momenta up to $41/2$ in both Fermi and Dirac sectors, the same range where the parameters of the Lagrangian have been fitted. The respective energy range of the two-quasiparticle excitations is confined by $\sim$250 MeV in the Fermi sector and by $\sim$-1950 MeV in its Dirac counterpart. Keeping the complete two-quasiparticle basis in Eqs. (116, 126) guarantees full self-consistency, in particular, the proper decoupling of the dipole translational mode from the intrinsic dipole excitations in R(Q)RPA Ring and Schuck (1980); Paar et al. (2003). This fact can be verified numerically, for instance, by calculating the isoscalar dipole strength distribution produced by the response to the operator of Eq. (131), where the radial form factor is corrected for the center of mass motion. Without this correction one typically sees a dominant peak located at zero energy, as we show in the left panel of Fig. 15 in comparison to the right panel, where the response to the corrected isoscalar operator is displayed, for 48Ca. One can see that in the latter case the zero-energy translational mode is suppressed. Moreover, this property is kept in the extended R(Q)TBA and EOM/R(Q)TBA3 models - indeed, the subtraction procedure of Eqs. (120, 126) leads to the purely R(Q)RPA kernel in the $\omega\to 0$ limit. This feature is known since early implementations of the time blocking method with the subtraction Litvinova and Tselyaev (2007). However, in R(Q)TBA and EOM/R(Q)TBA3 the translational mode may be fragmented because, like the physical states, it can be coupled to the phonons. Although, due to the subtraction procedure, the main peak of the translational mode remains at zero energy, its fragments may spread around it. As the excited states calculated with the uncorrected and corrected isoscalar dipole operators look different in both the R(Q)TBA and EOM/R(Q)TBA3, in the present implementations these models do not guarantee complete decoupling of the spurious mode from the physical states. A solution to this problem was proposed in Ref. Tselyaev (2014) in the form of a projection operator applied to the dynamical kernel, that prevents the coupling of complex configurations to the translational mode and, thus, removes its admixture to the physical states. Performing this transformation is beyond the scope of the present article, but will be considered in future work. The sensitivity of our present implementation to the two-quasiparticle basis incompleteness was inspected and revealed that an energy cut-off of this basis by $\sim$100 MeV in the Fermi sector and $\sim$-1800 MeV in its Dirac counterpart does not introduce noticeable changes in the excitation spectra, that can be used for more economical calculations.

The results of calculations for the electromagnetic dipole response in 42,48Ca are displayed in Fig. 16. The strength distribution obtained within EOM/RQTBA3 (red solid curves) is plotted against the results of RQRPA (black dot-dashed curves) and RQTBA (blue dashed curves) and compared to experimental data (green curves and circles) of Ref. nnd in terms of the dipole photoabsorption cross section

[TABLE]

i.e. with an additional energy factor in front of the strength distribution. These cross sections were investigated in the RQTBA framework in Ref. Egorova and Litvinova (2016) in order to establish the role of the $2q\otimes phonon$ configurations in the formation of the spreading width of the GDR in the chain of calcium isotopes. We found that these configurations indeed result in the significant spreading width improving considerably the agreement to data as compared to RQRPA, while the latter approach provides a very good description of the GDR’s position and the total strength. However, although we used a fairly large model space of the $2q\otimes phonon$ configurations, the total width of the RQTBA dipole strength distribution still underestimates its experimental value. Another major shortcoming was found on the high-energy shoulder of the GDR above its centroid, where the cross sections were systematically underestimated. A similar situation was reported in Ref. Tselyaev et al. (2016), where a systematic downshift of the non-relativistic QTBA strength distributions, as compared to the RPA ones, was revealed in calculations with various Skyrme forces for monopole, dipole and quadrupole resonances in both light and heavy nuclei. Now, when EOM/RQTBA3 with more complex $2q\otimes 2phonon$ configurations becomes available, we can see that it shows the potential of resolving those problems. Indeed, Fig. 16 shows that the new configurations present in EOM/RQTBA3 induce a stronger fragmentation of the GDR and its additional spreading toward both higher and lower energies. Another new effect is a relatively small, but a visible shift of the main peak toward higher energies. These features of changing the high-energy behaviour of the strength distributions can be directly related to the appearance of the new higher-energy complex configurations and, consequently, the higher-energy poles in the resulting response functions, that has the power to rearrange the energy balance of the overall strength distribution. This change is, however, modest enough to conserve the energy-weighted sum rule (EWSR) which varies quantitatively only within a few percent as compared to the RQRPA and RQTBA ones in the investigated energy region.

As it has been shown in Ref. Tselyaev (2007), in the resonant time blocking approximation, where the high-energy behaviour of the dynamical kernel is $\Phi(\omega)\sim 1/\omega$ at $\omega\to\infty$, the EWSR is exactly the same as in the (Q)RPA. In other words, the presence of the dynamical kernel of Eq. (119) with such a behaviour does not violate the EWSR of RQRPA. It can be seen from Eqs. (12,124) that the iterations of the dynamical kernel do not change its high-energy asymptotics as it is determined by the internal particle-hole propagators satisfying the model-independent Eq. (12). The presence of subtraction in Eq. (125), however, changes the static part of the kernel and, thus, introduces a small violation of the EWSR as in the conventional time blocking approximation of Eqs. (120,121) Litvinova and Tselyaev (2007); Litvinova et al. (2007a, 2008). As it is known in the literature, in extended models with the ground state correlations caused by PVC another sum rule violation occurs, as well as its restoration Adachi and Lipparini (1988); Tselyaev (2007). The formalism of Ref. Delion et al. (2016) proves that the self-consistent RPA (SCRPA) and possibly also the second SCRPA keeps all desirable qualities of the standard RPA intact, that includes also the sum rules. The complex ground state correlations are not included in the present work, but were recently addressed in other developments Robin and Litvinova .

The obtained change of the high-energy behaviour of the strength distributions in EOM/RQTBA3 may also provide some new arguments to the discussion of fitting the nuclear energy density functionals. In particular, using the position of the GDR as one of the reference observables may become difficult because of its model dependence.

The low-energy behaviour of the dipole strength distribution in nuclei has been a topic of an intense research during the last couple of decades. A strong astrophysics connection of this type of strength, in particular, to the r-process nucleosynthesis, attracted much of interest from both experimental and theoretical sides. In this context, both RQRPA and RQTBA models were examined for their performance and for their potential of describing the low-energy dipole strength associated with the pygmy dipole resonance (PDR), or the soft dipole mode Litvinova et al. (2007b, 2009a, 2009b); Endres et al. (2010, 2012); Massarczyk et al. (2012); Litvinova et al. (2013); Özel-Tashenov et al. (2014). Similarly to the GDR case, RQTBA provided an improved description of the PDR showing a considerably richer spectral structure because of the fragmentation effects induced by the $2q\otimes phonon$ and $2phonon$ configurations. However, in many cases it became clear that even calculations with quite large model spaces reveal a deficiency of configuration complexity that lead to too low level density in the discrete and quasicontinuum energy sectors below the particle emission threshold. Therefore, it is interesting to examine the newly-developed EOM/RQTBA3 for its performance in the low-energy regime.

Fig. 17 illustrates the performance of $R(Q)TBA$ and EOM/R(Q)TBA3 in the description of the isoscalar (top) and isovector (bottom) low-energy dipole response of 48Ca. The experimental particle emission threshold of 48Ca is $\sim$9.95 MeV nnd , so that the observed spectrum below this energy has a discrete character. Two kinds of experiments, namely the inelastic photon scattering ($\gamma,\gamma^{\prime}$) and ($\alpha,\alpha^{\prime}\gamma$) coincidence dominated by the strong interaction have been performed and discussed in Refs. Hartmann et al. (2004); Derya et al. (2014). It was found, in particular, that the sub-threshold dipole response of 48Ca does not show an isospin splitting as observed in heavier neutron-rich nuclei Endres et al. (2010), where the lowest-lying states are mostly isoscalar in nature and the states at higher energy are of the mixed isoscalar-isovector character. Here we do not intend to go into the details of the underlying structure of the states, which can be essentially model dependent, but rather look at the potential of our many-body models to reproduce the experimental strength distribution. For this purpose, we have extracted the reduced transition probabilities from the computed strength distributions by making use the simple relationship: $B_{\nu}=\pi\Delta S(\omega_{\nu})$, which follows from Eq. (127) and which is valid for sufficiently small values of $\Delta$ (here we used $\Delta$ = 20 keV) in the vicinity of the peak $\omega_{\nu}$. The first observation here is that R(Q)RPA does not produce any excited states in this energy region and gives its first solution just above 10 MeV. The inclusion of the 2q$\otimes$phonon configurations in R(Q)TBA improves the picture considerably: due to the fragmentation effects caused by the coupling of the R(Q)RPA modes to these configurations some strength comes down below the neutron threshold. The cumulative EME1 strength below 10 MeV agrees with the experimental data quite well as it was discussed in detail in Ref. Egorova and Litvinova (2016), however, the major fraction of the R(Q)TBA strength remains in the threshold area. On the other hand, the ISE1 strength distribution in R(Q)TBA shows the pattern which is rather close to the experimental one (here we plot the energy-weighted values in order to compare with Ref. Derya et al. (2014)). It can be also seen that adding more complex configurations within the EOM/R(Q)TBA3 further improves the description: it (i) slightly rearranges the ISE1 strength bringing it even closer to the observed one and (ii) gives two distinct EME1 states below 9 MeV, which can be approximately matched the experimental ones, and brings more strength from the higher-energy region to the threshold area, where also several experimental levels are found at about the same energies. The experimentally extracted BE1 values of the states in the threshold area may strongly depend on the methods employed in the data analysis Tonchev et al. (2010), so that those BE1 values may not be a good benchmark for the theory.

Another interesting case is the low-energy dipole response of the unstable neutron-rich 68Ni, which has been investigated both experimentally Wieland et al. (2009); Rossi et al. (2013); Wieland et al. (2018) and theoretically Litvinova et al. (2010, 2013); Papakonstantinou et al. (2015). The results of our calculations for the low-energy dipole response of 68Ni are displayed in Figs. 18, 19. As in the case of calcium isotopes, we compare the results of the three models, RQRPA, RQTBA and EOM/RQTBA3, with the same curve- and color-coding as in Fig. 16. Fig. 18 shows the strength functions calculated with a small value of the smearing parameter $\Delta=$20 keV, thus allowing for a clear illustration of the fragmentation mechanism in a parameter-free many-body theory with an effective static interaction. In particular, one can see how the spectrum of RQTBA emerges from a relatively poor one of RQRPA, which is essentially the single strong and relatively collective state at $~{}9.5$ MeV. The addition of $2q\otimes phonon$ configurations of RQTBA results in the blue curve which is obviously the fragmented major RQRPA state spread over a broader energy region. Remarkably, the RQTBA strength no longer shows any clear dominance of a single state, but is rather spread uniformly over the $7-15$ MeV energy interval. Finally, when we add the EOM/RQTBA3 strength distribution with an additional higher configuration complexity $2q\otimes 2phonon$, the fragmentation effect is reinforced again. One can notice, in particular, the appearance of excited states at lower energies and the overall even more uniform strength redistribution, compared to RQTBA. Thus, the three models with the increasing complexity of the dynamical kernel form a hierarchy which translates to the hierarchy of spectral functions with increasing richness of their fine structure, as it was predicted in Ref. Litvinova (2015). Now this purely theoretical conjecture is confirmed by direct calculations.

The bottom panel of Fig. 19 displays the same strength functions as Fig. 18, but computed with a larger smearing parameter $\Delta=$ 200 keV, in order to confront them with the experimental data of Ref. Wieland et al. (2018). In contrast to the case of 48Ca, the major part of the investigated strength below 15 MeV, which is often associated with the pygmy dipole resonance and the neutron skin oscillation, lies above the neutron emission threshold that is located experimentally at $\sim$7.8 MeV nnd . Thus, this strength forms mostly a continuous spectrum. It can be seen that RQRPA provides a too poor description of the observed strength: it gives a distinct peak at about 9.5 MeV while the experimental strength of Ref. Wieland et al. (2018) shows a nearly flat distribution with a slow growth toward higher energies up to 13 MeV, where it relatively sharply increases by a factor of two. RQTBA approach produces a significant improvement of the description of the strength: the main RQRPA peak is fragmented and the overall strength distribution comes out much flatter following better the experimental trend. The EOM/RQTBA3, in turn, smoothes the strength distribution further improving the agreement with data. The only remaining drawback is that the total strength between 6 and 10 MeV is somewhat overestimated. Further refinement of the model should clarify whether more spreading toward lower energies can be induced by more complex configurations and the exact treatment of the continuum Litvinova and Tselyaev (2007); Kamerdzhiev et al. (1998) or the static effective interaction of the NL3 type, the minimal RMF parametrization with only 6 parameters, employed for these calculations is responsible for the remaining discrepancy.

In the top panel of Fig. 19 we show the ISE1 counterpart of the low-energy dipole strength in 68Ni. Remarkably, the coarse-grain pattern of the isoscalar dipole strength is very similar to that of the electromagnetic one. A similar sharp peak appears in RQRPA at about 9.5 MeV and similar fragmentation effects are induced by the 2q$\otimes$phonon and 2q$\otimes$2phonon configurations. In the final EOM/RQTBA3 calculation a relatively distinct peak at approximately 7.5 MeV remains on the background of the flat isoscalar strength distribution, that is not the case for the EME1 strength. While there is no experimental data for the ISE1 strength in 68Ni, some theoretical studies are available. In particular, Ref. Papakonstantinou et al. (2015) provides RPA and QRPA calculations of the isoscalar dipole strength for a chain of nickel isotopes including 68Ni. In the low-energy region both QRPA based on the Gogny D1S forces and continuum RPA with the SLy4 Skyrme interaction give a dominant peak around 10.5 MeV, that agrees reasonably well with our RQRPA calculation. For the EME1 strength the authors of Ref. Papakonstantinou et al. (2015) obtain a two-peak structure at the energies corresponding to the major and a minor peaks of their ISE1 strength distribution. However, fragmentation effects, if they were added beyond R(Q)RPA, would, probably, change those patterns, as it typically occurs in various implementations of the PVC mechanism. The insights about the exact continuum effects provided in this work are very important and point out to the necessity of an accurate continuum treatment.

Other types of interactions may be also considered in a future work. Density-dependent parametrizations of the meson-exchange interaction Typel and Wolter (1999); Lalazissis et al. (2005) or point-coupling Bürvenich et al. (2002); Nikšić et al. (2008) should provide a better performance in the description of the modes related to the symmetry energy as they imply more careful fits of the isovector sector Meng (2016). Ideally, the realization of the presented approach should be based on a microscopic interaction, in order to increase the predictive power. Numerical implementations based on microscopic interactions should provide a reasonable approximation to the two-body density matrix at the starting point. There can be various strategies, such as the Similarity Renormalization Group Hergert et al. (2016), Brückner G-matrix McIlroy et al. (2018); Shen et al. (2017) or the Unitary Correlation Operator Method Papakonstantinou and Roth (2009, 2010) with subsequent solution of the RPA equations and extracting the two-body densities. The capabilities of various potentials describing the nucleon-nucleon scattering data to successfully perform within the presented approach will be also addressed by future effort.

VIII Summary and outlook

In this article we revisit, compare and advance non-perturbative approaches to the quantum many-body problem. The equation of motion method is reviewed for the one-fermion and two-time two-fermion Green functions in a strongly-correlated medium. The dynamical kernels of the final EOM’s containing three- and four-body propagators are approximated by the non-perturbative cluster expansions truncated on the two-body level. The resulting EOM’s form a closed set of equations for one- and two-fermion propagators, where the latter include the particle-hole, particle-particle, and hole-hole components.

This approach is confronted with another class of closely related methods developed originally as extensions of the Landau-Migdal Fermi-liquid theory by the particle-vibration coupling and time blocking techniques, PVC-TBA. We showed that, in fact, the latter methods employed the EOM’s dynamical mass operator, whose structure can be mapped to the coupling between the single-fermion and emergent phonon degrees of freedom. These phonons are built from the correlated fermionic pairs present in the EOM’s dynamical kernels and, in the simplest non-perturbative random phase approximation, acquire a character of harmonic vibrations. To address the description of the particle-hole response the PVC-TBA method starts, in contrast to the EOM, from the general Bethe-Salpeter equation for the four-time particle-hole propagator, however, the final PVC-TBA equation is reduced to the two-time dependence by a time projection method and, thereby, to the one-energy variable equation for the spectral image. Eventually, the dynamical kernel of the resonant PVC-TBA is found to be topologically equivalent to a part of the EOM’s kernel containing terms with single two-fermion correlation functions. However, while PVC-TBA is based on the effective description of the static part of the interaction kernel (which enters, in turn, its dynamical part), the EOM method provides an accurate derivation of both the static and dynamical kernels from the single underlying bare interaction.

The insights revealed throughout this work allowed for further developing the nuclear response theory beyond the previous content and capabilities of the PVC-TBA. For this purpose, we followed the opportunities offered by the EOM, first of all, in advancing its dynamical interaction kernel beyond the $2p-2h$ level and discussed a possible iterative algorithm on the way to a highly-accurate approach to the two-fermion correlation functions. We performed a numerical implementation of the approach with the $3p-3h$ dynamical kernel and investigated the dipole response of medium-mass nuclei 42,48Ca and 68Ni. The obtained results showed some important refinements in both high- and low-energy sectors of the dipole response and indicated that the approach is indeed systematically improvable and converging. The possibility of a continued iterative algorithm opens the way to a highly-accurate nuclear many-body approach of the shell-model quality, but without prohibitive limitations on the excitation energy and mass.

An important aspect of the EOM method, which is not a feature of the PVC-TBA, is its direct connection to the bare interaction between fermions. In fact, the EOM derivation of the response theory together with the theory for the one-fermion propagator is based solely on the knowledge about this interaction. In contrast, the existing PVC-TBA and its nuclear field theory analogs imply an assumption about the existence of the underlying energy density functional, which provides information about the static part of the interaction kernel. The parameters of this functional are fitted to data for bulk nuclear properties and nuclear matter and, thus, disconnected from the bare nucleon-nucleon interaction. This feature occurs to be rather a drawback, because it lowers the predictive power of the theory. Therefore, an approach based on the bare interaction, which would also include higher configurations in the dynamical kernel in a manner discussed in this work, appears as a desirable solution. Such a theory would further clarify the mechanisms of emergent collective phenomena, superfluidity and other dynamical aspects of strongly-correlated many-body systems.

Acknowledgements

Very fruitful discussions with Hans Feldmeier, Peter Ring, Caroline Robin, and Herlik Wibowo are gratefully acknowledged. We thank Oliver Wieland for providing insights and numerical data for the experimental low-energy dipole strength in 68Ni. This work is supported by the US-NSF Career Grant PHY-1654379.