Mean-Field Theory for Fermion Pairs and the ab initio   Particle-Vibration-Coupling Approach

Peter Schuck

arXiv:1905.04672·nucl-th·January 8, 2020

Mean-Field Theory for Fermion Pairs and the ab initio Particle-Vibration-Coupling Approach

Peter Schuck

PDF

TL;DR

This paper develops a mean-field theory for fermion pairs using a Dyson Bethe-Salpeter equation with frequency-dependent kernels, coupled with a particle-vibration coupling approach, showing promising numerical results for model systems.

Contribution

It introduces a coupled Dyson-BSE and particle-vibration coupling framework for fermion pairs, integrating static mean fields with dynamic correlation effects.

Findings

01

Numerical results for Lipkin model are promising.

02

Numerical results for 1D Hubbard chain are promising.

03

The approach effectively captures particle-vibration interactions.

Abstract

A Dyson Bethe-Salpeter equation (Dyson-BSE) for fermion pairs is presented whose kernel has a static and a one frequency dependent contribution, analogous to the self energy of the single particle Dyson equation with the (static) mean field term and the energy dependent correlation term. The static part of the Dyson-BSE is the self-consistent mean field for the vibrations. At the same time, for the correlated single particle self-energy a full particle-vibration coupling (PVC) scattering equation is established where the vibration is the same as obtained from the Dyson-BSE. Both equations, single particle Dyson equation and Dyson-BSE, are coupled through self-consistency. Numerical results for Lipkin and 1D Hubbard chain are very promising.

Tables1

Table 1. Table 1: The various momentum transfers in the 6 sites case.

$q = \pm \frac{2 π}{3}$	$q = \pm π$	$q = \pm \frac{π}{3}$
$51 \to q_{51} = + \frac{2 π}{3}$	$61 \to q_{61} = - π$	$42 \to q_{42} = - \frac{π}{3}$
$63 \to q_{63} = + \frac{2 π}{3}$	$52 \to q_{52} = + π$	$53 \to q_{53} = + \frac{π}{3}$
$41 \to q_{41} = - \frac{2 π}{3}$	$43 \to q_{43} = - π$
$62 \to q_{62} = - \frac{2 π}{3}$

Equations169

G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{t - t^{'}} = - i ⟨ 0∣ T (c_{k_{1}} c_{k_{2}})_{t} (c_{k_{2}^{'}}^{+} c_{k_{1}^{'}}^{+})_{t^{'}} ∣0 ⟩

G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{t - t^{'}} = - i ⟨ 0∣ T (c_{k_{1}} c_{k_{2}})_{t} (c_{k_{2}^{'}}^{+} c_{k_{1}^{'}}^{+})_{t^{'}} ∣0 ⟩

H

H

k k^{'} \sum

\overset{v}{ˉ}_{k_{1} k_{2} k_{3} k_{4}} = ⟨ k_{1} k_{2} ∣ v ∣ k_{3} k_{4} ⟩ - ⟨ k_{1} k_{2} ∣ v ∣ k_{4} k_{3} ⟩

\overset{v}{ˉ}_{k_{1} k_{2} k_{3} k_{4}} = ⟨ k_{1} k_{2} ∣ v ∣ k_{3} k_{4} ⟩ - ⟨ k_{1} k_{2} ∣ v ∣ k_{4} k_{3} ⟩

(

(

N_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{pp} = δ_{k_{1} k_{2}, k_{1}^{'} k_{2}^{'}} N_{k_{1} k_{2}}^{pp}; N_{k_{1} k_{2}}^{pp} = 1 - n_{k_{1}} - n_{k_{2}}

N_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{pp} = δ_{k_{1} k_{2}, k_{1}^{'} k_{2}^{'}} N_{k_{1} k_{2}}^{pp}; N_{k_{1} k_{2}}^{pp} = 1 - n_{k_{1}} - n_{k_{2}}

⟨ 0∣ c_{k_{1}}^{+} c_{k_{1}^{'}} ∣0 ⟩ = δ_{k_{1} k_{1}^{'}} n_{k_{1}}

⟨ 0∣ c_{k_{1}}^{+} c_{k_{1}^{'}} ∣0 ⟩ = δ_{k_{1} k_{1}^{'}} n_{k_{1}}

\tilde{e}_{k} = e_{k} + V_{k}^{MF}

\tilde{e}_{k} = e_{k} + V_{k}^{MF}

V_{k}^{MF} = ⟨ 0∣ {c_{k}, [H, c_{k}^{+}]} ∣0 ⟩ = k^{'} \sum \overset{v}{ˉ}_{k k^{'} k k^{'}} n_{k^{'}}

V_{k}^{MF} = ⟨ 0∣ {c_{k}, [H, c_{k}^{+}]} ∣0 ⟩ = k^{'} \sum \overset{v}{ˉ}_{k k^{'} k k^{'}} n_{k^{'}}

K_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{pp}

K_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{pp}

A_{k_{1} k_{2}} = c_{k_{1}} c_{k_{2}}

A_{k_{1} k_{2}} = c_{k_{1}} c_{k_{2}}

J_{k_{1} k_{2}}^{pp} = [A_{k_{1} k_{2}}, V] = j_{k_{1}} c_{k_{2}} + c_{k_{1}} j_{k_{2}}

J_{k_{1} k_{2}}^{pp} = [A_{k_{1} k_{2}}, V] = j_{k_{1}} c_{k_{2}} + c_{k_{1}} j_{k_{2}}

j_{k} = [c_{k}, V] = \frac{1}{2} k_{2} k_{3} k_{4} \sum \overset{v}{ˉ}_{k k_{2} k_{3} k_{4}} c_{k_{2}}^{+} c_{k_{4}} c_{k_{3}} .

j_{k} = [c_{k}, V] = \frac{1}{2} k_{2} k_{3} k_{4} \sum \overset{v}{ˉ}_{k k_{2} k_{3} k_{4}} c_{k_{2}}^{+} c_{k_{4}} c_{k_{3}} .

G_{k k^{'}}^{t - t^{'}} = - i ⟨ 0∣ T c_{k} (t) c_{k^{'}}^{+} (t^{'}) ∣0 ⟩

G_{k k^{'}}^{t - t^{'}} = - i ⟨ 0∣ T c_{k} (t) c_{k^{'}}^{+} (t^{'}) ∣0 ⟩

(i \partial_{t} - e_{k}) G_{k}^{t - t^{'}} = δ (t - t^{'}) + \int d t_{1} Σ_{k}^{t - t_{1}} G_{k}^{t_{1} - t^{'}}

(i \partial_{t} - e_{k}) G_{k}^{t - t^{'}} = δ (t - t^{'}) + \int d t_{1} Σ_{k}^{t - t_{1}} G_{k}^{t_{1} - t^{'}}

Σ_{k}^{t - t^{'}} = V_{k}^{MF} δ (t - t^{'}) - i ⟨ 0∣ T j_{k} (t) j_{k}^{+} (t^{'}) ∣0 ⟩_{irr.}

Σ_{k}^{t - t^{'}} = V_{k}^{MF} δ (t - t^{'}) - i ⟨ 0∣ T j_{k} (t) j_{k}^{+} (t^{'}) ∣0 ⟩_{irr.}

K_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{pp, 0}

K_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{pp, 0}

C_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}} = ⟨ 0∣ c_{k_{1}^{'}}^{+} c_{k_{2}^{'}}^{+} c_{k_{2}} c_{k_{1}} ∣0 ⟩ - n_{k_{1}} n_{k_{2}} δ_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}

C_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}} = ⟨ 0∣ c_{k_{1}^{'}}^{+} c_{k_{2}^{'}}^{+} c_{k_{2}} c_{k_{1}} ∣0 ⟩ - n_{k_{1}} n_{k_{2}} δ_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}

G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{r - ppRPA} = G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{0, r - ppRPA} + k_{3} k_{4} \sum G_{k_{1} k_{2} k_{3} k_{4}}^{0, r - ppRPA} \overset{v}{ˉ}_{k_{3} k_{4} k_{3}^{'} k_{4}^{'}}

G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{r - ppRPA} = G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{0, r - ppRPA} + k_{3} k_{4} \sum G_{k_{1} k_{2} k_{3} k_{4}}^{0, r - ppRPA} \overset{v}{ˉ}_{k_{3} k_{4} k_{3}^{'} k_{4}^{'}}

. G_{k_{3} k_{4} k_{1}^{'} k_{2}^{'}}^{r - ppRPA}

G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{0, r - ppRPA} = \frac{1 - n _{k_{1}} - n _{k_{2}}}{ω - e _{k_{1}} - e _{k_{2}}} δ_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}

G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{0, r - ppRPA} = \frac{1 - n _{k_{1}} - n _{k_{2}}}{ω - e _{k_{1}} - e _{k_{2}}} δ_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}

G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}

