format: text

text:
# Correlation energy of the pairing Hamiltonian

## I Introduction

An important goal of nuclear theory is to predict nuclear binding energies. The mean-field approximation certainly provides a good starting point, but the correlation energy associated with nearly degenerate configurations can be of the order of several MeV. Since one would like an order of magnitude better accuracy for global modeling, the correlations must be treated with some care. The correlations associated with broken mean-field symmetries, namely center-of-mass localization, deformations, and pairing, are especially important and need to be singled out for special attention. The study of such correlation effects has attracted much interest, see Refs.  for recent citations. Popular computational methods ***One may also ask whether simplified exact solutions may be of use for realistic applications. The Richardson solution for a pairing Hamiltonian is useful in this respect, but it requires a state-independent pairing interaction, which may not be realistic enough. That solution has recently been generalized to separable pairing interactions, but it is not yet clear whether the computational simplicity is preserved. include the generator coordinate method (GCM), variation-with-projection methods (e.g., VAP) , and the RPA .

In previous work, we have advocated the RPA approach for deformations and for center-of-mass localization. The correlation energy is calculated from the RPA formula

$$E_{corr}=\frac{1}{2}\left(\underset{i}{}\mathrm{}\omega _iTr(A)\right),$$
(1)
where $`\omega _i`$ is the frequency of the RPA phonon and $`A`$ is the $`A`$ matrix in the RPA equations. This method requires solving the RPA equations for each nucleus, which is computationally easy if the interaction is assumed to be separable. On other other hand, methods such as VAP and GCM are rather complicated to apply. One should also remember that a global theory requires a method that is not only systematic but also preserves the relative computational simplicity of the mean field approximation. In this respect, the Lipkin-Nogami method is a popular candidate for a computationally easy method to go beyond the mean field pairing theory .

In this paper, we continue our study of the accuracy of the RPA correlation formula (1), going on to pairing correlations. There is a long history of the RPA treatment of pairing correlations. An early study by Bang and Krumlinde  showed that the RPA formula reproduces the exact correlation energy rather well in a schematic model. Kyotoku et al.  compared the leading $`1/N`$ behavior of different methods including the Lipkin-Nogami method. They found that only RPA gave the exact coefficient in the condensed phase. The RPA method has in fact been used in realistic models of deformed nuclei. The RPA correlation in the normal phase was studied in Ref.  using the self-consistent version of RPA. But perhaps surprisingly in view of the large literature, we have not found any studies that specifically compare the RPA with the computationally attractive alternative methods, testing the behavior across shell closures and for odd $`N`$ systems. It is important that the method should not introduce spurious discontinuities when the mean field solution changes character; otherwise separation energies could not be reliably calculated .

To investigate the applicability of Eq. (1) and study the effects of symmetry breaking restoration, in this paper we employ a well-known schematic two-level pairing model, and compare its exact solutions to the approximations that are computationally attractive. The details of the RPA correlation energy are given in the next section. They include both the quasi-particle RPA extension of the BCS theory and the RPA for the pairing Hamiltonian when the strength is too weak for the mean field approximation to support the BCS solution (the “pairing-vibration” regime). The specific application of the various methods to the two-level model is given in the following section III.

## II Pairing Hamiltonian and the RPA

The pairing Hamiltonian and the BCS solution are well-known and we just summarize the equations to confirm the standard notation. In this paper we consider a Hamiltonian with an arbitrary single-particle term but a pairing interaction whose strength is state-independent,

$$H=\underset{j}{}ϵ_j\widehat{N}_jG\underset{j,j^{}}{}A_j^{}A_j^{},$$
(2)
where $`ϵ_j`$ is a single-particle energy. The number operator $`\widehat{N}_j`$ for the shell $`j`$ is given by

$$\widehat{N}_j=\underset{m>0}{}\left(a_{jm}^{}a_{jm}+a_{j\overline{m}}^{}a_{j\overline{m}}\right).$$
(3)
Here, $`a_{j\overline{m}}`$ is the annihilation operator for the time-reversal state and is given by $`a_{j\overline{m}}=()^{jm}a_{jm}`$. The pairing operator $`A_j^{}`$ in the Hamiltonian (2) is given by

$$A_j^{}=\underset{m>0}{}a_{jm}^{}a_{j\overline{m}}^{},$$
(4)
and $`A_j`$ is the Hermitian conjugate of $`A_j^{}`$. The exact solutions of the Hamiltonian were obtained a long time ago by Richardson and Sherman. In Appendix A, we solve the Richardson equation for a simple, but non-trivial case where there are two pairs in a single $`j`$-shell.

The BCS theory is the mean field solution to Eq. (2). We remind the reader that it is derived variationally from a ground state wave function of the form

$$|BCS=\underset{j,m>0}{}(u_j+v_ja_{jm}^{}a_{j\overline{m}}^{})|,$$
(5)
where $`v_j`$ is the usual pair occupation amplitude and $`u_j`$ satisfies $`u_j^2+v_j^2=1`$. This wave function is appropriate for even-$`N`$ systems, $`N`$ being the number of particles in a system. The extension to odd-$`N`$ systems is given at the end of this subsection. Important quantities in the variational solution are the chemical potential $`\lambda `$ and the pairing gap $`\mathrm{\Delta }=G_j\mathrm{\Omega }_ju_jv_j`$, $`\mathrm{\Omega }_j=(2j+1)/2`$ being the pair degeneracy of the $`j`$-shell. When the pairing strength $`G`$ is small, the variational equations only have the trivial solution, $`v_j,u_j=0`$ or 1 and $`\mathrm{\Delta }=0`$. The ground state wave function is thus given by $`|HF=_{j,m}a_{jm}^{}|`$ and the ground state energy is obtained as

$$E_{HF}=HF|H|HF=\underset{j:occupied}{}\mathrm{\Omega }_j(2ϵ_jG).$$
(6)
When the strength of the pairing interaction $`G`$ is larger than some critical value $`G_{crit}`$, the variational equations have a non-trivial solution $`\mathrm{\Delta }0`$ and the ground state energy is given by

