Intrinsic average particle number in the HFB+VAPNP method

The HFB state $\vert\Phi\rangle$ is a linear combination of eigenstates $\vert N\rangle$ of the particle-number operator, i.e.,

$$\vert\Phi \rangle = \sum_{N} a_N \vert N \rangle ,$$ (80)

where

$$\vert N\rangle = \frac{P^N \vert\Phi\rangle}{\sqrt{\langle\Phi\vert P^N\vert\Phi\rangle}},$$ (81)

and $\hat{N}\vert N\rangle=N\vert N\rangle$. The HFB+VAPNP method is based on the variation of the projected energy (26), which is the average value of the Hamiltonian on the state $\vert N\rangle$, $E^N=\langle N\vert\hat{H}\vert N\rangle$. Obviously, the projected energy does not depend on amplitudes $a_N$, although the intrinsic average number of particles,

$$\bar{N} = \langle\Phi\vert\hat{N}\vert\Phi\rangle = \sum_{N} \vert a_N\vert^2N = \rho ,$$ (82)

does depend on $a_N$.

The HFB+VAPNP variational procedure gives, in principle, the same value of the projected energy independently of the value of $\bar{N}$. This independence can, however, be subject to numerical instabilities whenever the amplitude $a_N$, corresponding to the projection on $N$ particles, is small. Therefore, for practical reasons, one is interested in keeping the average number of particles $\bar{N}$ as close as possible to $N$, which guarantees that amplitude $a_N$ is as large as possible.

In the standard HFB equations (21), the average number of particles is kept equal to $N$ by adjusting the Lagrange multiplier $\lambda$. However, in the HFB+VAPNP approach, $\lambda$ does not appear in the variational equations (65), because the variation of the constant term $\lambda\langle\Phi\vert\hat{N}P^N\vert\Phi\rangle/\langle\Phi\vert P^N\vert\Phi\rangle=\lambda N$ equals zero. Therefore, the HFB+VAPNP equations (65) do not allow for adjusting the average particle number $\bar{N}$, which, during the iteration procedure, may become vary different from $N$. Moreover, such uncontrolled changes of $\bar{N}$ from one iteration to another may preclude reaching the stable self-consistent solution.

In order to cope with these problems, one can artificially reintroduce a constant $\mu$, analogous to the Fermi energy $\lambda$, into the HFB+VAPNP equations (65), i.e.,

$$\left( \begin{array}{cc} h^{N}-\mu & \tilde{h}^{N} \\ \tilde{h}^{... ...( \begin{array}{c} \varphi^{N}_{1,k} \\ \varphi^{N}_{2,k} \end{array} \right) ,$$ (83)

provided it is equal to zero once the convergence is achieved. With this modification, the iterations proceed as follows. Suppose that at a given iteration, condition $\bar{N}$=$N$ is fulfilled. Then, no readjustment of $\bar{N}$ is necessary, and in the next iteration one proceeds with $\mu$=0.

Had such an ideal situation continued till the end, a nonzero value of $\mu$ would have never appeared, and the required solution would have been found. In practice, this situation never happens, and at some iteration one finds that Tr$\rho$, i.e., the sum of norms of the second components $\varphi^N_{2,nk}$ (68) of the HFB+VAPNP wave functions, is larger (smaller) than $N$. In such a case, in the next iteration one uses a slightly negative (positive) value of $\mu$, which decreases (increases) the norms of $\varphi^N_{2,nk}$, and decreases (increases) the average particle number in the next iteration. Since $\mu$ acts in exactly the same way as the Fermi energy does within the standard HFB method, the well-established algorithms of readjusting $\lambda$ can be used. Moreover, as soon as the iteration procedure starts to converge, the non-zero values of $\mu$ cease to be needed, and thus $\mu$ naturally converges to zero, as required. In practice, we find that the above algorithm is very useful, and it provides the same converged solution with any value of $\bar{N}=N\pm\Delta{N}$, $\Delta{N}$ being a small integer.