G_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}

K_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{pp, dyn.} = - i ⟨ 0∣ T J_{k_{1} k_{2}} (t) J_{k_{1}^{'} k_{2}^{'}}^{+} (t^{'}) ∣0 ⟩_{irr.}

K_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{pp, dyn.} = - i ⟨ 0∣ T J_{k_{1} k_{2}} (t) J_{k_{1}^{'} k_{2}^{'}}^{+} (t^{'}) ∣0 ⟩_{irr.}

(ω - \tilde{ϵ}_{k_{1}} + \tilde{ϵ}_{k_{2}}) R_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}} (ω) = N_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{p h}

(ω - \tilde{ϵ}_{k_{1}} + \tilde{ϵ}_{k_{2}}) R_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}} (ω) = N_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{p h}

+ k_{3} k_{4} \sum [K_{k_{1} k_{2} k_{3} k_{4}}^{p h, 0} + K_{k_{1} k_{2} k_{3} k_{4}}^{p h, dyn} (ω)] N_{k_{3} k_{4}}^{- 1} R_{k_{3} k_{4} k_{1}^{'} k_{2}^{'}} (ω) .

N_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{p h} = δ_{k_{1}, k_{2}} N_{k_{1} k_{2}}^{p h} \equiv δ_{k_{1}, k_{2}} (n_{k_{2}} - n_{k_{1}})

N_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{p h} = δ_{k_{1}, k_{2}} N_{k_{1} k_{2}}^{p h} \equiv δ_{k_{1}, k_{2}} (n_{k_{2}} - n_{k_{1}})

