text:
# Pair Fluctuations in Ultra-small Fermi Systems within Self-Consistent RPA at Finite Temperature

## I Introduction

Pairing properties of finite Fermi systems such as ultrasmall metallic grains have recently received a great deal of attention. This has been spurred by a series of spectacular experiments of Ralph, Black and Tinkham . In order to correctly describe pairing properties it has been recognized that the finiteness of the systems (grains) needs to consider quantum fluctuations, good particle number, number parity, etc. seriously, since the coherence length may be of the order of the system size. The situation for metallic grains has in the meanwhile been well described in several review articles  (see also ). Another system where the finiteness is at the forefront of the theoretical investigation since several decades is the superfluid atomic nucleus. As a matter of fact many of the theoretical tools such as particle number projection, even-odd effects, number parity, blocking effect, particle - particle Random Phase Approximation (pp-RPA), etc. have first been developed in nuclear physics  before finding their application to finite systems of condensed matter. Also the schematic pairing model with which we will mostly deal in this paper, namely the Picket Fence Model (PFM), whose exact solution has been found by Richardson and Sherman , has essentially been developed in the context of nuclear physics for the description of deformed superfluid nuclei. For finite condensed matter systems an early theoretical description was proposed by Mühlschlegel, Scalapino, and Denton  using the Static Path Approximation (SPA) to the partition function. This work stayed rather singular for a long time but the SPA has recently been applied successfully to the PFM both in the condensed matter  and nuclear  contexts. A further standard method to treat quantum fluctuations namely the well known RPA has quite extensively been used for nuclear systems  but equally for condensed matter problems .

In this work we will further elaborate on the RPA approach. We indeed have recently had quite remarkable success with a self consistent extension of the pp-RPA, which we called Self-Consistent RPA (SCRPA), by reproducing very accurately groundstate and excitation energies of the PFM  at zero temperature. This formalism was also developed independently by Röpke and collaborators who called it Cluster - Hartree - Fock (CHF) . Such type of generalization of the RPA theory grew out of the works of K.-J. Hara and D. Rowe  several decades ago. Shortly afterwards the theory was rederived using the method of many body Green functions . The success of the theory motivates us to develop the SCRPA formalism also for the finite temperature case and to study the thermodynamic properties of the BCS Hamiltonian using the PFM as an example. For the extension of SCRPA to finite temperature we use the Matsubara Green functions approach . It appeared that the approximation scheme is very effective in treating two-body correlations in the particle-particle (pp) channel as well as the Pauli principle effects. We should mention that we will work with real particles and not with quasiparticles what should limit our approach to temperatures above the critical temperature $`T_c`$ (i.e. to the normal phase). However, as we will see below, the definition of $`T_c`$ in SCRPA is not so clear and we will be able to continue our calculation quite deeply into the superfluid regime.

We organize the paper in the following way. In Section 2, the approach is outlined in general. Then, in Section 3, the formalism is applied to the PFM. A comparison with the exact solutions as well as with the results of other approximations is made in Section 4. Section 5 is devoted to comparison with other recent works. In section 6, we will summarize the results and draw some conclusions. In an appendix a variant for the calculation of the occupation numbers is proposed.

## II General formalism

In treating a finite many-body system at finite temperature, it is convenient to use the grand canonical ensemble although it violates the number conservation. With the definition

$`K=H\mu N`$

the grand partition function and statistical operator read

$`Z_G=e^{\beta \mathrm{\Omega }}=Tr(e^{\beta K})`$

$`\rho __G=Z_G^1e^{\beta K}=e^{\beta (\mathrm{\Omega }K)}`$

where $`\beta =1/T.`$ Then for any Schrödinger operator $`A_\alpha `$ the modified Heisenberg picture can be introduced

$`A_\alpha (\tau )=e^{K\tau }A_\alpha e^{K\tau }`$

and the temperature (or the Matsubara) Green’s Function (GF) is defined as

$`G_{\alpha \beta }^{\tau \tau ^{}}=T_\tau A_\alpha (\tau )A_\beta ^+(\tau ^{})=Tr\left[e^{\beta (K\mathrm{\Omega })}T_\tau e^{\tau K}A_\alpha e^{(\tau \tau ^{})K}A_\beta ^+e^{\tau ^{}K}\right]=`$

$$=Tr\left[\rho __GT_\tau e^{\tau K}A_\alpha e^{(\tau \tau ^{})K}A_\beta ^+e^{\tau ^{}K}\right]$$
$`(1)`$
Here, the brackets $``$ mean the thermodynamic average; $`T_\tau `$ is a $`\tau `$ ordering operator, which arranges operators with the earliest $`\tau `$ (the closest to $`\beta `$) to the right.

Let us consider the two-body Hamiltonian

$$H=\underset{12}{}t_{12}a_1^+a_2+\frac{1}{4}\underset{1234}{}\overline{v}_{1234}a_1^+a_2^+a_3a_4$$
$`(2)`$
where $`a,a^+`$ are fermion annihilation and creation operators; $`t_{12}`$ and $`\overline{v}_{1234}=v_{1234}v_{1243}`$ are the kinetic energy and the antisymmetrized matrix element of the two-body interaction. The Green’s function $`G_{\alpha \beta }^{\tau \tau ^{}}`$ for an arbitrary operator $`A_\alpha ^+`$ obeys the following equation of motion:

$`{\displaystyle \frac{}{\tau }}G_{\alpha \beta }^{\tau \tau ^{}}=\delta _{\tau \tau ^{}}[A_\alpha ,A_\beta ^+]T_\tau [A_\alpha ,K]^\tau A_\beta ^+(\tau ^{})=`$

$`=\delta _{\tau \tau ^{}}N_{\alpha \beta }+{\displaystyle \underset{\gamma }{}}{\displaystyle 𝑑\tau _1^{}_{\alpha \gamma }^{\tau \tau _1^{}}G_{\gamma \beta }^{\tau _1^{}\tau ^{}}}`$

In this expression it is possible to split the effective Hamiltonian $`_{\alpha \beta }^{\tau \tau ^{}}`$ into an instantaneous and a dynamic (frequency dependent) part

$`_{\alpha \beta }^{\tau \tau ^{}}`$ $`=`$ $`{\displaystyle \underset{\beta ^{}}{}}\left\{\delta _{\tau \tau ^{}}[[A_\alpha ,K],A_\beta ^{}^+]T_\tau [A_\alpha ,K]^\tau [K,A_\beta ^{}^+]^\tau ^{}_{irr}\right\}N_{\beta ^{}\beta }^1=`$
$`=`$ $`_{\alpha \beta }^{(0)}\delta _{\tau \tau ^{}}+_{\alpha \beta }^{(r)\tau \tau ^{}}`$
In the approximation of the instantaneous effective Hamiltonian i.e. neglecting $`_{\alpha \beta }^{(r)\tau \tau ^{}}`$, the Dyson equation for the two-body Matsubara GF $`G_{\alpha \beta }^{\tau \tau ^{}}`$ can be written as

$$\frac{}{\tau }G_{\alpha \beta }^{(0)\tau \tau ^{}}=\delta _{\tau \tau ^{}}N_{\alpha \beta }+\underset{\gamma }{}_{\alpha \gamma }^{(0)}G_{\gamma \beta }^{(0)\tau \tau ^{}}$$
$`(3)`$
In the treatment of two particle correlations let us specify the arbitrary operator $`A_\alpha `$ as $`A_{k_1k_2}=a_{k_1}a_{k_2}`$. In this case the Dyson equation (3) takes the following form in the frequency representation