$$E_{BCS}=BCS|H|BCS=2\underset{j}{}\mathrm{\Omega }_jϵ_jv_j^2\frac{\mathrm{\Delta }^2}{G}G\underset{j}{}\mathrm{\Omega }_jv_j^4.$$
(7)
For odd $`N`$ systems, one of the particles in a system does not form a pair and blocks a level. When the $`k`$-th level is blocked, the HF energy, the BCS energy, and the pairing gap are modified to

$`E_{HF}`$ $`=`$ $`{\displaystyle \underset{j:occupied}{}}\stackrel{~}{\mathrm{\Omega }}_j(2ϵ_jG)+ϵ_k,`$ (8)
$`E_{BCS}`$ $`=`$ $`2{\displaystyle \underset{j}{}}\stackrel{~}{\mathrm{\Omega }}_jϵ_jv_j^2{\displaystyle \frac{\mathrm{\Delta }^2}{G}}G{\displaystyle \underset{j}{}}\stackrel{~}{\mathrm{\Omega }}_jv_j^4+ϵ_k,`$ (9)
$`\mathrm{\Delta }`$ $`=`$ $`G{\displaystyle \underset{j}{}}\stackrel{~}{\mathrm{\Omega }}_ju_jv_j,`$ (10)
respectively. Here, $`\stackrel{~}{\mathrm{\Omega }}_j=\mathrm{\Omega }_j`$ for $`jk`$, and $`\stackrel{~}{\mathrm{\Omega }}_k=\mathrm{\Omega }_k\delta _{N,odd}`$. The chemical potential $`\lambda `$ is determined so that

$$\underset{j}{}\stackrel{~}{\mathrm{\Omega }}_jv_j^2+1=N,$$
(11)
is satisfied. This formalism is referred to as the blocked-BCS theory.

### A Random phase approximation

We next introduce the random phase approximation to compute the correlation energy. Let us first consider the RPA in the superfluid phase (QRPA). We refer to Refs.  for the formulation. The result is the well-known RPA matrix equation

$$\left(\begin{array}{cc}A& B\\ B& A\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right)=\mathrm{}\omega \left(\begin{array}{c}X\\ Y\end{array}\right),$$
(12)
where the matrices $`A`$ and $`B`$ are given by

$`A_{ij}`$ $`=`$ $`2E_i\delta _{i,j}G\sqrt{\stackrel{~}{\mathrm{\Omega }}_i}\sqrt{\stackrel{~}{\mathrm{\Omega }}_j}(u_i^2u_j^2+v_i^2v_j^2),`$ (13)
$`B_{ij}`$ $`=`$ $`G\sqrt{\stackrel{~}{\mathrm{\Omega }}_i}\sqrt{\stackrel{~}{\mathrm{\Omega }}_j}(u_i^2v_j^2+v_i^2u_j^2),`$ (14)
respectively. Here, $`\stackrel{~}{\mathrm{\Omega }}_j`$ is defined in the previous subsection, and $`E_j=\sqrt{(ϵ_j\lambda Gv_j^2)^2+\mathrm{\Delta }^2}`$ is the quasi-particle energy. The RPA excitation operator $`Q^{}`$ is given in terms of quasi-particle operators

$`\alpha _{jm}^{}`$ $`=`$ $`u_ja_{jm}^{}v_ja_{j\overline{m}},`$ (15)
$`\alpha _{j\overline{m}}^{}`$ $`=`$ $`u_ja_{j\overline{m}}^{}+v_ja_{jm},`$ (16)
as

$$Q^{}=\underset{j}{}\left(X_j\underset{m}{}\alpha _{jm}^{}\alpha _{j\overline{m}}^{}Y_j\underset{m}{}\alpha _{j\overline{m}}\alpha _{jm}\right)/\sqrt{\stackrel{~}{\mathrm{\Omega }}_j}.$$
(17)
The QRPA correlation energy is given by Eq. (1) with the $`A`$ matrix given by Eq. (13). Note that the second term on the right hand side of Eq. (1) is not $`_i2E_i`$ but $`Tr(A)`$. This point was overlooked in the integral approaches to the correlation energy in Refs. .

We next consider the RPA in a normal-fluid phase (pp-RPA), which describes pairing vibrations. The equation for the pp-RPA can be obtained from the QRPA equation by setting $`\mathrm{\Delta }=0`$ and $`v_h=u_p=1`$, where $`p`$ and $`h`$ denote particle and hole states, respectively. As for the chemical potential $`\lambda `$, there is no definite value for it in the normal fluid phase since it can be anywhere between the highest occupied and the lowest unoccupied levels. Notice, however, that the pairing vibration describes the ground state in the $`N\pm 2`$ systems and the role of the chemical potential is therefore just to shift the energies by an amount $`\pm 2\lambda `$. After removing these trivial energy shifts, one obtains

$`A_{pp^{}}`$ $`=`$ $`2ϵ_p\delta _{p,p^{}}G\sqrt{\stackrel{~}{\mathrm{\Omega }}_p}\sqrt{\stackrel{~}{\mathrm{\Omega }}_p^{}},`$ (18)
$`A_{hh^{}}`$ $`=`$ $`2(ϵ_hG)\delta _{h,h^{}}G\sqrt{\stackrel{~}{\mathrm{\Omega }}_h}\sqrt{\stackrel{~}{\mathrm{\Omega }}_h^{}},`$ (19)
$`B_{ph}`$ $`=`$ $`G\sqrt{\stackrel{~}{\mathrm{\Omega }}_p}\sqrt{\stackrel{~}{\mathrm{\Omega }}_h},`$ (20)
$`A_{ph}`$ $`=`$ $`B_{pp^{}}=B_{hh^{}}=0.`$ (21)
This $`2(N_p+N_h)\times 2(N_p+N_h)`$ dimensional matrix equation, where $`N_p`$ is the number of unoccupied shells and $`N_h`$ is that of occupied shells, can be decoupled into two $`(N_p+N_h)\times (N_p+N_h)`$ matrix equations as

