Particle-Number Projection After Variation

Introducing the particle-number projection operator for $N$ particles,

$$P^{N}={\textstyle{\frac{1}{2\pi }}}\int d\phi \ e^{i\phi (\hat{N}-N)},$$ (68)

where $\hat{N}$ is the number operator, the average HFB energy of the particle-number projected state can be expressed as an integral over the gauge angle $\phi$ of the Hamiltonian matrix elements between states with different gauge angles [37,38]. In particular, for the Skyrme-HFB method implemented here, the particle-number projected energy can be written as [26,25]

$$\textsf{E}^{N}[\rho,\tilde{\rho}]=\frac{\left\langle \Phi \v... ...\int d\phi ~y(\phi )\int d^3{\bf r}~{\cal H}({\bf r},\phi )~,$$ (69)

where the gauge-angle dependent energy density ${\cal H}({\bf r},\phi )$ is derived from the unprojected energy density ${\cal H}({\bf r})$ (10) by simply substituting the particle and pairing local densities $\rho ({\bf r})$, $\tilde{\rho}({\bf r})$, $\tau ({\bf r})$, and ${\bf J}_{ij}({\bf r})$ by their gauge-angle dependent counterparts $\rho ({\bf r},\phi )$, $\tilde{\rho}({\bf r},\phi )$, $\tau ( {\bf r},\phi )$, and ${\bf J}_{ij}({\bf r},\phi )$, respectively. The latter densities are calculated from the gauge-angle dependent density matrices as

$$\begin{array}{c} \displaystyle{\rho ({\bf r}\sigma ,{\bf r^{... ...,\sigma^{\prime })\Phi _{\alpha}({\bf r},\sigma )}, \end{array}$$ (70)

where the gauge-angle dependent matrix elements read

$$\begin{array}{c} \displaystyle{\rho_{\alpha^{\prime }\alpha}... ...beta}(\phi )\tilde{\rho}_{\beta\alpha^{\prime }}} , \end{array}$$ (71)

and depend on the unprojected matrix elements (47) and on the gauge-angle dependent matrix

$$C(\phi ) =e^{2i\phi }\left[ 1+\rho (e^{2i\phi }-1)\right]^{-1}.$$ (72)

Function $y(\phi )$ appearing in Eq. (69) is defined as

$$y(\phi ) =\frac{x(\phi )}{\int d\phi^{\prime }\,x(\phi^{\pri... ...}}}\frac{e^{-i\phi N}\det (e^{i\phi}I)}{\sqrt{\det C(\phi )}},$$ (73)

where $I$ is the unit matrix.

Since the gauge-angle dependent matrices (70) and (71) are all diagonal in the same canonical basis that diagonalizes the unprojected density matrices (47), all calculations are very much simplified when they are performed in the canonical basis. In particular, in the canonical basis the matrices (71) read

$$\rho_{\mu}(\phi) =\displaystyle \frac{e^{2\imath \phi}v_{\mu... ...phi} u_{\mu}v_{\mu}}{u_{\mu}^{2}+e^{2\imath \phi}v_{\mu}^{2}},$$ (74)

while the function $x(\phi)$ can be calculated as

$$x(\phi) =\displaystyle{e^{-\imath N\phi}\prod_{\mu >0}\left( u_{\mu}^{2}+e^{2\imath \phi}v_{\mu}^{2}\right) },$$ (75)

where $v_{\mu}$ and $u_{\mu}$ ( $v_{\mu}^2+u_{\mu}^2=1$) are the usual canonical basis occupation amplitudes.

All the above expressions apply to independently restoring the proton and neutron numbers, so, in practice, integrations over two gauge angles have to be simultaneously implemented. In practice, these integrations are carried out by using a simple discretization method, which amounts to approximating the projection operator (68) by a double sum [39], i.e.,

$$P^{NZ}=\frac{1}{L}\sum_{l_{n}=0}^{L-1}\ e^{i\phi _{n}(\hat{N... ...}\frac{1}{L} \sum_{l_{p}=0}^{L-1}\ e^{i\phi _{p}(\hat{Z}-Z)},$$ (76)

where

$$\phi_{q}=\frac{\pi }{L}l_{q},~~~~q=n,p.$$ (77)

Usually no more than $L=9$ points are required for a precise particle number restoration.