$$i\omega _nG_{k_1k_2k_1^{}k_2^{}}^{SCRPA}=N_{k_1k_2k_1^{}k_2^{}}+\underset{p_1p_2}{}_{k_1k_2p_1p_2}^{(0)}G_{p_1p_2k_1^{}k_2^{}}^{SCRPA},$$
$`(4)`$
where, in supposing that the single particle density matrix is diagonal in the basis used (this is for example the case inhomogeneous matter):

$$N_{k_1k_2k_1^{}k_2^{}}=[a_{k_1}a_{k_2},a_{k_2^{}}^+a_{k_1^{}}^+]=\delta _{k_1k_2k_1^{}k_2^{}}\left(1n_{k_2^{}}n_{k_1^{}}\right)$$
$`(5)`$
and

$$_{k_1k_2k_1^{}k_2^{}}^{(0)}=\underset{p_1p_2}{}[[a_{k_1}a_{k_2},K],a_{p_2}^+a_{p_1}^+]N_{p_1p_2k_1^{}k_2^{}}^1$$
$`(6)`$
Here $`\delta _{k_1k_2k_1^{}k_2^{}}`$ is the antisymmetrized Kronecker symbol and $`n_k=a_k^+a_k`$ are the single particle occupation numbers which can be found from the single-particle Matsubara GF

$`G_{kk^{}}^\tau =T_\tau a_k(\tau )a_k^{}^+(0)`$

as

$$n_k=a_k^+a_k=\underset{\tau 0^{}}{lim}G_k^\tau $$
$`(7)`$
In general the single-particle Matsubara Green’s function $`G_k^{\tau \tau ^{}}`$ obeys the following Dyson equation:

$$\left(\frac{}{\tau }+\epsilon _k\right)G_{kk^{}}^\tau =\delta (\tau )+𝑑\tau _1M_k^{\tau \tau _1}G_k^{\tau _1}$$
$`(8)`$
or in the frequency representation

$$G_k^{\omega _n}=G_k^0+G_k^0M_kG_k^{\omega _n}$$
$`(9)`$
where

$$G_k^0=\frac{1}{i\omega _n\epsilon _k}$$
$`(10)`$
Here $`\epsilon _k`$ contains already the usual (instantaneous) mean field so that $`M_k`$ denotes only the dynamical part of the mass operator.

Now the problem is to find an approximation for the mass operator $`M_k`$ consistent with the SCRPA. A solution to this problem has been proposed in , which goes via the two body $`T`$ – matrix representation of the single particle mass operator , evaluating $`T`$ within SCRPA.

Let us add at this point a word of physical interpretation of the mean-field operator (6) . A quick look allows to realize that it contains no higher than two body correlation functions and therefore for their determination, with (4), one obtains a selfconsistency problem. Furthermore one can consider the nucleus as a gas of zero point pair fluctuations. These fluctuations create their own mean field, i.e. one pair fluctuation moves in the average potential created by all the other pair fluctuations. This average pair fluctuation field is graphically represented in Fig. 1. It gives rise, as usual, to a nonlinear problem. Of course, the single particle mean field introduced in (9) and further developed below is coupled to the selfconsistent pair potential. This is the deeper meaning of $`^{(0)}`$ of (6).

## III Application to the Picket Fence Model

The model consists of an equidistant multilevel pairing Hamiltonian with each level two fold degenerate, i.e. only spin up/down fermions of one kind can occupy one level. The corresponding Hamiltonian is given by

$$H=\underset{k=1}{\overset{\mathrm{\Omega }}{}}e_kN_kG\underset{i,k=1}{\overset{\mathrm{\Omega }}{}}P_i^+P_k$$
$`(11)`$
with

$$N_k=c_k^+c_k+c_{\overline{k}}^+c_{\overline{k}},$$

$$P_k^+=c_k^+c_{\overline{k}}^+$$
$`(12)`$
where $`\overline{k}`$ means the time reversed of $`k`$, single particle energies are $`e_k=k\epsilon \lambda ,`$ with level spacing $`\epsilon `$ chosen to be equal to $`1`$, and $`\mathrm{\Omega }`$ stands for the number of levels. The chemical potential $`\lambda `$ will be chosen such as to conserve the average number of particles $`N=\mathrm{\Omega }`$ of the system. The operators defined in (12) form an SU(2) algebra for each level $`j`$ and obey the following commutation relations

$$[P_j,P_k^+]=\delta _{jk}(1N_j),$$

$$[P_j,N_k]=2\delta _{jk}P_j,$$
$`(13)`$
### A SCRPA equations

To study the model at finite temperature we define in analogy to (1) the following set of two-body Matsubara GFs

$`G_{ji}^\tau =T_\tau \overline{P}_j(\tau )\overline{P}_i^+(0),`$

where

$$\overline{P}_j=\frac{P_j}{\sqrt{<\left|1N_j\right|>}}$$

Applying the instantaneous approximation for the mass operator we obtain the expressions for the two body SCRPA GF’s:

$$i\omega _nG_{ji}^{SCRPA}=\delta _{ji}+\underset{k}{}_{jk}^{(0)}G_{ki}^{SCRPA},$$
$`(14)`$
with

$$_{jk}^{(0)}=2\delta _{jk}\left(e_j+\frac{G}{<1N_j>}\underset{j^{}}{}<P_j^+P_j^{}>\right)G\frac{<(1N_j)(1N_k)>}{\sqrt{\left|<1N_j><1N_k>\right|}}$$
$`(15)`$
To find the correlation functions of the form $`<(1N_j)(1N_k)>`$ we will use the following approximation:

when $`jk`$ it is a simple factorization procedure, which has turned out to be accurate in the zero temperature limit:

$$<(1N_j)(1N_k)>=<1N_j><1N_k>$$
$`(16a);`$
but when $`j=k`$ we use the following exact relation

$$<(1N_j)^2>=<1N_j>+2<P_j^+P_j>$$
$`(16b)`$
which can easily be obtained taking into account that $`n_j^2=n_j`$ and $`n_jn_{\overline{j}}=P_j^+P_j`$ (here $`n_j=c_j^+c_j`$). It should be noted that in  the factorization (16a) was also used for the diagonal part (16b) and quite accurate results were obtained. We will show below that with (16b) one obtains still improved results. As shown in  it is possible to avoid above approximation. However, this is at the cost of a considerable numerical complication. We refrain from this here because it brings only very little improvement of results.

With this ansatz a particle-particle RPA-like equation is obtained

$$G_{ji}^{SCRPA}=\delta _{ji}\frac{1}{zC_j}\frac{G\sqrt{\left|D_jD_i\right|}}{\left(zC_j\right)\left(zC_i\right)}\times \left[1+G\underset{k}{}\frac{D_k}{zC_k}\right]^1,$$
$`(17)`$
where

$$z=i\omega _n,$$

$$D_i=<1N_i>$$

and

$`C_j=2\left(e_jGn_j+{\displaystyle \frac{G}{D_j}}{\displaystyle \underset{j^{}j}{}}<P_j^+P_j^{}>\right)`$

From this one easily can find the excitation spectrum of the model in equating the denominator of (17) to zero

$$1+G\underset{k}{}\frac{D_k}{zC_k}=0$$
$`(18)`$
Knowing the poles of the Green’s function (17), one can write down its spectral representation (we here give it as a function of imaginary time), with the corresponding residua:

$`G_{p_1p_2}^\tau =\theta \left(\tau \right)\sqrt{D_{p_1}D_{p_2}}\left[X_{p_1}^\mu X_{p_2}^\mu e^{E_\mu \tau }\left(1+n_B(E_\mu )\right)+Y_{p_1}^\mu Y_{p_2}^\mu e^{E_\mu \tau }n_B(E_\mu )\right]+`$