$$\left(\begin{array}{cc}A& B\\ B& C\end{array}\right)\left(\begin{array}{c}X^a\\ Y^a\end{array}\right)=\mathrm{}\omega _a\left(\begin{array}{c}X^a\\ Y^a\end{array}\right),$$
(22)
and

$$\left(\begin{array}{cc}C& B\\ B& A\end{array}\right)\left(\begin{array}{c}X^r\\ Y^r\end{array}\right)=\mathrm{}\omega _r\left(\begin{array}{c}X^r\\ Y^r\end{array}\right).$$
(23)
Here, the matrix $`C`$ is defined as $`C_{hh^{}}=A_{hh^{}},C_{ph}=0`$, and RPA amplitudes are given by $`X_p^a=X_p,Y_h^a=Y_h,X_h^r=X_h`$, and $`Y_p^r=Y_p`$. The first equation (22) describes the ground state of the $`N+2`$ system and is referred to as the addition mode, while the second equation (23) describes the $`N2`$ system and is called the removal mode. Noticing that Eq. (1) is a half of a sum of the difference between the RPA and the Tamm-Dancoff approximation (TDA) frequencies for each mode of excitation , we find the correlation energy for the addition mode to be

$$E_{corr}^a=\frac{1}{2}\left(\underset{i}{}\mathrm{}\omega _{ai}Tr(A)\right),$$
(24)
while that for the removal mode is

$$E_{corr}^r=\frac{1}{2}\left(\underset{i}{}\mathrm{}\omega _{ri}Tr(C)\right).$$
(25)
The total correlation energy is the sum of these, $`E_{corr}=E_{corr}^a+E_{corr}^r`$.

Typically, one finds that an RPA mode goes to zero frequency at the mean-field phase transition. This is not the case for Eqs. (22) and (23), thus giving a discontinuity in the frequencies at the transition. However, as we have mentioned above, this is an artifact of an awkward choice of the chemical potential $`\lambda `$, and it will be the case that the sum of an addition and removal mode goes to zero.

### B Lipkin-Nogami method

An alternative way to restore the broken gauge symmetry of the BCS approximation is to carry out the number projection of the BCS wave function. The Lipkin-Nogami method provides an approximate way for number projection. It was first invented by Lipkin , and was developed by Nogami and his collaborators . Because of its relative simplicity, it has been widely applied .

In the Lipkin-Nogami method, the expectation value $`H\lambda \widehat{N}\lambda _2\widehat{N}^2_{BCS}`$ is varied. The resultant equations to be solved are given by

$`{\displaystyle \frac{2}{G}}={\displaystyle \underset{j}{}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Omega }}_j}{\sqrt{\stackrel{~}{ϵ}_j+\mathrm{\Delta }^2}}},`$ (26)
$`2{\displaystyle \underset{j}{}}\stackrel{~}{\mathrm{\Omega }}_jv_j^2+\delta _{N,odd}=N,`$ (27)
$`\lambda _2={\displaystyle \frac{G}{4}}\left\{{\displaystyle \frac{\left(_i\stackrel{~}{\mathrm{\Omega }}_iu_i^3v_i\right)\left(_i\stackrel{~}{\mathrm{\Omega }}_iu_iv_i^3\right)_i\stackrel{~}{\mathrm{\Omega }}_iu_i^4v_i^4}{\left(_i\stackrel{~}{\mathrm{\Omega }}_iu_i^2v_i^2\right)^2_i\stackrel{~}{\mathrm{\Omega }}_iu_i^4v_i^4}}\right\},`$ (28)
where $`\stackrel{~}{ϵ}_i,v_i`$, and $`u_i`$ are defined as

$`\stackrel{~}{ϵ}_i`$ $`=`$ $`ϵ_i+(4\lambda _2G)v_i^2\lambda ,`$ (29)
$`v_i^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{\stackrel{~}{ϵ}_i}{\sqrt{\stackrel{~}{ϵ}_i^2+\mathrm{\Delta }^2}}}\right),`$ (30)
$`u_i^2`$ $`=`$ $`1v_i^2,`$ (31)
respectively. We use the blocked-Lipkin-Nogami prescription for odd $`N`$ systems. One of the characteristic features of the Lipkin-Nogami method is that the pairing gap $`\mathrm{\Delta }`$ has a finite value even in the weak $`G`$ limit where $`\mathrm{\Delta }`$ is zero in the BCS approximation. The ground state energy in the Lipkin-Nogami method is given by

$$E_{LN}=2\underset{j}{}\stackrel{~}{\mathrm{\Omega }}_jϵ_jv_j^2\frac{\mathrm{\Delta }^2}{G}G\underset{j}{}\stackrel{~}{\mathrm{\Omega }}_jv_j^44\lambda _2\underset{j}{}\stackrel{~}{\mathrm{\Omega }}_ju_j^2v_j^2+ϵ_k\delta _{N,odd}.$$
(32)
## III Two-level pairing model

We now apply the above equations to a schematic two-level model. This model was first introduced in Ref., and has been used in the literature to test several approximations . We label the lower and the higher levels 1 and 2, taking $`ϵ_1=ϵ/2,ϵ_2=ϵ/2`$. The Hamiltonian can be numerically diagonalised using the quasi-spin formalism . The basis states are denoted by $`|S_1S_{10};S_2S_{20}`$, where $`S_i`$ and $`S_{i0}`$ are defined by $`S_i=(\mathrm{\Omega }_i\nu _i)/2`$ and $`S_{i0}=(N_i\mathrm{\Omega }_i)/2`$, respectively. The latter takes a value of $`S_i,S_i+1,\mathrm{},S_i`$. $`\nu _i`$ is the seniority quantum number. For the ground state of even-$`N`$ systems, the total seniority $`\nu _1+\nu _2`$ is 0, while it is 1 for odd-$`N`$ systems. The matrix elements of the Hamiltonian read