R_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}} (t - t^{'}) = - i ⟨ 0∣ T {c_{k_{2}}^{†} (t) c_{k_{1}} (t) c_{k_{1}^{'}}^{†} (t^{'}) c_{k_{2}^{'}} (t^{'})} ∣0 ⟩

R_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}} (t - t^{'}) = - i ⟨ 0∣ T {c_{k_{2}}^{†} (t) c_{k_{1}} (t) c_{k_{1}^{'}}^{†} (t^{'}) c_{k_{2}^{'}} (t^{'})} ∣0 ⟩

K_{k_{1} k_{2} k_{3} k_{4}}^{p h, 0} = N_{k_{1} k_{2}}^{p h} \overset{v}{ˉ}_{k_{1} k_{4} k_{2} k_{3}} N_{k_{3} k_{4}}^{p h}

K_{k_{1} k_{2} k_{3} k_{4}}^{p h, 0} = N_{k_{1} k_{2}}^{p h} \overset{v}{ˉ}_{k_{1} k_{4} k_{2} k_{3}} N_{k_{3} k_{4}}^{p h}

\displaystyle\quad\Big{[}-\frac{1}{2}\sum_{ll^{\prime}l^{\prime\prime}}(\delta_{k_{2}k_{4}}\bar{v}_{k_{1}ll^{\prime}l^{\prime\prime}}C_{l^{\prime}l^{\prime\prime}k_{3}l}+\delta_{k_{1}k_{3}}\bar{v}_{ll^{\prime}k_{2}l^{\prime\prime}}C_{k_{4}l^{\prime\prime}ll^{\prime}})

+ l l^{'} \sum (\overset{v}{ˉ}_{k_{1} l k_{3} l^{'}} C_{k_{4} l^{'} k_{2} l} + \overset{v}{ˉ}_{k_{4} l k_{2} l^{'}} C_{k_{1} l^{'} k_{3} l})

\displaystyle\quad-\frac{1}{2}\sum_{ll^{\prime}}(\bar{v}_{k_{1}k_{4}ll^{\prime}}C_{ll^{\prime}k_{2}k_{3}}+\bar{v}_{ll^{\prime}k_{2}k_{3}}C_{k_{1}k_{4}ll^{\prime}})\Big{]},

K_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{p h, dyn} (t - t^{'})

K_{k_{1} k_{2} k_{1}^{'} k_{2}^{'}}^{p h, dyn} (t - t^{'})

J_{k_{1} k_{2}}^{p h}

J_{k_{1} k_{2}}^{p h}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

11institutetext: Institut de Physique Nucléaire Orsay 22institutetext: Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France

Mean-Field Theory for Fermion Pairs and the ab initio Particle-Vibration-Coupling Approach

Peter Schuck Present address: Institut de Physique Nucléaire d’Orsay, Université Paris-Sud, CNRS–IN2P3. 15, Rue Georges Clemenceau, 91406 Orsay Cedex, France.11

(Received: date / Revised version: date)

Abstract

A Dyson Bethe-Salpeter equation (Dyson-BSE) for fermion pairs is presented whose kernel has a static and a one frequency dependent contribution, analogous to the self energy of the single particle Dyson equation with the (static) mean field term and the energy dependent correlation term. The static part of the Dyson-BSE is the self-consistent mean field for the vibrations. At the same time, for the correlated single particle self-energy a full particle-vibration coupling (PVC) scattering equation is established where the vibration is the same as obtained from the Dyson-BSE. Both equations, single particle Dyson equation and Dyson-BSE, are coupled through self-consistency. Numerical results for Lipkin and 1D Hubbard chain are very promising.

pacs:

PACS-keydiscribing text of that key and PACS-keydiscribing text of that key

1 Introduction

Nuclei are a paradigmatical case where single nucleon properties, response functions, and pair propagators are strongly influenced by collective modes. In order to cope with this situation, the so-called particle-vibration coupling (PVC) approach has been put forward by P. F. Bortignon and collaborators Colo . The PVC approach is based on effective nucleon forces of the Skyrme or Gogny type in the non-relativstic scheme or on effective density dependent coupling constants in the relativistic theory. In this sense PVC is still largely phenomenological.

The objective of this work is to present a PVC approach which can be used in ab initio calculations where one starts from a given bare force. Such systems are, e.g., electron systems with the Coulomb force and the coupling of single particle (s.p.) motion to the plasmon resonance but one can also think to start for nuclear systems with a ’bare’ force which reproduces very well nucleon-nucleon phase shifts as does for example the Bonn force CD-Bonn . In order to achieve this, we will make use of a mean field approach for correlated fermion pair operators of the type $(a^{+}a)$ or $(aa)$ and $(a^{+}a^{+})$ with $a^{+},a$ fermion creation and destruction operators. We will show that the corresponding correlated ground state wave function is given to good approximation by the so-called Coupled-Cluster-Doubles (CCD) wave function jemai13 . The collective modes obtained in this way can be used to set up a particle-vibration scattering equation entering the self-energy of the single fermion Dyson equation, so that via the s.p. occuaption numbers the single particle and two particle propagators get selfconsistently coupled. The performance of this approach will be demonstrated with model cases where a comparison with exact solutions is possible.

2 The formalism

The mean-field approach for fermion pairs is, in principle, not new and we will recapitulate it here only very succinctly. We will present it, establishing a formally exact Bethe-Salpeter equation (BSE) whith a one frequency kernel which has a definite expression. The point we want to make here is that in others than the nuclear field a one frequency BSE is very uncommon and mostly propagators with four times (or three times after translational invariance in time is taken into account) are considered implying that also the integral kernel depends on three frequencies. However, since this leads to numerically inextractable equations, mostly the three time (or three frequency) kernel is taylored down either to a static one or to a one frequency kernel making use of some suitable approximations. We have recently published a longer article where this problematic is discussed at length in the particle-hole (ph) channel in the context of electronic systems Julien . In this work we will discuss, additionally, the particle-particle (pp) channel with its possible pairing instability. However, we also will expose specific aspects in the ph channel and in any case pp and ph channels are coupled, e.g. via screening of the bare interaction in both channels. As mentioned in the Introduction, a novelty will also be that we set up a single particle self-energy which is consistent with the mean field approach of fermion pairs. Notably a particle-vibration scattering equation will be established.

Let us start with what we want to call a formerly exact Dyson-Bethe-Salpeter equation (Dyson-BSE) for the two time pair propagator

[TABLE]

Here T is the time ordering operator, and $|0\rangle$ stands for the exact ground state and the fermion pair operators turn in time with the two body Hamiltoinan

[TABLE]

where $e_{kk^{\prime}}$ is the single particle matrix comprising kinetic energy and external potential, and

[TABLE]

is the antisymmetrized matrix element of the two body interaction.

The two-body propagator obeys the following exact Dyson-BSE Julien ; NPA628

[TABLE]

where with $\delta_{k_{1}k_{2},k^{\prime}_{1}k^{\prime}_{2}}=\delta_{k_{1}k^{\prime}_{1}}\delta_{k_{2}k^{\prime}_{2}}-\delta_{k_{1}k^{\prime}_{2}}\delta_{k_{2}k^{\prime}_{1}}$

[TABLE]

and we supposed that we work in the canonical basis where the density matrix is diagonal, that is

[TABLE]

Furthermore the s.p. energies in (4) are given by

[TABLE]

where the mean field shift is included

[TABLE]

and where $\{..\}$ stands for the anticommutator. We assumed that mean-field energies and density matrix can be diagonalized simultaneously. The integral kernel is given by

[TABLE]

where we abreviated

[TABLE]

and

[TABLE]

with

[TABLE]

Please note that the $K$ matrix which after Fourier transform depends only on one frequency can be interpreted as a self-energy for the motion of a fermion pair. As usual the self energy is split into a frequency independent, static part and a truely frequency dependent, dynamic part. The latter must be two line irreducible, hence the index ’irr.’. This partition is in complete analogy to the case of the s.p. self-energy appearing in the s.p. Dyson equation for the s.p. Green’s function

[TABLE]

The Dyson equation then reads

[TABLE]

with the self-energy expressed as

[TABLE]

where the index ’irr.’ again indicates that the corresponding correlation function should be one-line irreducible.

It is well known that the self-energy of the s.p. motion is closely related to the optical potential of elastic nucleon-nucleon scattering Bell-Squ . Now suppose one wants to consider elastic deuteron-nucleus scattering (or any elastic scattering of composite bosons) how would one define a corresponding optical potential which by definition will depend only on the incoming deuteron energy ? Of course quite naturally the deuteron optical potential can directly be derived from the one frequency $K$ -matrix of the Dyson-BSE defined above, see an early work on this in Ref. opt-pot .

2.1 Static part of the BSE kernel

Let us now discuss the $K^{pp,0}$ term of the BSE kernel. To establish an explicit form for $K^{pp.0}$ , we have to evaluate the double commutator contained in the pair mean field part of $K^{pp,0}$ NPA628

[TABLE]

where

[TABLE]

which is the fully correlated two body or cumulant form of the density matrix.

We see that $K^{pp,0}$ involves, besides occupation numbers, static two-body correlation functions. They are of two types: there are single line corrections with one of the two s.p. motions uneffected by the correlations (those with the Kronecker symbols) and there are exchange terms where a two-body correlation is exchanged between the two particles. Since our starting point is a two-body propagator, a self-consistent scheme can be established. This is similar to the selfconsistency involved with the s.p. mean-field, only here, naturally, two-body correlation functions have to be iterated rather than s.p. densities in the case of the s.p. mean-field. We, therefore, call $K^{pp,0}$ the ’fermion pair mean field’. Of course, there appear also s.p. densities in $K^{pp,0}$ and we will later show how they can be consistently obtained from the s.p. Green’s function.

A closer investigation of the exchange kernel, however, shows that the exchange is rather of the ph type. At least to lowest order, that is to second order, the exchange is given by a static ph exchange. This ph exchange is well known that it screens the pairing force by almost a factor of two as has first been evaluated by Gorkov, Melik-Barkhudarov (GMB) GMB ; Strinati . Actually for systems like the nuclear ones where there are more than two species of fermions (that is, four), the screening can also become anti-screening Pethick ; Urban . Also GMB did not use strict second order but replaced the vertices by the scattering length, that is the vertices have been dressed to $T$ -matrices in the low energy limit, that is, by the scattering length. Though this resummation can be understood easily by graphical analysis, how this can be derived more analytically will be discussed below in sect. II.B.

The first term on the r.h.s. of (15) is the usual two body matrix element of the interaction modified with correlated occupation numbers ( the standard particle-particle RPA as described in RS uses HF occupation numbers). One can presum this term what leads to the so-called renormalized pp-RPA Hirsch

[TABLE]

with

[TABLE]

The Dyson-BSE can then be written in the following way

[TABLE]

2.2 Dynamic part of the BSE kernel

Let us now discuss the time-dependent, dynamic part $K^{\rm dyn.}$ of the interaction kernel

[TABLE]

From (10) we see that this expression involves four different contributions: two contributions contain the two interaction vertices on the same line and two contributions on opposite lines. The latter, therefore, contain exchange processes while the former are responsible for s.p. self energy corrections. Approximating the the 3p-1h propagator involved in (20) by a product of a hole propagator and the three body propagator in second order $T$ -matrix approximation, we give a schematic grapphical representation of the term in Fig.1. This illustrates how one can replace in the second order screening term discussed above, the bare vertices by ladder $T$ -matrices and then eventually by the scattering lengths as done in GMB . We see that the exchange contributions are of the screening (or anti-screening) type whereas the other two contributions renormalize the s.p. by particle-vibration couplings. Of course, in general, all four lines are correlated.

As a matter of fact the 3p-1h propagator in (20) lends itself to several “natural” approximations other than the one we just discussed. For example instead of considering an uncorrelated ph-propagator exchange, one could take the ph-response function which will be presented in sect. II.C Lombardo ; Urban .

2.3 The ph-channel

Since we see that ph and pp channels are coupled, we immediately also give the Dyson-BSE in the ph-channel Julien

[TABLE]

With

[TABLE]

and the two-time response function defined by (with $k_{1}\neq k_{2}$ and $k_{1}^{\prime}\neq k_{2}^{\prime}$ )

[TABLE]

The static part of the integral kernel is given by

[TABLE]

and the dynamic part

[TABLE]

with

[TABLE]

As in the pp-channel, we can introduce a renormalized ph-propagator

[TABLE]

with

[TABLE]

and where we introduced the ladder $T$ -matrix

[TABLE]

which is important when dealing with systems with a hard-core potential.

Keeping from the ph-kernel only the remaining instantaneous part, one arrives at a self-consistent mean-field equation for the ph-propagation

[TABLE]

where $\tilde{K}^{ph,0}$ is the part of (24) where only the second, third, fourth, and fifth terms are kept.

The eigenvalue form of this equation, given below, is known as the self-consistent RPA (SCRPA) equation jemai13 ; NPA628 . A graphical representation of Eq.(30) is given in Fig.2

We now want to give the spetral representation of the response function, since it may be important for the following when we discuss the particle-vibration coupling (PVC).

[TABLE]

On the pole we get with (31) the following eigenvalue equation (SCRPA equation jemai13 )

[TABLE]

where

[TABLE]

Please let us note that the same equations can be derived from a different perspective. Considering the average excitation energy given by the normalised energy weighted sum rule

[TABLE]

with

[TABLE]

and its inversion

[TABLE]

With the excitation creation operator (35) we suppose that the excited state is given by $Q^{{\dagger}}_{\mu}|0\rangle=|\mu\rangle$ and that there exists the so-called killing condition

[TABLE]

Minimizing the sum rule with respect to the $X,Y$ amplitudes, again leads to the eigenvalue problem (32).

Evaluating the double commutators with the help of (36,37), we find that $A,B=A,B[X,Y,n_{k}]$ , that is, the SCRPA matrix is a functional of the RPA amplitudes $X,Y$ and the occupation probabilities $n_{k}$ where we made the approximation that $\langle\hat{n}_{k}\hat{n}_{k^{\prime}}\rangle\simeq n_{k}n_{k^{\prime}}$ what usually is a very good approximation. This approximation is not necessary but avoiding it complicates the formalism quite a bit. We again can see that the ph-channel is coupled via the exchange to the pp-channel.

Further details of the procedure and approximation schemes are presented in Julien ; jemai13 . For the moment let us come to the important question how to incorporate a consistent evaluation of the s.p. occupation numbers $n_{k}$ , a question which has remained open so far.

3 Determination of the Single Particle Occupation Numbers

The single particle occupation numbers being a s.p. quantity, it is natural to determine them via the s.p. propagator (12). The point will be to find an approximation of the self energy which is consistent with the ground state implicitly used in our two fermion mean field equations to which we want to restrict ourselves from now on. As was already remarked earlier, the ground state corresponding to the fermion pair mean field equations is to good approximation given by the CCD ground state wave function

[TABLE]

As has been shown in jemai13 , this state is to good approximation the vacuum to the creation operators where this time not all indices are allowed but only those lying in the ph configuration space.

[TABLE]

and the corresponding inversion (36).

These matters have amply been discussed in earlier literature where also the relation between the $z$ and $X,Y$ amplitudes is given jemai13 , with references in there. Notably it has been shown how the amplitudes $X,Y$ are solutions of the self-consistent RPA (SCRPA) equations (32) which we here also call fermion ph-pair mean field equations. Since this is well documented, we will not further comment on this point.

We will, however, elaborate a s.p. selfenergy entering the Dyson equation (13) which is consistent with the SCRPA solution. A self-consistent system of equations can then be established for the calculation of the occupation numbers and the amplitudes $X,Y$ , see Mohsen . For this, let us consider the dynamic part of the s.p. mass operator given in (14). This contains the following 3-body propagator

[TABLE]

We want to establish an integral equation for this propagator. As usual, we employ the equation of motion (EOM) and approximate the integral kernel by the static part. We obtain

[TABLE]

where

[TABLE]

and

[TABLE]

In order to be consistent with our ground state wave function (38), we have to restrict the indices of $K^{0}$ . One namely can show that if one wants to calculate the hole occupation numbers, then the following destruction operator kills $|Z\rangle$ of (38) exactly o-RPA

[TABLE]

where the amplitudes are related by

[TABLE]

and the amplitudes $x,U$ are obtained from the eigenvalue equation corresponding to (41) together with (13)

[TABLE]

For the explicit form of the matrix elements, see sect. V.A below.

So the hole destructor couples to a $2p-1h$ configuration. Inversly the particle destructor couples to a 2h-1p configuration. We will only discuss the former, the latter being treated analogously. In the double commutator of (43), each of the two triples of 2p-1h fermion operators contracts a particle state to the interaction. Naturally from each triple then only a ph pair remains. We can express those ph pairs via the $Q^{{\dagger}},Q$ operators of (35) using the inverse relation (36). Commuting then the destructors to the right, we exploit the killing property (37) and then $K^{0}$ is expressed by (self-consistent) RPA amplitudes $X,Y$ and occupation numbers. A graphical representation of this PVC vertex is shown in Fig.3.

In turn the occupation numbers are directly related to the s.p. Green’s function (11) and then via the dynamical part of the s.p. self energy which is related to the solution of (41), we have a closed system of equations for the SCRPA amplitudes, via the SCRPA equations, and the s.p. occupation probabilities. We want to call this sytem of equations the eo-SCRPA. It has been solved for the Lipkin and 1D Hubbard model with very good success as can be seen in the figures presented in sect.V. From the s.p. Green’s function we also can calculate the ground state energy in the usual way

[TABLE]

Of course, the solution of the $2p-1h(2h-1p)$ equations will in general be quite demanding because of the eventually large configuration space. However, presently in nuclear physics quite routinly in the so-called second RPA huge configuration $2p-2h$ spaces are considered, so that a $2p-1h$ space should be a less important problem. We should also point out that the $2p-1h$ integral equation (41) can be interpreted as a particle-vibration scattering eqation with full respect of the Pauli principle. This is schematically shown in Fig.4. The vibrations (wiggly lines) here are the solutions of the SCRPA equation.

4 Conservation laws and Ward Identities

4.1 Goldstone modes

It is well established that standard HF-RPA (BCS-QRPA) approach exhibits a Goldstone (zero) mode if the HF solution corresponds to a continuously broken symmetry. For finite systems, one mostly talks about a zero or spurious mode. For nuclei and other selfbound systems like 3He droplets, HF always breaks translational invariance and the corresponding RPA shows a zero mode RS which corresponds to a coherent displacement of the whole system. In trapped cold atom gases there exists the so-called Kohn mode where the atom cloud oscillates coherently in an external harmonic trapping potential without internal excitations of the wave packet Tohy-symm . Within BCS-QRPA, one obtains in infinite matter, because of the broken particle number $U(1)$ symmetry, a Goldstone mode, the so-called Bogoliubov-Anderson mode Combescot . Also in finite superfluid nuclei zero modes appear Rabhi .

As mentioned, these Goldstone modes reflect basic principles of quantum mechanics and it is very important not to destroy this property in theories which go beyond the HF-RPA scheme. As we see from (35), the crucial point is that the $Q^{{\dagger}}$ operator can represent the symmetry operator (let us call it $\hat{S}$ ) in question as a paticular solution of the SCRPA eqauations and that the relation $[H,\hat{S}]=0$ is not destroyed in the course of applying the formalism. In SCRPA all components of the symmetry operator, besides the diagonal ones, are present and it is important to keep them all Tohy-symm ; Delion-Schuck-Tohy ; DSD05 ; SchuToh ; tddm-scrpa . Therefore, one may think that if the diagonal matrix elements of the symmetry operator $\hat{S}_{kk}=0$ , then in any case the Goldstone mode will be present in the solution. Many symmetry operators have this property. This is the case for the linear total momentum operator because of its odd parity. In deformed nuclei, the rotational symmetry is broken. The angular momentum operator has no diagonal elements either because it is not time reversal invariant.

More subtle is the question of pairing which is one of the broken symmetries often encountered in Fermi systems. In this case the symmetry operator is the particle number operator which contains a Hermitian diagonal piece which cannot be included into the (quasi-particle) RPA operator $Q^{{\dagger}}$ as discussed already. However, in infinite systems the Bogoliubov-Anderson mode comes nevertheless because the diagonal piece of the number operator has zero weight. For finite systems like nuclei one eventually has to extend the theory to cope with the problem, see Ref. jemai13 .

Therefore in practically all situations SCRPA keeps with the Goldstone property. We want to underline the importance of this statement because it is extremely rare to find beyond standard RPA approximations which satisfy the Goldstone theorem besides the multi time (energy) approach for Green’s functions discussed by Hedin Hedin which, however, is numerically untractable beyond standard RPA.

4.2 Sum rules

We show that the energy weighted sum rule, given by the following identity

[TABLE]

is fulfilled within SCRPA. Here $|\nu\rangle$ is a complete set of eigenstates and $F$ is a one body operator

[TABLE]

One can show that this identity is automatically fulfilled if one considers that $|\nu\rangle$ is the set of SCRPA or r-RPA eigenstates. By using the inverse transformation of the fermion pair operators $c^{{\dagger}}_{k}c_{k^{\prime}}$ , one obtains

[TABLE]

with $M_{kk^{\prime}}=n_{k^{\prime}}-n_{k}$ .

Using the general system of RPA equations with excitation energies $\Omega_{\nu}=E_{\nu}-E_{0}$ , one can eliminate the amplitudes $X,Y$ in favor of the RPA matrices $A,B$ . From there it is only a small step to show that the sum-rule (48) is fullilled Delion-Schuck-Tohy .

4.3 Gauge invariance

Another important property of standard RPA which is fullfilled by SCRPA is gauge invariance (or Ward identity). Gauge invariance of standard RPA is, e.g., nicely demonstrated by Feldman and Fulton FF . The extra terms containing the two body correlation functions in (29) cancel in the limit where the two open legs are put on the same spot in position space. Actually, gauge invariance of standard RPA *as well as * SCRPA can easily be verified from (33), (35). If in these equations the operator $\delta Q^{+}_{k_{1}k_{2}}=c^{{\dagger}}_{k_{1}}c_{k_{2}}/\sqrt{n_{k_{2}}-n_{k_{1}}}$ is transformed into $r$ -space and the diagonal element is taken, as demanded to verify gauge invariance (see FF , eq, (3.69)), we immediately realize that this diagonal operator commutes with the remainder (also written in $r$ -space, once the Hamiltonian $H$ is replaced by its interaction part $V$ , that is, the Coulomb interaction. Therefore, gauge invariance is fulfilled. This argument is valid discarding spin but, as shown in FF , this does not invalidate the general proof. These considerations also entail that the so-called “velocity-length” equivalence in the dipole transition is preserved Delion-Schuck-Tohy . Please note that all properties mentioned in III.A,B,C are also fullfilled with the renormalized RPA of (27), see Schaefer ; SchuToh ; tddm-scrpa ; Rowe68 ; Catara .

5 Applications

5.1 The Lipkin Model

The Lipkin model, see, e.g., RS , is one of the most frequently used models in nuclear physics to test approximation schemes.

The single-particle space of the Lipkin model consists of two fermion levels, each of which has a N-fold degeneracy RS . The upper (lower) level has the energy of $\frac{e}{2}$ ( $-\frac{e}{2}$ ). The Hamiltonian of the Lipkin model is given by

[TABLE]

with $e$ is the inter-shell spacing, $V$ is the coupling constant, and

[TABLE]

with $2J_{0}=\hat{n}_{1}-\hat{n}_{0}$ , $\hat{n}_{i}=\sum c^{\dagger}_{im}c_{im}$ and $N$ is the number of particles equivalent to the degeneracies of the shells. We consider the odd excitation operator as (44), that is

[TABLE]

with the eigenvalue equation of (46). Based on the solution of the SCRPA equations jemai13 with the definition of the ph-pair excitation operator as $Q^{+}=(XJ_{+}-YJ_{-})/d_{0}$ (with $d_{0}=\sqrt{\langle-2J_{0}\rangle}$ and $J_{+}=(XQ^{+}-YQ)d_{0}$ ), we obtain the $X,~{}Y$ amplitudes as being the solutions of SCRPA equations (32). The matrix elements of (46) are then given by

[TABLE]

and the corresponding secular equation

[TABLE]

The occupations numbers are given by

[TABLE]

where $\lambda_{\pm}$ are the eigenvalues of the $2\times 2$ matrix problem,

[TABLE]

with $\beta=e-VXY(N-4)-VXY(N-4)(1+2Y^{2})\frac{\langle J_{0}\rangle}{n_{11}}$ . Thus,

[TABLE]

In the above equations (55) the correlation functions are expressed by the RPA amplitudes $X,Y$ in the way it is described in section III. The correlation functions which contain quadratic forms of occupation number operators as $\langle J_{0}J_{0}\rangle$ in above equation can in principle be expressed by the RPA amplitudes as well but leading to heavier expressions. Usually, we, therefore will employ the factorisation approximation leading in the present case to $\langle J_{0}J_{0}\rangle\simeq\langle J_{0}\rangle^{2}$ what mostly turns out to be quite satisfactory. However, in the case of the Lipkin model one also can use the Casimir relation

[TABLE]

to close the systme of equations. The results are shown in Figs.5 and in Mohsen . They concern in the order: i) the expectation value $\langle J_{0}\rangle$ of the difference of populations in upper and lower level, ii) the first excitation energy of the even system, iii) the excitation energies of the odd system. The correlation energy (not shown) is reproduced with the same quality. All quantities are very well reproduced throughout couplings up to the critical value $\chi=\chi_{\rm crit.}$ where the standard RPA breaks down and the system wants to change to the ’deformed’ basis. However, even values slightly beyond $\chi_{\rm crit.}=1$ are still quite acceptable. All quantities for $N=2$ are reproduced exactly. By some lucky accident the occupancies even for $N=4$ come out to be exact (as shown in Fig.5, Fig.8 and in Mohsen ).

5.2 The Hubbard Model

The Hubbard model is widely used to deal with the physics of strongly correlated electrons. Since the model can be solved exactly in one dimension (1D) and for small cluster sizes, it is very useful for theoretical investigations jemai05 ; Hubbard . To be precise, our ”Hubbard model” is a 6-site system at half filling with periodic boundary condition, described by the usual Hamiltonian jemai05 ; Hubbard :

[TABLE]

Here, $\hat{n}_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ , $c_{i\sigma}^{\dagger}$ and $c_{i\sigma}$ are the creation and annihilation operators for an electron at site $i$ with spin $\sigma$ , $U$ is the on-site (spin-independent) interaction, $-t$ is the hopping term of the kinetic energy. The eigenstates of the system will be linear combinations of Slater determinants, which are denoted by the kets $|1\ldots 6\rangle$ , with occupations of the sites $1\ldots 6$ . The hamiltonian is rewritten in plane wave basis,

[TABLE]

with the transformation

[TABLE]

where $\hat{n}_{\bf{k},\sigma}=a^{\dagger}_{\bf{k},\sigma}a_{\bf{k},\sigma}$ , $\varepsilon_{\bf{k}}=-2t\cos\left(ka\right)$ , which are, respectively, the number operator of particles of the mode $({\bf{k}},\,\sigma)$ and the energies of one particle on a lattice with a the parameter of the lattice which is taken as $a=1$ . For a problem with $N$ sites, the condition of periodicity is given by $c_{N+1,\sigma}=c_{1,\sigma}$ . This implies that $e^{-ik\,N}=1$ , hence the values taken by $k$ will be $k=\frac{2\,\pi}{N}\,n$ . In addition, the first Brillouin zone is defined on the field where $-\pi\leq k<\pi$ , which gives us the values of $n$ as $\frac{-N}{2}\leq n<\frac{N}{2}$ .

For the six sites, we have the possible states with the following wave vectors:

[TABLE]

and with the kinetic energies (see Fig.9), respectively,

[TABLE]

The transfer wave vector( $q_{ph}=k_{p}-k_{h}$ ) takes the possible values as shown in the Table 1.

At this point we proceed exactly as in the case of the Lipkin model: The excitation operator for the even system is given by

[TABLE]

with $~{}K^{\pm}_{ph\sigma}=J^{\pm}_{ph\sigma}/\sqrt{N_{ph\sigma}}$ , $~{}J^{+}_{ph\sigma}=a^{{\dagger}}_{p\sigma}a_{h\sigma}$ , $~{}N_{ph\sigma}=n_{h\sigma}-n_{p\sigma}$ . With the inversion

[TABLE]

we can calculate the mean values needed for the matrix elements of the SCRPA equations for the even particle number case

[TABLE]

and similar expressions for other expectation values of this type where we replaced the ”ph” operators by the RPA creation and destruction operators from the inversion (36) and then commute the $Q$ operators to the right until they kill the ground state. All matrices become functionals of the occupancies $n_{h}$ and $n_{p}$ and $X,Y$ amplitudes in analogy to what was the case in the Lipkin model and, thus, the diagonalisation process implies at the same time an iteration on the occupancies and the amplitudes.

For the odd particle number case, we make again the following ansatz

[TABLE]

From there, we can, as outlined in the general section II, and as just now for the case of the Lipkin model, calculate the occupation numbers. The results for the occupation numbers are again excellent, see Fig. 10. Also the excitation energies of the even particle number system, see Fig. 11 are very well reproduced.

6 Conclusion

In this work, we coupled even and odd particle numbers RPA selfconsistently. Both systems are based on the same correlated RPA ground state. From the odd system, we get the occupation numbers, odd particle excitation energies, and the ground state energies whereas from the even SCRPA equations we get the excitation energies of the even system and transition probablities. Both even and odd systems are coupled through non-linear equations. We called this system of equations ’even-odd’ SCRPA (eo-SCRPA). Applications to the Lipkin model and a six sites Hubbard ring at half filling gave excellent results for all quantities. The equations are relatively complex due to their non-linearity but they should be solvable with modern computers for realistic problems. The coupling of even and odd RPA’s has a couple of advantages: it gives richer results, i.e., excitation energies of even and odd particle number systems; there is a natural way how to obtain the ground state energy via the s.p. Green’s function and, last but not least, all qualities of standard RPA, as there are, the Goldstone theorem, sum-rules, and gauge invariance, respectively Ward identities are maintained. Since there is no phenomenological input in the eo-SCRPA equations and the hard core is tamed by the $T$ -matrix, our scheme can be qualified as an “ab initio” PVC approach. The results of the Lipkin and Hubbard models seem to be very promising. Work for realistic applictions is planned for the future Litvinova .

7 Acknowledgements

I am very greatful to M. Jemai and Zhou Bo for their help with the preparation of the figures and M. Jemai for collaboration on this subject. I appreciate longstanding collaboration on SCRPA with D. Delion, J. Dukelsky, and M. Tohyama. Thanks are also due to G. Röpke for discussions concerning the pp-channel. V. Olevano and J. Toulouse are collaborators of a very fruitful recent work and publication in the ph-channel Julien .

Appendix A Equation of Motion for odd particle number operator for Hubbard Model

For the Hubbard model (60) we define the odd excitation operator as in Eq.(44),

[TABLE]

with $J^{+}_{ph-}=a^{\dagger}_{p-}a_{h-}$ and $\sigma=\uparrow,\downarrow=+,-$ . Remembering the notations for the occupation probabilities

[TABLE]

we have $n_{k_{2}\sigma}=n_{k_{3}\sigma}$ , $n_{k_{4}\sigma}=n_{k_{5}\sigma}$ , $n_{k_{2}\sigma}=1-n_{k_{3}\sigma}$ and $n_{k_{1}\sigma}=1-n_{k_{6}\sigma}$ . This gives

[TABLE]

The term without interaction $H_{0}=\sum_{k\sigma}\varepsilon_{k}\hat{n}_{k\sigma}$ is given by

[TABLE]

with ${{\cal N}}_{p^{\prime}ph}=\langle(1-\hat{n}_{p^{\prime}+})(-2J^{0}_{ph,-})\rangle+\langle J^{+}_{ph,-}J^{-}_{hp,-}\rangle$ . The term in the Hamiltonian for the transfer $q=0$ , $H_{q=0}=\frac{U}{6}\sum_{kk^{\prime}}\hat{n}_{k+}\hat{n}_{k^{\prime}-}$ leads to

[TABLE]

with $\sum_{k}\hat{n}_{k\sigma}=\sum_{p}\hat{n}_{p\sigma}+\sum_{h}\hat{n}_{h\sigma}=3$ in the half-filled case. Now let us calculate the elements ${{\cal C}}_{p^{\prime}ph}$ for the first row (or column) as

[TABLE]

The elements of the matrix except the first row (or column) are given as follows

[TABLE]

In the following, as already discussed several times, we retain from (74) only those terms where the particle states of the left and right triple operators in ${{\cal D}}$ connect to the interaction. The remaining density operator from the interaction is approximated by its diagonal form. This leads to expressions evaluated in (75) below. First let us discuss what kind of terms we are neglecting in this way. It should be noted that the terms of type $\langle J^{\pm}_{ph}J^{\pm}_{p^{\prime}h^{\prime}}J^{\pm}_{p^{\prime\prime}h^{\prime\prime}}\rangle=0$ , $\langle J^{\pm}_{ph}S_{p_{1}p_{2}}J^{\pm}_{p^{\prime\prime}h^{\prime\prime}}\rangle$ are probably small (with $S_{p_{1}p_{2}}=a^{\dagger}_{p_{1}}a_{p_{2}}$ for $p_{1}\neq p_{2}$ ) and $\langle J^{\pm}_{ph}S_{h_{1}h_{2}}J^{\pm}_{p^{\prime\prime}h^{\prime\prime}}\rangle$ also small (with $S_{h_{1}h_{2}}=a^{\dagger}_{h_{1}}a_{h_{2}}$ for $h_{1}\neq h_{2}$ ) in eq.(74). As shown in jemai13 , the term $\langle SJ\rangle=0$ and $\langle SS\rangle$ are small. Only the terms nonzero in eq.(74) like $\langle J^{\pm}_{ph+}n_{k\pm}J^{\pm}_{p^{\prime}h^{\prime}-}\rangle$ which can be calculated as shown in (75). With the short hand notation $ph\sigma\equiv i$ , $k\sigma\equiv k$ , $\hat{N}_{i}=\hat{n}_{h\sigma}-\hat{n}_{p\sigma}$ and $N_{i}=n_{h\sigma}-n_{p\sigma}$ , we can evaluate the following terms

[TABLE]

and other mean values of similar type.

Bibliography31

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) G. Colo, P. F. Bortignon, H. Sagawa, K. Moghrabi, M. Grasso, Nguyen Van Giai, J. Phys.: Conf. Ser. 321, 012018 (2011); G. Colo, Yifei Niu, E. Vigezzi, P. F. Bortignon, EPJ Web of Conferences 107 , 06001 (2016); G. Colo, H. Sagawa, P. F. Bortignon, Phys. Rev. C 82 , 064307 (2010); G. Colo, P. F. Bortignon, Progr. Theor. Phys. Supplement 146 : 543 (2002).
2(2) R. Malfliet, Phys. Rev. C 63 , 024001 (2001).
3(3) M. Jemai, D. S. Delion and P. Schuck, Phys. Rev. C 88, 044004 (2013).
4(4) V. Olevano, J. Toulouse, P. Schuck, J. Chem. Phys. 150 , 084112 (2019).
5(5) J. Dukelsky, G. Röpke, P. Schuck, Nucl. Phys. A 628 , 17 (1998).
6(6) J. S. Bell, J. Squires, Phys. Rev. Lett. 3, 96 (1959).
7(7) P. Schuck, Z. Physik 241 , 395 (1971).
8(8) L. P. Gor’kov, T. K. Melik-Barkhudarov, Sov. Phys.-JETP 13, 1018 (1961).