$$+\theta \left(\tau \right)\sqrt{D_{p_1}D_{p_2}}\left[X_{p_1}^\mu X_{p_2}^\mu e^{E_\mu \tau }n_B(E_\mu )+Y_{p_1}^\mu Y_{p_2}^\mu e^{E_\mu \tau }\left(1+n_B(E_\mu )\right)\right]$$
$`(19a)`$
$`G_{h_1h_2}^\tau =\theta \left(\tau \right)\sqrt{D_{h_1}D_{h_2}}\left[Y_{h_1}^\mu Y_{h_2}^\mu e^{E_\mu \tau }\left(1+n_B(E_\mu )\right)+X_{h_1}^\mu X_{h_2}^\mu e^{E_\mu \tau }n_B(E_\mu )\right]+`$

$$+\theta \left(\tau \right)\sqrt{D_{h_1}D_{h_2}}\left[Y_{h_1}^\mu Y_{h_2}^\mu e^{E_\mu \tau }n_B(E_\mu )+X_{h_1}^\mu X_{h_2}^\mu e^{E_\mu \tau }\left(1+n_B(E_\mu )\right)\right]$$
$`(19b)`$
$`G_{ph}^\tau =G_{hp}^\tau =\theta \left(\tau \right)\sqrt{\left|D_pD_h\right|}\left[X_p^\mu Y_h^\mu e^{E_\mu \tau }\left(1+n_B(E_\mu )\right)+Y_p^\mu X_h^\mu e^{E_\mu \tau }n_B(E_\mu )\right]+`$

$$+\theta \left(\tau \right)\sqrt{\left|D_pD_h\right|}\left[X_p^\mu Y_h^\mu e^{E_\mu \tau }n_B(E_\mu )+Y_p^\mu X_h^\mu e^{E_\mu \tau }\left(1+n_B(E_\mu )\right)\right]$$
$`(19c)`$
$`n_B\left(E_\mu \right)={\displaystyle \frac{1}{e^{\beta E_\mu }1}},`$

where the index $`p`$ refers to the states above Fermi level and the index $`h`$ to the ones below. The following amplitudes were introduced in these formulas:

$$X_p^\mu =\frac{\sqrt{GD_p}}{\left|C_p\right|E_\mu }F_\mu ,Y_h^\mu =\frac{\sqrt{GD_h}}{\left|C_h\right|+E_\mu }F_\mu $$
$`(20a)`$
$$X_h^\mu =\frac{\sqrt{GD_h}}{\left|C_h\right|E_\mu }F_\mu ,Y_p^\mu =\frac{\sqrt{GD_p}}{\left|C_p\right|+E_\mu }F_\mu $$
$`(20b)`$
with

$$F_\mu ^2=\frac{}{z}\left[1+G\underset{k}{}\frac{D_k}{zC_k}\right]_{z=E_\mu }$$
$`(20c)`$
These amplitudes obey the usual normalization conditions

$$\underset{p}{}X_p^\mu X_p^\mu ^{}+\underset{h}{}Y_h^\mu Y_h^\mu ^{}=\delta _{\mu \mu ^{}}$$

$$\underset{h}{}X_h^\mu X_h^\mu ^{}+\underset{p}{}Y_p^\mu Y_p^\mu ^{}=\delta _{\mu \mu ^{}}$$
$`(21)`$
Two-body correlation functions can be obtained from the Green’s function (19) as follows

$`<P_{p_1}^+P_{p_2}>=G_{p_1p_2}^{\tau 0^{}}=\sqrt{D_{p_1}D_{p_2}}\left({\displaystyle \underset{\mu }{}}X_{p_1}^\mu X_{p_2}^\mu n_B(E_\mu )+{\displaystyle \underset{\mu }{}}Y_{p_1}^\mu Y_{p_2}^\mu \left(1+n_B(E_\mu )\right)\right)`$

$$<P_{h_1}^+P_{h_2}>=G_{h_1h_2}^{\tau 0^{}}=\sqrt{D_{h_1}D_{h_2}}\left(\underset{\mu }{}Y_{h_1}^\mu Y_{h_2}^\mu n_B(E_\mu )+\underset{\mu }{}X_{h_1}^\mu X_{h_2}^\mu \left(1+n_B(E_\mu )\right)\right)$$
$`(22)`$
$`<P_p^+P_h>=G_{ph}^{\tau 0^{}}=\sqrt{\left|D_pD_h\right|}\left({\displaystyle \underset{\mu }{}}X_p^\mu Y_h^\mu n_B(E_\mu )+{\displaystyle \underset{\mu }{}}Y_p^\mu X_h^\mu \left(1+n_B(E_\mu )\right)\right)`$

### B Occupation numbers in the SCRPA

In order to close the set of the SCRPA equations, it is necessary to find the so far unknown occupation numbers $`n_k=<c_k^+c_k>`$ . For this, we should find the single particle Green’s function $`G_k^\tau `$ consistent with the SCRPA scheme. As discussed in sect. II, the single particle mass operator $`M_k`$ has in general the exact representation in terms of the two body $`T`$-matrix  and then an appropriate approximation for the $`G_k^\tau `$ can be obtained. It consists in using the mass operator $`\stackrel{~}{M_k}`$ calculated through the $`T`$-matrix found in the framework of SCRPA. As the relation between the $`T`$-matrix and the sum of the all irreducible Feynman graphs in the pp-chanel is also known then the following expression for the single particle mass operator can be obtained :

$$\stackrel{~}{M_k}=G\underset{k_1k_2}{}G_{\overline{k}}^{0(\tau _1^{}\tau _1)}G_{k_1k_2}^{(\tau _1^{}\tau _1)}\stackrel{~}{}_{k_2k}^{(0)}$$
$`(23)`$
$`\stackrel{~}{}_{k_2k}^{(0)}`$ is expressed through the effective Hamiltonian (15) without the disconnected part

$$\stackrel{~}{}_{kk}^{(0)}=_{kk}^{(0)}2\delta _{kk}\epsilon _k$$
$`(24)`$
where $`\epsilon _k`$ is defined below in (30).

In addition to this transparent scheme there also exists an additional consistency requirement . It follows from the possibility to calculate the average value of the Hamiltonian $`<H>`$ in two ways. On the one hand one has the following relation between the single particle Green’s function (9) and $`<H>`$ :

$$<H>=\frac{1}{2}\underset{\tau ^{}\tau 0^+}{lim}\underset{k}{}\left[\frac{}{\tau }e_k\right]\left(G_k^{(\tau \tau ^{})}+G_{\overline{k}}^{(\tau \tau ^{})}\right)$$
$`(25)`$
On the other hand there exists the straightforward calculation of $`<H>`$ through the two body Green’s functions (22)

$$<H>=\underset{k}{}e_k<N_k>G\underset{k_1k_2}{}<P_{k_1}^+P_{k_2}>$$

$$=\frac{N^2}{4}+2\underset{p}{}e_p<N_p>2\underset{p}{}\underset{\mu }{}\sqrt{GD_p}F_\mu \left(n_B(E_\mu )X_p^\mu +(1+n_B(E_\mu ))Y_p^\mu \right)$$
$`(26)`$
where we used the particle-hole symmetry of the model  and reduced sums over p and h only to the one over the particle states.

The additional consistency condition lies in the requirement that both expressions (25) and (26) should give exactly the same results. This only is satisfied if the single particle GF is expanded to first order in the renormalized single particle mass operator (23) (one may verify that this is in analogy to the standard RPA scheme, i.e. the standard RPA average energy is obtained via (25) using a single particle GF with only perturbative renormalization from the RPA-modes):

$$G_p=G_p^0+G_p^0M_p^{SCRPA}G_p^0$$
$`(27)`$
with