$`S_1^{}S_{10}^{};S_2^{}S_{20}^{}|H|S_1S_{10};S_2S_{20}=\delta _{S_1,S_1^{}}\delta _{S_2,S_2^{}}\{ϵ(S_{20}S_{10}{\displaystyle \frac{\mathrm{\Omega }_1\mathrm{\Omega }_2}{2}})\delta _{S_{10},S_{10}^{}}\delta _{S_{20},S_{20}^{}}`$ (33)
$`G\left(S_1(S_1+1)S_{10}(S_{10}+1)+S_2(S_2+1)S_{20}(S_{20}+1)\right)\delta _{S_{10},S_{10}^{}}\delta _{S_{20},S_{20}^{}}`$ (34)
$`G\sqrt{(S_1(S_1+1)S_{10}(S_{10}1)}\sqrt{S_2(S_2+1)S_{20}(S_{20}+1)}\delta _{S_{10},S_{10}^{}+1}\delta _{S_{20},S_{20}^{}1}`$ (35)
$`G\sqrt{(S_1(S_1+1)S_{10}(S_{10}+1)}\sqrt{S_2(S_2+1)S_{20}(S_{20}1)}\delta _{S_{10},S_{10}^{}1}\delta _{S_{20},S_{20}^{}+1}\}.`$ (36)
We first consider a symmetric two-level problem i.e., $`\mathrm{\Omega }_1=\mathrm{\Omega }_2=\mathrm{\Omega }`$, and assume the number of fermionic particle is $`N=2\mathrm{\Omega }`$. Equations for several approximations can be solved analytically for such a system. The gap equation in the BCS theory leads to a pairing gap of

$$\mathrm{\Delta }=\sqrt{G^2\mathrm{\Omega }^2\frac{\stackrel{~}{ϵ}^2}{4}},$$
(37)
together with

$`v_1^2`$ $`=`$ $`u_2^2={\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{\stackrel{~}{ϵ}}{2G\mathrm{\Omega }}}\right)v^2,`$ (38)
$`u_1^2`$ $`=`$ $`v_2^2={\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{\stackrel{~}{ϵ}}{2G\mathrm{\Omega }}}\right)u^2.`$ (39)
$`\lambda `$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{ϵ}ϵ}{2}}Gv^2,`$ (40)
where $`\stackrel{~}{ϵ}`$ is defined as $`\stackrel{~}{ϵ}=2\mathrm{\Omega }ϵ/(2\mathrm{\Omega }1)`$.

From Eq. (37), one can find the critical strength of the phase transition to be $`G_{crit}=ϵ/(2\mathrm{\Omega }1)`$. For $`G`$ larger than $`G_{crit}`$, the system is in the superfluid phase, and the matrices $`A`$ and $`B`$ for the QRPA equation are given by (see Eqs. (13) and (14))

$`A_{11}`$ $`=`$ $`A_{22}=G\mathrm{\Omega }{\displaystyle \frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }^2}}+{\displaystyle \frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}},`$ (41)
$`A_{12}`$ $`=`$ $`A_{21}={\displaystyle \frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}},`$ (42)
$`B_{11}`$ $`=`$ $`B_{22}={\displaystyle \frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }^2}}+{\displaystyle \frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}},`$ (43)
$`B_{12}`$ $`=`$ $`B_{21}=G\mathrm{\Omega }{\displaystyle \frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}},`$ (44)
The solutions of the QRPA equation are $`\mathrm{}\omega =0`$ and $`\sqrt{4\mathrm{\Delta }^22\mathrm{\Delta }^2/\mathrm{\Omega }}`$, from which the correlation energy is computed as

$$E_{corr}=\sqrt{\mathrm{\Delta }^2\mathrm{\Delta }^2/2\mathrm{\Omega }}G\mathrm{\Omega }+\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }^2}\frac{\mathrm{\Delta }^2}{2G\mathrm{\Omega }}.$$
(45)
When the strength of the pairing interaction $`G`$ is smaller than $`G_{crit}`$, the system is in the normal fluid phase. The mean field energy is then evaluated as $`E_{HF}=ϵ\mathrm{\Omega }G\mathrm{\Omega }`$. For the present two-level model, the $`A,B`$, and $`C`$ matrices for the pp-RPA are just numbers and are given by $`A=ϵG\mathrm{\Omega }`$, $`B=G\mathrm{\Omega }`$, and $`C=ϵG\mathrm{\Omega }+2G`$, respectively. The frequencies for the addition and the removal modes are then found to be

$`\mathrm{}\omega _a`$ $`=`$ $`G+\sqrt{ϵ+G}\sqrt{ϵ+G2G\mathrm{\Omega }},`$ (46)
$`\mathrm{}\omega _r`$ $`=`$ $`G+\sqrt{ϵ+G}\sqrt{ϵ+G2G\mathrm{\Omega }},`$ (47)
respectively. The total correlation energy is thus obtained as

$$E_{corr}=\sqrt{ϵ+G}\sqrt{ϵ+G2G\mathrm{\Omega }}(ϵG\mathrm{\Omega }+G).$$
(48)
Figure 1 shows the RPA frequencies for each modes of excitation as a function of $`G`$, for $`\mathrm{\Omega }`$ = 8. As we noted in the previous section, we see that because of the chemical potential the pp-RPA frequencies do not match with the QRPA frequencies at the critical point of phase transition from a normal-fluid to a superfluid phases. Figure 2 compares the ground state energy obtained by several methods. The solid line is the exact solution obtained by numerically diagonalizing the Hamiltonian. The dashed line is the ground state energy in the mean field BCS approximation. It considerably deviates from the exact solution through the entire range of $`G`$ shown in the figure. The dot-dashed line takes into account the RPA correlation energy in addition to the mean field energy. It reproduces very well the exact solutions, except in the vicinity of the critical point of the phase transition. Around the critical point, one would need to compute the correlation energy with some care, using e.g., the self-consistent RPA discussed in Ref.  which removes the cusp behaviour of the RPA frequency around the critical strength.

It is interesting to compare the present mean field plus RPA approach with the Lipkin-Nogami method. The equations for the Lipkin-Nogami method for the two-level model were solved in Ref. . The results are given by

$`\mathrm{\Delta }^2`$ $`=`$ $`(G\mathrm{\Omega })^2(1\stackrel{~}{\kappa }^2),`$ (49)
$`v_1^2`$ $`=`$ $`u_2^2=(1+\stackrel{~}{\kappa })/2,`$ (50)
$`u_1^2`$ $`=`$ $`v_2^2=(1\stackrel{~}{\kappa })/2,`$ (51)
$`4\lambda _2G`$ $`=`$ $`{\displaystyle \frac{2G\mathrm{\Omega }\stackrel{~}{\kappa }^2}{(2\mathrm{\Omega }1)(1\stackrel{~}{\kappa }^2)}},`$ (52)
where $`\stackrel{~}{\kappa }`$ is the physical solution, which satisfies $`0\stackrel{~}{\kappa }1`$, of an equation

$$2(1\mathrm{\Omega })\stackrel{~}{\kappa }^3+(2\mathrm{\Omega }1)\kappa \stackrel{~}{\kappa }^2+(2\mathrm{\Omega }1)\stackrel{~}{\kappa }(2\mathrm{\Omega }1)\kappa =0,$$
(53)
and $`\kappa `$ is defined as $`ϵ/2G\mathrm{\Omega }`$. The ground state energy in the Lipkin-Nogami method given by Eq. (32) is denoted by the thin solid line in Fig. 2. Although it reproduces the exact results for large values of $`G`$, it deviates significantly from them for small values. This behaviour is consistent with the numerical observation in Ref.  as well as the result of Ref. where it was found that the Lipkin-Nogami method is only correct in the limit of a strong pairing force. In marked contrast, the ground state energy in the BCS plus QRPA coincides with the exact solution at the leading order of an expansion in 1/$`\mathrm{\Omega }`$ for any strength of the pairing interaction as long as it is larger than the critical value .

From Fig. 2, it is unclear whether the Lipkin-Nogami method or RPA are more accurate for purposes of a global theory of binding. So we now consider a more realistic situation, varying the particle number $`N`$ rather than the interaction strength $`G`$. We consider the paring energy in oxygen isotopes, taking the neutron 1p and 2s-1d shells as the lower and higher levels of the two-level model. The pair degeneracy $`\mathrm{\Omega }`$ thus reads $`\mathrm{\Omega }_1=3`$ and $`\mathrm{\Omega }_2=6`$, and the number of particle in a system is given by $`N=A82`$ for the AO nucleus. We assume that the energy difference between the two levels $`ϵ`$ is given by $`ϵ=41A^{1/3}`$ and the pairing strength $`G=23/A`$. The upper panel of Fig. 3 shows the ground state energy as a function of $`A`$. In order to match with the experimental data for the 16O nucleus, we have added a constant $`72.8`$ MeV to the Hamiltonian for all the isotopes. The exact solutions are denoted by the filled circles. The deviation from the BCS approximation (the dashed line) is around 2 MeV for even $`A`$ systems and it is around 1.2 MeV for odd $`A`$ systems. This value varies within about 0.5 MeV along the isotopes and shows relatively strong $`A`$ dependence. One can notice that the RPA approach (the dot-dashed line) reproduces quite well the exact solutions. On the contrary, the Lipkin-Nogami approach (the thin solid line) is much less satisfactory and shows a different $`A`$ dependence from the exact results. The pairing gaps $`\mathrm{\Delta }`$ in the BCS approximation and in the Lipkin-Nogami method are shown separately in the lower panel of Fig. 3. For the Lipkin-Nogami method, we show $`\mathrm{\Delta }+\lambda _2`$, which is to be compared with experimental data . The closed shell nucleus 16O and its neighbour nuclei 15,17O have a zero pairing gap in the BCS approximation, and the Lipkin-Nogami method does not work well for these nuclei, as can be casted in Fig. 2. On the other hand, the RPA approach reproduces the correct $`A`$ dependence of the binding energy. Evidently, the RPA formula provides a better method to compute correlation energies than the Lipkin-Nogami method, especially for shell closures.

## IV Summary

Returning to our initial motivation, we seek a computationally tractable way to include pairing effects in a global model of nuclear binding energies, going beyond the BCS theory. Two attractive possibilities are the Lipkin-Nogami method and the RPA, in particular if the pairing interaction has a separable form. In this paper, we used a solvable two-level pairing Hamiltonian to show that the RPA formula for the correlation energy reproduces well the exact solutions both in a normal-fluid and a superfluid phases. On the contrary, the Lipkin-Nogami method is considerably less accurate for weak pairing, and therefore is not suitable in transition regions. As a consequence, the Lipkin-Nogami method fails to reproduce the correct mass number dependence of the binding energy around a shell closure. The correlation energy associated with the number fluctuation for the neutron mode in O isotopes was shown to be order of 2 MeV and has a relatively strong mass dependence. Although the correlation energy is small compared with the absolute value of the ground state energy, this suggests that including the correlation energy in the RPA provides a promising way to develop a better microscopic systematic theory for nuclear binding energies.

## Acknowledgments

The authors thank P.-G. Reinhard and W. Nazarewicz for useful discussions. This work was supported by the U.S. Department of Energy under Grant DOE-ER-40561.

## A Richardson solution

For a system where the number of pairs is $`N_{pair}`$, i.e., a 2$`N_{pair}`$-fermion system, the ground state energy of the Hamiltonian (2) is given by

$$E_{gs}=\underset{\lambda =1}{\overset{N_{pair}}{}}z_\lambda ,$$
(A1)
where $`z_\lambda `$ are $`N_{pair}`$ solutions of $`N_{pair}`$ coupled equations given by

$$\underset{j}{}\frac{\mathrm{\Omega }_j}{2ϵ_jz_\lambda }\stackrel{}{}_\lambda ^{}\frac{2}{z_\lambda ^{}z_\lambda }=\frac{1}{G}.$$
(A2)
The prime in the summation of $`\lambda ^{}`$ means to take only those $`\lambda ^{}`$ which are different from $`\lambda `$, and $`\mathrm{\Omega }_j=(2j+1)/2`$ is the pair degeneracy of the $`j`$-shell. In general, $`z_\lambda `$ are complex.

For a two pair system ($`N_{pair}`$=2) in a single $`j`$-shell, setting $`ϵ_j=0`$ and $`\mathrm{\Omega }_j=\mathrm{\Omega }`$, the Richardson equation reads

$`{\displaystyle \frac{\mathrm{\Omega }}{z_1}}{\displaystyle \frac{2}{z_2z_1}}`$ $`=`$ $`{\displaystyle \frac{1}{G}}`$ (A3)
$`{\displaystyle \frac{\mathrm{\Omega }}{z_2}}{\displaystyle \frac{2}{z_1z_2}}`$ $`=`$ $`{\displaystyle \frac{1}{G}}.`$ (A4)
The solutions of these equations are found to be

$$z=G(\mathrm{\Omega }1)\pm i\sqrt{G^2(\mathrm{\Omega }1)},$$
(A5)
from which the ground state energy reads $`E_{gs}=2G(\mathrm{\Omega }1)`$. This result coincides with the solution obtained using the seniority scheme,

$$E(N_{pair})=G\left\{\frac{\mathrm{\Omega }}{2}\left(\frac{\mathrm{\Omega }}{2}+1\right)\frac{1}{4}(2N_{pair}\mathrm{\Omega })^2+\frac{1}{2}(2N_{pair}\mathrm{\Omega })\right\}.$$
(A6)


text:
# Will the Population of Humanity in the Future be Stabilized?

## I Introduction

The problem of mankind population growth is the one of the global problems concerning the development of the mankind and its future. Will the demographic explosion existing now at the mankind population of the whole world stop as already it has stopped in the most developed countries? Will the population of the Earth will be stabilized as it follows from S.Kapitza’s theory (see )near 12 billions, or it will continue its growth at slower rate? What are the driving parameters that govern the development of mankind (such as presents in non-linear open system (see)? Are these parameters genetically predetermined or can they be changed and controlled by means of human activity? Stabilization of the mankind population of earth in the future is a sad prospect for mankind, because the absence of the numerical increasing of any biological population almost always leads, early or late, to cessation of any development (the examples are many species of insects, e.g. termites, frozen in development for millions years). Hence, the appearance of more active biological species becomes in that case quite probable. These active form of life will dominate the mankind and may force it out from its present ecological niche. The aim of the present paper is to introduce a new parameter in the phenomenological theories of humankind growth describing the development of mankind population if regard humankind as a large non-linear multifractal system. Namely, the fractal dimension of the whole mankind number in arbitrary moment of time - the local fractional dimension $`d(t)`$ (the fractional dimension for whole population of mankind). The introducing of this parameter allows to receive, as a special cases, the results of theories (see ),-), and several new scenarios of the mankind’s future as well. Alongside with probable increasing of the mankind population, the scenario of ruin of the present civilization - diminution of its number down to zero - is shown to be possible. We note, that correct analysis of dynamics of multifractal sets (see) requires the introduction of the mathematical concept of fractional derivatives (see -), which allow to take into account partly the memory of system about the past (in this case it is the memory, that includes genetic memory of mankind about its past development).

## II Fractional Dimension Quantity of Mankind on Time axis

We assume that the dynamics of the human population can be described within the framework of fractal geometry concepts and mathematical formalism of fractional differential equations (see -). For this purpose let us consider the set of all people as a multifractal set $`N(t)`$ $`(0t\mathrm{})`$ consisting of $`N(t)`$ of elements at the given time $`t`$. Following the Kapitza theory (see ), assume for a certain interval of time the existence of a self-similarity of this set characterized by its fractal dimension and introduce a local fractional dimension (LFD) $`d(t)`$ which describes the fractal properties of the set at the time $`t`$. This local fractional dimension is determined by those variables and their functions (these variables will be defined later) , that are treated as the driving (control) parameters of the human community development. Among such parameters there can be parameters of genetic origin (for example, probably, density of the population in cities per unit urban area, etc.) and ”external” parameters (e.g., a possibility of supplying the mankind population of earth by necessary amounts of food, water and energy or quick climate changes or pollution of the environment and the atmosphere, etc.). We shall characterize an alteration of the mankind population $`N(t)`$ over a small time interval by the generalized fractional Riemann - Liouville derivative  (which coincides with the usual Riemann - Liouville derivative if $`d`$=const.)

$$D_{+,t}^{1+\nu (t)}N(t)=\frac{^\alpha }{t^\alpha }\underset{0}{\overset{t}{}}\frac{N(t^{})dt^{}}{\mathrm{\Gamma }(\alpha d(t^{}))(tt^{})^{d(t^{})\alpha +1}},$$
(1)
$`d(t)=1+\nu (t)>0`$ (2)
$$\nu (t)\nu (X_1(t),X_2(t)\mathrm{}X_i(t)),i=1,2,\mathrm{},\alpha =\{d\}+1$$
(3)
In (1) $`\nu (t)`$ is the fractional quantity and defines the differences between the derivatives of integer order and fractional derivatives (1) thus being the driving parameter for the growth of mankind as a whole, $`\{d\}`$ is equals to the integer part of $`d(t)>0`$ $`(\alpha 1d(t)<\alpha )`$ , $`\alpha =0`$ for $`d<0`$ , and the set $`X_i`$ are the control parameters determining all external and internal influences on the mankind population growth .

The explicit information about the function $`\nu `$ and, hence, about LFD can be obtained only after a careful investigations and processing the statistical data of different events, circumstances and situations impact on the development of a human population. The fractional derivative using, as defined in (1), instead of the integer first derivative allows to introduce and take into account an obvious thing - a certain kind of mankind’s memory of the past and memory about the development rates over the past years (integration with a weight function over all times till $`t`$ beginning from a fixed moment). It gives a way of considering of different parameters $`X_i`$ that influence at the mankind’s development by means of LFD’s dependence upon them. The present theory, as well as (), is a phenomenological theory and the exact definition of function $`\nu (t)`$ form is beyond its scope. We note that it has sense to consider $`\nu (t)1`$, apparently, only for times greater than $`T_1^{}`$ because before a certain time, introduced in (, -) (the time of demographic transition $`T=T_1^{}`$), the theories mentioned describe the empirical data about the number and progress of the mankind population quite well. Moreover, we shall restrict ourselves in this paper to analyzing growth of the mankind population for three special cases of $`\nu (t):\nu (t)=0,\nu (t)=1`$ and $`\nu (t)=1`$.

## III Foerster,Hoerner, Kapitza Theories

It is known, that growth of the population of earth from ancient times untill now is well approximated by an empirical relation suggested by Foerster Von H. (see ) and improved by H. von Xoerner (see )

$$N(t)=\frac{C}{T_1T^{}},$$
(4)
$$\frac{d(N(t))}{dt}=\frac{C}{(T_1T)^2}$$
(5)
$`C=2^.10^{11},T_1=2025`$

The reasonable generalization (4),(5) for future time (suitable for $`T>T_1^{}`$ that not gives as result an infinite value $`N(t)`$ for values $`T=T_1`$) was given in the phenomenological theory of the population S.Kapitza (see ) with the help of introducing of the mean people’s lifetime $`\tau `$ ($`\tau `$ =42 of years). In this theory the relation for $`N(T)`$ (with C’=1,86$`{}_{}{}^{.}10_{}^{11},T_1^{}=2007`$) takes place

$$\frac{N(t)}{dt}=\frac{C^{}}{(T_1^{}T)^2+\tau ^2}$$
(6)
From (6) the restriction of the mankind population of earth by quantity $`14^.10^9`$ follows. Unfortunately, the theory of the population S.Kapitza does not take into account neither exterior nor interior control parameters (even in simple form) basing only at the self-similarity of growth of a population of the people.

## IV Generalization of Foerster, Hoerner and Kapitza Theories for Multifractal Set of a Quantity of Mankind $`N(t)`$

Let us assume the hypothesis about a fractal nature of set $`N(t)`$ (maintaining assumption about selfsimilar the sets $`N(t)`$). In that case derivative on time $`\frac{}{t}`$ in equation (6) should be substitute for fractional derivative $`D_{+,t}^{1+\nu (t)}`$. This operation take into account the memory of mankind about the past development. The right part (6) must be changed too for including in (6) an influence of the fractal dimension. So, instead of (6), obtain an equation

$$D_{+,t}^{1+\nu (t)}N(t)=\frac{C^{}}{|T_1^{}T|^{2+\nu (t)}+\frac{2+\nu (t)}{2}\tau ^{2+\nu (t)}}$$
(7)
The equation (7) can be considered as the basic equation of the phenomenological theory of development mankind’s population offered in this paper. The selection of different functions of fractional corrections for $`\nu (t)`$ allows to estimate character of changing $`N(t)`$ as functions of time.

## V Forecasts of the Future Development of Mankind Population for Special Cases of the Fractional Dimension Choice

Some simple special cases forecasting of growth of the mankind’s population of earth are considered below at the basis of the equation (7). The meanings a fractal dimension, for simplicity , chosen as integer.

a. Case $`\nu (t)=0`$

At $`\nu (t)=0`$ the fractional derivative $`D`$ coincides with $`\frac{}{t}`$ derivative and the equation (7) coincides with equation (6) (hence, (7) includes S.Kapitza theory (see ) as a special case). It corresponds, probably, a compensation the positive and the negative control parameters $`X_i`$ (including exterior and ”interior” parameters), driving population of mankind.

b. Case $`\nu (t)1`$

The selection $`\nu (t)1`$ corresponds to a dominance of the negative tendencies in the future development of mankind and it is stipulated, for example, occurrence of irreversible changes in molecules DNA (owing to epidemic AIDS, etc.), irreversible cosmic cataclysmic (for example, drop on earth of meteorites of huge mass), impossibility for mankind to cope with negative factors of biogeozinos results in by effects of development of a technical civilization. In this case the equation (7) reduce to equation

$$N(t)=\frac{2C^{}}{2|T_1^{}T|+\tau }$$
(8)
From (8) the presence of the maximum of number of mankind follows at $`T=T^{}`$ (i.e. in 2007) and it is equal $`8,86^.10^9`$ (if $`\tau `$ is not changed). After transition through a maximum $`N(t)`$ the number of mankind decreases (if the scenario will not vary) and in 2107 year $`N(t)`$ will be equal $`1,54^.10^9`$. By the year 3007 the number of mankind will decrease down to $`182^.10^6`$, i.e. the mankind population of all earth will be equal to number of the people is occupying a dozen of large modern cities. The complete degenerating of mankind, as a result of decreasing of mankind population in the considered scenario and subsequent leaving mankind from biological arena may be expected through millions years ( not large time from biological point of view). By then population the mankind of earth will decrease up to several thousand.

c. Case $`\nu (t)=2.`$

The equation (7) for this case is transformed into an integral equation at maintenance of the general tendency to accumulation of the negative factors resulting in to negative value $`d(t)`$ . At negative values $`d(t)`$ the integral equation gives the prompt diminution of mankind population and extinction of existence of mankind as a species is follows. So, at $`d(t)1`$ equations (7) transforms to a

$$\underset{T^{}}{\overset{T}{}}N(t)𝑑t=C^{}$$
(9)
supposing the absence the mankind at the earth (in that case equation (9) has no solution for $`N(t)0`$ thou for $`\nu >2`$ solution of equation (7) is exist (so $`N(t)0`$ for $`\nu 2`$). The time interval (necessary for disappearance mankind) is determined in this case by time for which $`d(t)`$ will transfer from value equal to unity (status of mankind now) to value equal to negative unity. This time interval can be very short: from several years (cosmic cataclysmic) to about several centuries (virus pandemia with lethal change of a heredity, increasing of the mean temperature of earth on some degrees owing to throw out of carbonic gas etc.). It is necessary to pay attention on possibility of change of the negative scenario of mankind development stipulated by appearance and dominance of the positive control factors $`X_i`$ (including those factors that due to conscious activity of mankind). In that case the inevitability of diminution of mankind population and destruction of mankind are not inevitable.

d. case $`\nu =1.`$

We shall consider the optimistic scenario of change of mankind population. It relevant to a dominance of positive driving parameters: $`d(t)>1.`$ So, e.g. for $`d(t)=2`$ $`(\nu =1)`$ from (7) the equation follows

$$\frac{^2}{t^2}N(t)=\frac{C^{}}{|T^{}T|^3+1,5\tau ^3}$$
(10)
The precise solution of an initial value problem of this equation is unwieldy, so we shall note , that at $`T>>\tau `$ the quantity $`N(t)`$ will increase faster then first degree of time ($`N(t)(TT_1)`$). The mankind population is increase and it characterizes in this case by following form (if there will not be includes appearance of a factors of conscious mankind activity which change the scenario and restricts unlimited growth of population)

$$N(t)|^t\mathrm{}\frac{C^{}(TT_1^{})}{2,29\tau ^2}\mathrm{ln}\frac{(TT_1^{})}{\tau }$$
(11)
So, at 3000 years, at the rate of increasing of the population defined by (11) (with the allowances that the corrections to $`N(t)`$ of value $`ln[\frac{(TT_1^{})}{\tau }]`$ are dropped) population of mankind will increase up to 150 billions. That is improbable large, though, but it is may be not unreasonable because of future technical possibilities of mankind (probably, this number is the upper number for existence mankind population occupying the earth). For fractional values LFD $`d(t)`$, increasing or decreasing of a population of the mankind will be characterized by intermediate dependencies between the received for the whole values $`d(t)`$. In case of a periodic dependencies $`d(t)`$ from time the population of the world will change periodically depending on a concrete choice of $`d`$ and $`N(t)`$ and will not be a monotonous function of time.

The examples are considered allow to determine the interval of change of fractional dimension $`d(t)`$ in reason boundary for number of quantity mankind in future: $`1<d<2`$. The boundary values $`(d=1,d=2)`$ are result in or to ruin of mankind, or to so large mankind population that Earth may not endure. The last case must result in to change of scenario and it consist in the change in correlation between positive and negative components of control parameters $`X_i`$ towards increasing of influence of negative parameters.

Let the basic equation (7) is replaced by generalization of the S.Kapitza equation (6). Basic equation in that case reads

$$\frac{^{1+\nu (t)}}{t^{1+\nu (t)}}n(t)=K\mathrm{sin}^2\frac{n(t)}{K}=\frac{1}{K},$$
(12)
$`n={\displaystyle \frac{N(t)}{K}},t={\displaystyle \frac{TT_1^{}}{\tau }},K=\sqrt{C^{}\tau ^1}`$ (13)
It gives for $`N(t)`$ the qualitative effects analogous to results obtained from the equation (7). So, in this connection, a selection for describing the future increasing the population of mankind by equation (7), as the basis, or equation (12), containing, as well as the equation (7), driving parameters $`X_i`$ in fractional dimension $`d(t)`$, is not simple. As one of advantages of use of the equations (7) or (12) for describing demographic problems (with some of them the mankind already has confronted now) we shall stress an opportunity of insert and account in the theory many factors (such as , incurable illnesses, natural cataclysmic etc.) defining future of mankind as a result of influences the control parameters $`X_i`$ (included in the dimension $`d(t)`$ of fractal set for number’s distribution of the people in the time axis).

## VI Conclusions

1.The main purpose of this paper was to analyze possibility of introducing mankind’s population driving parameters $`X_i`$ in the phenomenological theories of the mankind population of the earth (, , ) by method of attributing to set of the people $`N(t)`$ the fractional dimension $`d(t)`$, depending from these parameters. At a choice $`d(t)=1`$ the numbered theories are a special case of this phenomenological theory. The consideration of examples with integer meanings of $`d(t)`$ is caused only by their mathematical simplicity and gives a reasonable meaning of fractal dimensions ($`1<d<2`$) for describing the time dependence of population of mankind.

2. In case when the interpretation the fractal properties of set of the homosapience given in this paper corresponds to a reality, ( more real cases correspond to fractional meanings $`d(t)`$) the future of mankind is not so mournful as in the case of the S.Kapitza theory (see ) and the exposition of change of number of mankind within the framework of the phenomenological theories of the population can be reduced to a selection of control parameters and filling them by the concrete contents.

3. The change of number of mankind (described in the framework a phenomenological theories of the population) can be adjusted by such choice of control parameters (and filling the fractal dimension $`d(t)`$ by the concrete contents of dependencies of them) at which the population of mankind will grow so slowly, that overpopulation and the problems connected with it will do not arise in the foreseen future. Last will allow the theory be more realistic for predicting and menaging the future growth population of mankind as one of biological species occupying our world.

4. The growth of mankind population regulation (included in the parameters $`X_i`$) will allow to avoid degenerating of mankind and to keep as much as long time the main ecological niche at Earth occupied by mankind. Last will give time for more realistic forecasting of the future mankind as one of biological kinds occupying our world.


text:
# THE BEST UNBIASED ESTIMATOR FOR THE CMB ANGULAR BISPECTRUM

## References