$$M_p^{SCRPA}=\frac{\sqrt{D_p}}{D_p^0}\stackrel{~}{M_p},$$
$`(28)`$
$$D_p^0=1f_pf_{\overline{p}}$$

$$f_p=\frac{1}{1+e^{\epsilon _p\beta }},$$
$`(29)`$
and

$$\epsilon _p=e_pGf_p\frac{f_p}{n_p}\frac{D_p}{D_p^0}.$$
$`(30)`$
Finally we get the following expression for the SCRPA single particle mass operator $`M_p^{SCRPA}`$:

$$M_p^{SCRPA}=\frac{\sqrt{GD_p}}{D_p^0}\underset{\mu }{}F_\mu \left[\frac{S_{\mu p}^{(1)}(f_p+n_B(E_\mu ))}{i\omega _n+\epsilon _pE_\mu }+\frac{S_{\mu p}^{(2)}(1f_p+n_B(E_\mu ))}{i\omega _n+\epsilon _p+E_\mu }\right]$$
$`(31)`$
where

$$S_{\mu p}^{(1)}=\underset{p^{}}{}X_p^{}^\mu \stackrel{~}{}_{pp^{}}+\underset{h^{}}{}Y_h^{}^\mu \stackrel{~}{}_{ph^{}},$$

$$S_{\mu p}^{(2)}=\underset{p^{}}{}Y_p^{}^\mu \stackrel{~}{}_{pp^{}}+\underset{h^{}}{}X_h^{}^\mu \stackrel{~}{}_{ph^{}}$$

The corresponding single particle occupation numbers $`n_p`$, found from (7), is the following

$$n_p=<c_p^+c_p>=f_p+\frac{\sqrt{GD_p}}{D_p^0}\underset{\mu }{}F_\mu \left[X_p^\mu \frac{n_B(E_\mu )D_p^0f_p^2}{\left(2\epsilon _pE_\mu \right)}+Y_p^\mu \frac{(n_B(E_\mu )+1)D_p^0+f_p^2}{\left(2\epsilon _p+E_\mu \right)}\right]$$

$$f_p(1f_p)\beta \frac{\sqrt{GD_p}}{D_p^0}\underset{\mu }{}F_\mu [X_p^\mu (f_p+n_B(E_\mu )+Y_p^\mu (1f_p+n_B(E_\mu )]$$
$`(32)`$
Let us demonstrate now that using the single particle Green’s function (27) with the mass operator (31) and occupation numbers (32) in the calculation of the average energy (25) indeed leads to the equation (26). At first one finds the derivative of the single particle GF (27)

$$\frac{G_p^{(\tau \tau ^{})0^+}}{\tau }=\epsilon _pf_p+\frac{\sqrt{GD_p}}{D_p^0}\underset{\mu }{}F_\mu \{(1f_p)[S_{\mu p}^{(1)}\frac{(\epsilon _pE_\mu )n_B(E_\mu )}{(2\epsilon _pE_\mu )^2}+S_{\mu p}^{(2)}\frac{(\epsilon _p+E_\mu )(1+n_B(E_\mu ))}{(2\epsilon _p+E_\mu )^2}]$$

$$+\epsilon _pf_p\left[S_{\mu p}^{(1)}\frac{f_p+n_B(E_\mu )}{(2\epsilon _pE_\mu )^2}+S_{\mu p}^{(2)}\frac{1f_p+n_B(E_\mu )}{(2\epsilon _p+E_\mu )^2}\right]$$

$$+\beta \epsilon _pf_p(1f_p)\left[S_{\mu p}^{(1)}\frac{f_p+n_B(E_\mu )}{2\epsilon _pE_\mu }+S_{\mu p}^{(2)}\frac{1f_p+n_B(E_\mu )}{2\epsilon _p+E_\mu }\right]$$

$$f_p[S_{\mu p}^{(1)}\frac{f_p+n_B(E_\mu )}{2\epsilon _pE_\mu }+S_{\mu p}^{(2)}\frac{1f_p+n_B(E_\mu )}{2\epsilon _p+E_\mu }]\}=$$

$$=\epsilon _pn_p+\underset{\mu }{}\frac{\sqrt{GD_p}F_\mu }{D_p^0}\left[S_{\mu p}^{(1)}\frac{n_B(E_\mu )D_p^0f_p^2}{2\epsilon _pE_\mu }+S_{\mu p}^{(2)}\frac{n_B(E_\mu )D_p^0+f_p^2}{2\epsilon _p+E_\mu }\right]$$

Inserting this expression in (25) we obtain

$$<H>=\frac{1}{2}\underset{\tau ^{}\tau 0^+}{lim}\underset{k}{}\left[\frac{}{\tau }e_k\right]\left(G_k^{(\tau \tau ^{})}+G_{\overline{k}}^{(\tau \tau ^{})}\right)=$$

$$=\frac{N^2}{4}+2\underset{p}{}(e_p+\epsilon _p)n_p2\underset{p\mu }{}\sqrt{GD_p}F_\mu \left(n_B(E_\mu )X_p^\mu +(1+n_B(E_\mu ))Y_p^\mu \right)+G\underset{p}{}\frac{D_p}{D_p^0}f_p^2=$$

$$=\frac{N^2}{4}+4\underset{p}{}e_pn_p2\underset{p}{}\underset{\mu }{}\sqrt{GD_p}F_\mu \left(n_B(E_\mu )X_p^\mu +(1+n_B(E_\mu ))Y_p^\mu \right)$$

This is exactly equal to (26).

The system of the SCRPA equations is fully closed now. Together with (17), (19) and (22) this represents a self-consistent problem for pair fluctuations.

We want to indicate at this point that the above way to determine the single particle occupancies is not the only possibility. In the Appendix we will give another variant which, however, yields results close to the ones with the method of this section. The non uniqueness of the occupation numbers reflects the fact that with the truncated ansatz (3), at zero temperature, no corresponding ground state wave function can be found, as explained in . For a wave function to exist, the ansatz (3) must be extended. It can, however, be shown that the correction terms are small .

### C Exact statistical treatment of the PFM

For an exact statistical treatment of the Picket Fence Model we have to find all exact eigenvalues and eigenstates of the Hamiltonian (11). Since singly-occupied levels do not participate in the pair scattering, eigenstates can be classified according to the number of unpaired particles $`S`$ (seniority). There are $`C_{\mathrm{\Omega }S}^\mathrm{\Omega }=\frac{\mathrm{\Omega }!}{S!(\mathrm{\Omega }S)!}`$ different multiplets of this type, each of dimension $`C_{\frac{NS}{2}}^{\mathrm{\Omega }S}`$ and degeneracy $`2^S`$. If we define the following set of basis states for each multiplet:

$$|\{s_i,N_i\}>$$
$`(33)`$
where $`s_i=N_i=1`$ for singly-occupied levels and $`s_i=0,N_i=0`$ or $`2`$ for remaining levels ($`_is_i=S`$), the Hamiltonian matrix will have the following diagonal and off-diagonal elements

$$\{s_i,N_i\}\left|H\right|\{s_i,N_i\}=\underset{kS}{}\left(e_k\lambda \right)+\underset{k\mathrm{\Omega }S}{}\left[\left(e_k\lambda \right)\frac{G}{4}\left(4N_k\right)\right]N_k$$
$`(34)`$
$$s_{j_1}N_{j_1},\mathrm{},s_{j_k}2,\mathrm{},s_{j_i}0,\mathrm{},s_{j_n}N_{j_n}\left|H\right|s_{j_1}N_{j_1},\mathrm{},s_{j_k}0,\mathrm{},s_{j_i}2,\mathrm{},s_{j_n}N_{j_n}=G$$
$`(35)`$
The exact eigenvalues and eigenstates can be calculated by diagonalization of this matrix in each multiplet. The exact grand canonical average $`<A>`$ of any operator can then be obtained with the help of the grand partition function $`Z_G=Tre^{\beta \widehat{H}}`$ and the statistical operator $`\widehat{\rho _G}=Z_G^1e^{\beta \widehat{H}}`$ as

$$<A>=Tr[\widehat{A}\widehat{\rho _G}]$$
$`(36)`$
## IV Results and Discussion

In order to check the accuracy of our theory and of the different approximations schemes we first calculate the average energy of the system $`<H>`$ as a function of the particle number $`N`$ and temperature $`T`$. The results of the different calculations are presented in Figures 2 – 4. Calculations were made for a value of the pairing constant $`G`$ which is smaller but close to the critical value $`G_{cr}`$ at $`T=0`$. The phase transition from the normal to the superfluid phase occurs in the system when $`GG_{cr}`$ . We compare the SCRPA results with the exact ones for the Grand Canonical Ensemble (GCE) as well as results of the standard Thermal RPA (TRPA) and Thermal MFA (TMFA). One can see that when the number of levels $`\mathrm{\Omega }`$ (and number of particles $`N`$) increases the description of the intrinsic energy becomes better and at $`\mathrm{\Omega }=N=10`$ the TSCRPA results practically coincide with the exact ones. Especially the last case will be considered below more carefully.

Let us now come to the discussion of the behavior of the excitation energies. The dependence of the excitation energies of the addition mode (see ) as a function of $`G`$ is shown in Fig. 5 at zero temperature. The SCRPA (solid lines) is compared with the standard RPA calculations (dashed lines) and the exact ones (open circles). Increasing the interaction constant, the lowest energy in the RPA goes to zero and at $`G_{cr}0.33`$ the collapse takes place which is connected with the transition from the normal to the superfluid phase. At finite temperature (see Fig. 6 where the dependence of the lowest addition mode is presented as a function of $`G`$ at $`T=0.5`$) this collapse occurs at a higher value of the interaction constant ($`G_{cr}0.43`$) what is due to the reduced intensity of the residual interaction because of the thermal factors. This collapse is absent in the exact calculations at zero temperature and also in the SCRPA calculations at zero and finite temperatures. It is also remarkable that the SCRPA yields a rise of all excitation energies with increasing $`G`$ in contrast with RPA and in very good agreement with the exact results. This comes from the fact that in the PFM with the Kramer’s degeneracy of levels the Pauli repulsion is extremely strong overruling the original attractive interaction. In this model, therefore, standard RPA gives qualitatively wrong result.

We next consider the behavior of the system near the phase transition point. To make distinctions between different results more apparent we not only show the full intrinsic energy $`<H>`$ but also the correlation energy $`E_{corr}`$ which is defined as

$$E_{corr}=<H><H>_0$$
$`(37)`$
where $`<H>_0`$ is the average energy calculated in Mean Field Approximation. In Figures 7 and 8, the average energy $`<H>`$ and correlation energy $`E_{corr}`$ as a function of $`T`$ are displayed for the interaction constant $`G=0.4`$ (at $`T=0`$ this value of $`G`$ is larger than $`G_{cr}0.33`$). With increasing $`T`$ the mean field rearrangement occurs and the system goes from the superfluid phase to the normal one at $`T_c0.38`$. Note, that within the TRPA the lowest excitation energy $`\omega _1`$ vanishes when $`TT_c`$, whereas within the TSCRPA $`\omega _1`$ stays finite. For both the correlation energy and the intrinsic energy the TSCRPA gives more precise results as compared to the other approximations. It is remarkable that the TSCRPA results are accurate down to practical zero temperature, in spite of the fact that within standard BCS theory one enters the superfluid regime. A quasiparticle formulation of SCRPA  will only be necessary for stronger $`G`$ values driving the system more deeply into the symmetry broken phase.

To analyze the region near the phase transition point in more detail, the heat capacity is calculated as a partial derivative of the intrinsic energy with respect to $`T`$

$$C_v=\frac{<H>}{T}$$
$`(38)`$
The results are shown in Fig. 9. The TRPA and TMFA give discontinuities of $`C_v`$ at $`T_c`$ (we recall again that our results are obtained using a normal fluid approach and not transforming to quasiparticles). The heat capacity calculated in the TSCRPA has some kink near $`T_c`$ but has no discontinuities and is quite similar to the exact result through out the whole range of temperature.

Nevertheless, the TSCRPA and also the exact solution feel the phase transition to the superfluid phase. It can be seen in Fig. 8 where both the TSCRPA and the exact correlation energies show a depression near $`T=0.38`$. This originates from strong pair fluctuations leading to the BCS phase transition in TMFA with the critical temperature $`T_{cr}^{BCS}=0.38`$ for $`G=0.4`$. However, one notices (see also ) that the sharp phase transition of mean field is in reality completely smeared out and only a faint, though clearly visible, signal survives.

To investigate the formation of such fluctuations in more detail, it is useful to consider the spectral function

$$A(k,\omega )=\frac{1}{i}\left[G_k^{\omega _n}|_{i\omega _n\omega i\eta }G_k^{\omega _n}|_{i\omega _n\omega +i\eta }\right]$$
$`(39)`$
which includes two-body correlations through the self-consistent treatment of the mass operator $`M_k^{SCRPA}(\omega _n)`$. Based on the spectral function the density of states can be calculated as

$$N(\omega )=\underset{k}{}A(k,\omega )$$
$`(40),`$
The results of the calculations of $`N(\omega )`$ with $`G_k`$ of (27) in (39) are shown in Fig. 10 at different values of $`T`$ and for $`G`$ = 0.4. It is clearly seen that the distance between the two quasiparticle peaks around the Fermi energy $`\epsilon _F`$ ($`\omega =0`$) increases with decreasing temperature. This process sets in even above the BCS transition temperature $`T_c^{BCS}=0.38`$. This rarefaction of the level density around $`\epsilon _F`$ above $`T_c`$ is not avoid of similarity with the situation in high $`T_c`$ – superconductors where a so-called ’pseudo gap’ in the level density appears already above $`T_c`$ . This ’pseudo gap’ also is often attributed to a decrease in the level density around $`\epsilon _F`$ due to pair fluctuations . Apparently it is a quite generic feature that pair correlations diminish the density of levels around $`\epsilon _F`$ whereas particle-hole correlations give rise to an increase.

In order to make the temperature dependence of the gap more transparent let us introduce an effective (or canonical) gap which recently was proposed in Eq. (22) of Ref.

$$\mathrm{\Delta }=G\sqrt{\underset{ik}{}\left(P_i^+P_kc_i^+c_kc_{\overline{i}}^+c_{\overline{k}}\right)}$$
$`(41)`$
In the BCS approximation, the effective gap $`\mathrm{\Delta }`$ coincides with the usual grand canonical BCS gap

$$\mathrm{\Delta }_0=G\underset{k}{}P_k_{BCS}$$

The dependencies of the effective gap on the interaction constant $`G`$ at zero temperature and on temperature $`T`$ at $`G=0.4`$ are shown in Figures 11 and 12. The SCRPA results give a very good description of the gap at zero and non zero temperatures. It is clearly seen that the SCRPA and exact calculations do not display the phase transition at the point where the BSC gap vanishes. Notice that in Fig. 11 the SCRPA result deteriorates for values of $`G`$ well beyond the critical value. In this regime a quasiparticle generalization of the SCRPA is necessary  because one enters deeply in the superfluid region.

## V Comparison with other works

The PFM has recently widely been used for the study of quantum pair fluctuations at finite temperature both in the context of nuclear physics  and of ultrasmall metallic grains . In both fields the SPA approach where one additionally takes into account number parity projection and quantum (RPA) fluctuation around mean field was employed.

To compare our results with the above mentioned formalism let us introduce the relevant energy scales. These are the average level spacing $`\delta `$ and the BCS energy $`\mathrm{\Delta }=\mathrm{\Omega }\delta /[2\mathrm{sinh}(\delta /G)]`$ . The properties of the system described by the pairing Hamiltonian can be calculated as universal functions of the single scaling parameter $`\delta /\mathrm{\Delta }`$. As long as the grain is not too small $`\delta <<\mathrm{\Delta }`$, the fluctuation region around $`T_c`$ is narrow, and the mean field (BCS) description of superconductivity is appropriate. It is the case in ref. , where the PFM is investigated with characteristic values of $`\delta /\mathrm{\Delta }0.2`$. The result of that work shows that the domain where the parity projection is important lies in a small region near $`T_c`$. At temperatures higher than $`T_c`$ the role of the fluctuation is decreased and usual SPA becomes rather good in reproducing the exact canonical results.

When the size of the system is decreased, fluctuations start to smear the normal – superconducting transition. The finite level spacing suppresses the BCS gap and when $`\delta `$ becomes of the order of $`\mathrm{\Delta }`$, the fluctuation region becomes of order $`T_c`$ and the BCS description of superconductivity breaks down even at zero temperature. Returning to our calculations, we can see that SCRPA yields the best results for $`\delta /\mathrm{\Delta }1`$ (see Fig. 11, 12). This region corresponds to ultrasmall grains where strong pairing fluctuations are dominant. In this sense our results can be compared (at least qualitatively) with the results of . To make this comparison more accurate we calculated average energy and specific heat for a system with 50 fermions and interaction constant $`G=0.127`$ and $`G=0.255`$, what corresponds to $`\delta /\mathrm{\Delta }=50`$ and $`\delta /\mathrm{\Delta }=1`$ respectively. In general our results correlate well with . From Fig. 13, where $`C_v`$ is displayed as a function of $`T/\delta `$, we can see that when $`\delta /\mathrm{\Delta }\mathrm{}`$ the specific heat approaches a linear dependence with temperature while when $`\delta \mathrm{\Delta }`$ a bump structure arises at low temperature which is a sign of the presence of strong pairing fluctuations. As it has been seen from our previous discussion, the SCRPA gives a better description with increasing particles number. And while we did not perform exact GCE calculations for the system with 50 particles, we can expect from our studies above and in  that our result should be very close to the exact one.

In conclusion of this section we can say that the results of TSCRPA are at least as good as the ones of SPA with extensions. However, contrary to SPA, TSCRPA has no problem at low temperature and excitation energies and correlation functions can be calculated directly as a function of temperature.

## VI Conclusion

In this work we generalized our recent work  on the multilevel pairing model (PFM) within the SCRPA approach to finite temperature (TSCRPA). The PFM has been recognised to account in many respects for the physics of (superconducting) metallic nano-grains. In our context SCRPA can be viewed as a self-consistent mean field theory for pair fluctuatuations. The results at $`T=0`$ in  are in very close agreement with the exact ones obtained from the Richardson procedure . It is therefore an important issue to also exploit SCRPA at finite $`T`$ and to assess its accuracy with respect to exact results in this case. Our comparison is mostly done for the case of ten levels with $`10`$ particles where it is still of some ease to establish the exact partition function. We, however, also considered with TSCRPA the case of 50 particles in $`50`$ levels assuming that the results are of equal quality or even better as the ones obtained for $`10`$ particles. We base our studies on the Matsubara $`1`$ and $`2`$ particle Green’s functions which allows us to calculate correlation and excitation energies, specific heat, level densities, etc. It can be considered as a general advantage of our approach that all these quantities are directly accessible in the whole range of temperatures and coupling constants. For the latter this holds in this work only true for interaction values not driving the system deeply into the superfluid regime, since in this work we only have been working with normal particles and not with quasiparticles. For $`G>>G_{cr}`$ we have to employ the Self Consistent Quasiparticle RPA (SCQRPA). It has recently been demonstrated in the two level pairing model that also SCQRPA gives very promising results . The quality of our results for the above mentioned quantities turns out to be excellent and it does not fail in any qualitative nor quantitative aspect. Most of the time the agreement with the exact solution is within the couple of percent level. One particularly interesting feature of our investigation is the fact that we achieved to calculate the single particle Green’s function consistently within TSCRPA. This enabled us to give, for the first time, for the PFM the evolution of the single particle level density with temperature. The construction of the exact solution for this quantity is very cumbersome and we refrained from working this out here. However, backed with the positive experience for all other above mentioned quantities we believe that also the level density is reasonably accurate. The interesting aspect of our calculation is that with decreasing temperature the density of single particle states around the Fermi level decreases even above the critical temperature as defined by BCS – theory. It is suggestive to see this feature in analogy to the appearance of a so called pseudogap in high $`T_c`$ superconductors where also a depression in the level density is observed approaching $`T_c`$ from above . It would be interesting to attempt an experimental verification with metallic nanograins of our prediction that indeed the density of levels rarefies with decreasing temperature already in the non superfluid regime.

We also gave a short comparison of TSCRPA with results at finite $`T`$ obtained with other approaches like the static path approximation (SPA) to the partition function. Though the results seem generally comparable, we think that TSCRPA is more versatile, giving direct access to correlation functions, level densities, excitation energies, etc. in the whole temperature range, quantities which are otherwise difficult to obtain.

In the present work we restricted ourselves to values of the coupling which are below or slightly above the critical value. In the future we shall elaborate on the SCRPA for quasiparticles at finite $`T`$ (TSCQRPA) which will allow us to consider the system deeply in the superfluid phase and to study the transition from one phase to the other in more detail.

ACKNOWLEDGMENTS

We acknowledge useful discussions with P. Bozek, F. Hekking and J. Hirsch.

APPENDIX: OCCUPATION NUMBER – VARIANT

Let us give another possible way to find the occupation numbers at finite temperature. Before to complete this task we firstly derive the SCRPA equations with the Equation of Motion method and find some useful relations between phonon amplitudes at zero temperature. To attain this let us introduce the pair addition and removal operators (phonons)  as

$$Q_\mu ^+=\underset{p}{}X_p^\mu \overline{P}_p^+\underset{h}{}Y_h^\mu \overline{P}_h^+$$
$`(A1)`$
$$Q_\lambda ^+=\underset{h}{}X_h^\lambda \overline{P}_h\underset{p}{}Y_p^\lambda \overline{P}_p$$
$`(A2)`$
and apply the variational procedure where all expectation values are found with respect to the vacuum of phonons (A1,A2).

$$<[\delta Q_\nu ,[,Q_\nu ^+]]>=\omega _\nu <[\delta Q_\nu ,Q_\nu ^+]>$$
$`(A3)`$
If we use the factorisation procedure (16) the following system of equations for the phonon amplitudes is obtained

$$\underset{p^{}}{}X_p^{}^\mu A_{pp^{}}\underset{h^{}}{}Y_h^{}^\mu B_{ph^{}}=E_\mu X_p^\mu $$

$$\underset{p^{}}{}X_p^{}^\mu B_{hp^{}}\underset{h^{}}{}Y_h^{}^\mu C_{hh^{}}=E_\mu Y_h^\mu ,$$
$`(A4)`$
where

$$A_{pp^{}}=2\delta _{pp^{}}[e_p]+2G\delta _{pp^{}}\left[\underset{lp}{}\frac{<P_l^+P_p>}{D_p}n_p\right]G\sqrt{D_pD_p^{}}=2\delta _{pp^{}}\epsilon _p+\stackrel{~}{A}_{pp^{}}$$

$$B_{ph}=G\sqrt{D_pD_h}$$

$$C_{hh^{}}=2\delta _{hh^{}}[e_hG]+2G\delta _{hh^{}}\left[\underset{lh}{}\frac{<P_l^+P_h>}{D_h}+1n_h\right]G\sqrt{D_hD_h^{}}=\left(2\delta _{hh^{}}\epsilon _h+\stackrel{~}{C}_{hh^{}}\right)$$

$$B_{hp}=G\sqrt{D_hD_p}$$
$`(A5)`$
where the single-particle energies are introduced as

$$\epsilon _p=e_p,\epsilon _h=e_h+G$$
$`(A6)`$
Thereafter the expressions for the phonon amplitudes are obtained as

$$X_p^\mu (2\epsilon _pE_\mu )=\underset{p^{}}{}X_p^{}^\mu \stackrel{~}{A}_{pp^{}}+\underset{h^{}}{}Y_h^{}^\mu B_{ph^{}}X_p^\mu =\frac{_p^{}X_p^{}^\mu \stackrel{~}{A}_{pp^{}}+_h^{}Y_h^{}^\mu B_{ph^{}}}{2\epsilon _pE_\mu }$$
$`(A7a)`$
$$Y_h^\mu (2\epsilon _h+E_\mu )=\underset{p^{}}{}X_p^{}^\mu B_{hp^{}}+\underset{h^{}}{}Y_h^{}^\mu \stackrel{~}{C}_{hh^{}}Y_h^\mu =\frac{_p^{}X_p^{}^\mu B_{hp^{}}+_h^{}Y_h^{}^\mu \stackrel{~}{C}_{hh^{}}}{2\epsilon _h+E_\mu }$$
$`(A7b)`$
Due to the particle-hole symmetry, the removal mode satisfy exactly the same equations. It means that

$$X_p^\mu =X_{h=Np+1}^{\lambda =\mu },Y_h^\mu =Y_{p=Nh+1}^{\lambda =\mu }$$
$`(A8)`$
Below we will use the notation $`\mu `$ for both modes.

Returning to the single particle occupation numbers $`n_k`$ let us remind that in the picket fence model the following exact relation between one body operator $`N_k`$ and two body operator $`P_k^+P_k`$ is verified for any non singly-occupied level $`k`$:

$$N_k=2P_k^+P_k$$
$`(A9)`$
Taking the average of both parts of this relation with respect to the SCRPA vacuum state we obtain

$$<N_p>=2D_p\underset{\mu }{}(Y_p^\mu )^2$$
$`(A10)`$
On the other hand, using relations (A7) we can rewrite this expression as

$$<N_p>=2D_p\underset{\mu }{}\frac{1}{\left(E_\mu +2ϵ_p\right)^2}\left[\underset{h}{}X_h^\mu B_{ph}\underset{p^{}}{}Y_p^{}^\mu \stackrel{~}{C}_{pp^{}}\right]\left[\underset{h}{}X_h^\mu B_{hp}\underset{p^{}}{}Y_p^{}^\mu \stackrel{~}{C}_{p^{}p}\right]$$
$`(A11)`$
It is easy to show that the fraction in this expression can be expressed through the GFs as

$$\frac{1}{\left(E_\mu +2ϵ_p\right)^2}=\underset{t^{}t0^+}{lim}𝑑t_1^{}𝑑t_1G_p^{0(tt_1)}\theta (t_1^{}t_1)e^{i(E_\mu +\epsilon _p)(t_1^{}t_1)}G_p^{0(t_1^{}t^{})}$$
$`(A12)`$
where

$$G_p^{0(tt_1)}=i<T_tc_p(t)c_p^+(t_1)>=i\theta (tt_1)e^{i\epsilon _p(tt_1)}$$
$`(A13)`$
If we now take into account the spectral representation of the pair operator $`P_p`$

$$P_p(t)=D_p\left[\underset{\mu }{}e^{iE_\mu t}Q_\mu X_p^\mu +\underset{\mu }{}e^{iE_\mu t}Q_\mu ^+Y_p^\mu \right]$$
$`(A14)`$
and define the antichronological (time reversed) single particle GF as

$$G_{\overline{p}}^{0(tt_1)}=i<T_tc_{\overline{p}}^+(t)c_{\overline{p}}(t_1)>=i\theta (t_1t)e^{i\epsilon _p(tt_1)}$$
$`(A15)`$
we can pass to the expression for the sp GF which gives occupation numbers which are consistent within the frame of the SCRPA with exact relation (A9)

$$G_p^{\left(tt^{}\right)}=G_p^{0(tt^{})}+𝑑t_1𝑑t_1^{}G_p^{0(tt_1)}M_p^{SCRPA}G_p^{0(t_1^{}t^{})}$$
$`(A16)`$
where $`M_p^{SCRPA}`$ is a single particle mass operator

$$M_p^{SCRPA}=\underset{k_1k_2}{}G_{\overline{p}}^{0(t_1^{}t_1)}\overline{}_{pk_1}^{(0)}G_{k_1k_2}^{(t_1^{}t_1)}\overline{}_{pk_2}^{(0)}$$
$`(A17)`$
and $`\overline{}_{pk}^{(0)}`$ is the renormalised effective Hamiltonian which has the following form in terms of RPA matrixes $`A_{pp^{}}`$ and $`B_{ph}`$

$$\overline{}_{pp^{}}^{(0)}=\sqrt{D_p}\stackrel{~}{A}_{pp^{}}\frac{1}{\sqrt{D_p^{}}}$$
$`(A18a)`$
$$\overline{}_{ph}^{(0)}=\sqrt{D_p}B_{ph}\frac{1}{\sqrt{D_h}}$$
$`(A18b)`$
From (A16) and (A17) we then can define the occupation numbers as usual. In general, the occupation numbers from (32) and (A16), (A17) have slightly different values. This is due to the fact that the ansatz (A1), (A2) for the RPA operators is too restricted for a groundstate fulfilling $`Q|0>=0`$ to exist. It can, however, be shown  that the necessary corrections to (A1), (A2) are small. Anyhow, the differences in occupation numbers obtained from the two methods are a measure of the importance of the terms neglected in (A1), (A2). The difference of the mass operators (23) and (A17) is that in the latter both vertices are dressed, whereas in (23) one vertex remains at the unrenormalized value (G).

To find occupation numbers $`<n_p>`$ at finite temperature we adopt the expressions obtained for single particle GF at zero temperature (eq.(A16)-(A18)). To do it, it is necessary to change all zero temperature GFs to the Matsubara ones and use for the vertices the renormalised effective Hamiltonian (15)

$$\overline{}_{pk}^{(0)}=\sqrt{D_p}\stackrel{~}{}_{pk}^{(0)}\frac{1}{\sqrt{\left|D_k\right|}}$$
$`(A19)`$
where

$$\stackrel{~}{}_{pk}^{(0)}=_{pk}^{(0)}2\delta _{pk}\epsilon _k$$

and

$$\epsilon _k=e_kGf_k,$$

$$f_k=\frac{1}{1+e^{\epsilon _k\beta }}$$

Then we get for the single particle Matsubara GF

$$G_p^{\left(\tau \tau ^{}\right)}=G_p^{0(\tau \tau ^{})}+𝑑\tau _1𝑑\tau _1^{}G_p^{0(\tau \tau _1)}M_p^{SCRPA}G_p^{0(\tau _1^{}\tau ^{})}$$

$$=G_p^{0(\tau \tau ^{})}+𝑑\tau _1𝑑\tau _1^{}G_p^{0(\tau \tau _1)}\left[\underset{k_1k_2}{}G_{\overline{p}}^{0(\tau _1^{}\tau _1)}\overline{}_{pk_1}^{(0)}G_{k_1k_2}^{(\tau _1^{}\tau _1)}\overline{}_{pk_2}^{(0)}\right]G_p^{0(\tau _1^{}\tau ^{})}$$
$`(A20)`$
where

$$G_p^{0(\tau \tau ^{})}=e^{\epsilon _p(\tau \tau ^{})}\left[\theta \left(\tau \tau ^{}\right)(1f_p)\theta \left(\tau ^{}\tau \right)f_p\right],$$

Taking this integral in the limit $`\tau \tau ^{}0^{}`$ we obtain the occupation numbers $`n_p`$ in the SCRPA at finite temperature

$`n_p=<c_p^+c_p>=f_p+D_p{\displaystyle \underset{\mu }{}}\left(Y_p^\mu \right)^2\left[n_B\left(E_\mu \right)\left(12f_p\right)+\left(1f_p\right)^2\right]`$

$`+D_p{\displaystyle \underset{\mu }{}}\left(X_p^\mu \right)^2\left[n_B\left(E_\mu \right)\left(12f_p\right)\left(f_p\right)^2\right]`$

$`f_p\left(1f_p\right)\beta D_p\left[{\displaystyle \underset{\mu }{}}\left(Y_p^\mu \right)^2\left(1+n_B\left(E_\mu \right)f_p\right)\left(2\epsilon _p+E_\mu \right)\right]`$

$$f_p\left(1f_p\right)\beta D_p\left[\underset{\mu }{}\left(X_p^\mu \right)^2\left(n_B\left(E_\mu \right)+f_p\right)\left(2\epsilon _pE_\mu \right)\right]$$
$`(A21)`$
One should note here one thing about the direct use of the exact relation (A9) for the definition of the $`<N_p>`$. In reality the identity (A9) is only valid for the collective subspace (spanned by the non singly-occupied levels) or, in other words, for levels which have partial seniority $`s_k=0`$. But when we work in the grand canonical ensemble, seniority of the level $`s_k`$ is not a good quantum number, since averaging procedure over GCE mixes all seniorities. This fact is reflected in the above expression (A21) where we can see that level occupation numbers $`<N_p>=n_p+n_{\overline{p}}`$ due to thermal factors $`f_p`$ (Fermi-Dirac distribution) are not equal to $`<2P_p^+P_p>`$.

Before coming to a numerical example let us make a further comment. Knowing the occupation numbers as a function of the amplitudes $`X`$ and $`Y`$ as in (A21) or (32), a natural idea would be to use this to express the ground state energy entirely as a function of $`X`$ and $`Y`$. Then minimizing under the constraint of the normalization conditions (21) leads to an equation determining $`X,Y`$ amplitudes. Such a procedure has been proposed in the past by Jolos and Rybarska . However, already in our example where the occupation numbers are exactly known (at zero temperature) as a function of $`X,Y`$ amplitudes, we have checked that this leads to worse results than with the approach advanced in the main text. In other examples the functional $`n_p[X,Y]`$ may only be known approximately and then the minimization of the ground state leads to deteriorated results, as we also have checked. In the Green’s function approach where one first calculates excitation energies, i.e. energy differences, before calculating the ground state energy, such uncertainties like the precise knowledge of $`n_p[X,Y]`$ are minimized and the solution of the SCRPA equations is therefore much more stable.

Let us now compare the results based on the two different definitions of the occupation numbers (eq. (32) and (A21)). We will denote results with eq. (32) by TSCRPA and results obtained with the second definition (eq. (A21)) by TSCRPA1. All calculations are made for the system with 10 particles on 10 levels. Results for the correlation energy (37) as a function of the coupling constant $`G`$ are given in Tab. 1 and Tab. 2 for $`T=0`$ and $`T=1`$, respectively. In these tables we denote results obtained with (A21) and (26) by TSCRPA1-1 and the ones obtained with (A21) and (25) by TSCRPA1-2. We can see that TSCRPA systematically gives underestimated results with respect to the exact ones. Using definition (A21) we then can see that TSCRPA1-1 at zero temperature slightly overestimates the exact ground state energy (this corresponds to the results given in ) while TSCRPA1-2 always underestimates it. At finite temperature all approximations (TSCRPA and TSCRPA1) give close results which underestimate the exact ones. In Tab. 3 and 4 we show excitation energies of the first addition mode for $`T=0`$ and $`T=1`$. In both cases, TSCRPA and TSCRPA1, excitation energies grow with increasing $`G`$ for zero and nonzero temperatures reproducing quite well the exact results. We can see that for the calculated quantities (correlation and excitation energies) TSCRPA1-1 is slightly closer to the exact results than TSCRPA.

Table 1. Correlation energy as function of $`G`$ at zero temperature.

| G | Exact | TSCRPA1-1 | TSCRPA1-2 | TSCRPA |
| --- | --- | --- | --- | --- |
| 0.10 | -0.036 | -0.037 | -0.036 | -0.036 |
| 0.20 | -0.167 | -0.169 | -0.153 | -0.159 |
| 0.30 | -0.435 | -0.445 | -0.342 | -0.379 |
| 0.33 | -0.551 | -0.564 | -0.408 | -0.461 |
| 0.34 | -0.593 | -0.608 | -0.430 | -0.489 |
| 0.35 | -0.638 | -0.654 | -0.453 | -0.518 |
| 0.36 | -0.685 | -0.702 | -0.476 | -0.548 |
| 0.40 | -0.896 | -0.917 | -0.569 | -0.670 |

Table 2. Correlation energy as function of $`G`$ at $`T=1`$.

| G | Exact | TSCRPA1-1 | TSCRPA1-2 | TSCRPA |
| --- | --- | --- | --- | --- |
| 0.10 | -0.030 | -0.029 | -0.018 | -0.029 |
| 0.20 | -0.142 | -0.126 | -0.082 | -0.122 |
| 0.30 | -0.372 | -0.299 | -0.208 | -0.295 |
| 0.33 | -0.476 | -0.367 | -0.259 | -0.365 |
| 0.34 | -0.515 | -0.392 | -0.278 | -0.391 |
| 0.35 | -0.556 | -0.418 | -0.297 | -0.417 |
| 0.36 | -0.600 | -0.444 | -0.317 | -0.445 |
| 0.40 | -0.803 | -0.559 | -0.403 | -0.569 |

Table 3. Excitation energy of the first addition mode as function of $`G`$ at $`T=0`$.

| G | Exact | TSCRPA1 | TSCRPA |
| --- | --- | --- | --- |
| 0.10 | 1.001 | 1.001 | 1.001 |
| 0.20 | 1.005 | 1.012 | 1.012 |
| 0.30 | 1.014 | 1.049 | 1.053 |
| 0.33 | 1.018 | 1.068 | 1.074 |
| 0.34 | 1.019 | 1.075 | 1.081 |
| 0.35 | 1.022 | 1.082 | 1.089 |
| 0.36 | 1.023 | 1.089 | 1.098 |
| 0.40 | 1.031 | 1.123 | 1.136 |

Table 4. Excitation energy of the first addition mode as function of $`G`$ at $`T=1.0`$.

| G | TSCRPA1 | TSCRPA |
| --- | --- | --- |
| 0.10 | 1.017 | 1.030 |
| 0.20 | 1.075 | 1.126 |
| 0.30 | 1.175 | 1.281 |
| 0.40 | 1.299 | 1.476 |
| 0.41 | 1.312 | 1.497 |
| 0.42 | 1.324 | 1.518 |
| 0.43 | 1.337 | 1.539 |
| 0.44 | 1.349 | 1.560 |
| 0.45 | 1.362 | 1.581 |