Deformation energy within the HFB+PNP method

Results in this section were obtained with the SIII Skyrme force whose density-dependent terms do not create additional problems (cf. Sec. 4). Let us first analyze the simpler case of $^{18}$O. Figure 10 presents the deformation energy $E(\beta )$ as a function of the quadrupole deformation $\beta$. As a reference curve, we show the unprojected deformation energy emerging from the HFB calculations and the associated PNP energy curve (PAV; solid squares; the contour radius $r_0$=1). Near the ground state ($\beta=0$), the projected energy is lowered by about 1.5MeV due to additional PN correlations. At larger deformations, the correlation energy decreases due to the stronger static neutron pairing.

Figure 10: (color online) Deformation energy $E(\beta )$ as a function of quadrupole deformation $\beta$, calculated for $^{18}$O with the SIII force and volume pairing interaction. Results of HFB (open triangles) are compared with different variants of PAV (squares and circles) and VAP PNP (stars) (see text for details).

At deformations $\beta\simeq-0.12$ and $\beta\simeq+0.12$, the projected energy curve exhibits unphysical jumps. By comparing with Fig. 7, one concludes that at these deformations the neutron 1d$_{5/2}$ poles cross the integration contour. Obviously, the residue contributions of these poles cause the sudden jumps in the deformation energy. The 1d$_{5/2,5/2}$ pole introduces a positive contribution at $\beta\simeq-0.12$, while the 1d$_{5/2,1/2}$ pole introduces another positive contribution at $\beta\simeq+0.12$. Based on this observation, two other sets of PAV calculations were carried out. The first calculation (open circles) was done by excluding contributions from the 1d$_{5/2}$ poles, as is the case for the ground state configuration. This can be accomplished by reducing the integration radius from $r_0=1$ to a smaller value of about $r_0=0.1$, cf. Fig. 7. At small deformations, $-0.12\leq\beta\leq0.12$, the new results are identical to those obtained with the unit circle, and at the larger deformations (prolate or oblate), the energy curve smoothly continues without any jump. Thus, in this example, an appropriate shift of the integration contour allows us to obtain smooth and unique projected energy. The second PAV curve (open squares) has been obtained by always including the lowest 1d$_{5/2}$ pole, i.e., by continuously varying the integration radius as a function of $\beta$, to ensure that it always stays between the first and second 1d$_{5/2}$ pole (cf. Fig. 7). The resulting energy curve coincides at large deformations with the standard PAV result, and then smoothly continues to $\beta=0$, where the 1d$_{5/2}$ poles ($q$=3) disappear.

Figure 10 also presents the fully self-consistent VAP results. Similar to the PAV case, two sets of calculations were performed. The solid (open) stars correspond to including (excluding) the contributions from the lowest 1d$_{5/2}$ poles. In both cases, one obtains smooth curves, which, beyond the spherical point, differ from one another.

In this rather simple case of $^{18}$O, both in the PAV and VAP calculations, one can avoid unphysical jumps of the projected energy curve by making a specific selection of ``active" poles that are considered during contour integration. Such selection of residues can, in principle, become a part of the definition of the projected energy. The variational principle can then be invoked to pick the selection that yields the lowest projected energy. In the discussed case of $^{18}$O, the PAV and VAP energies obtained by excluding the 1d$_{5/2}$ poles are the lowest, and they are smooth functions of deformation. Therefore, such a selection can be adopted for the final PNP energy in this nucleus. It is clear, however, that one cannot a priori tell which selection of poles leads to the lowest projected energy. For example, in heavier oxygen isotopes, the lowest energy is obtained by including some of the 1d$_{5/2}$ poles.

Let us now consider a more complicated case of the HFB+LN calculations for the neutron rich nucleus $^{32}$Mg. The total HFB energy (without the corrective $\lambda^{(2)}$ LN term) is shown in Fig. 11 as a function of $\beta$. Solid squares denote the result of PAV PNP calculations on the top of HFB+LN. At $\beta\simeq+0.1$, the PAV curve exhibits a small jump, after which its behavior changes character. This is clearly related to the proton 1d$_{5/2,5/2}$ pole crossing the unit circle, cf. Fig. 8. Otherwise, the PAV deformation energy is quite smooth as a function of deformation, despite the fact that three pairs of neutron poles cross the unit circle in the deformation range considered. This apparent lack of sensitivity to neutron poles can be traced back to the fact that they cross the integration contour at or near points where they pairwise cross one another. Therefore, the increasing degeneracy factor $q$ makes the poles disappear at the crossing points; hence, the total PAV curve behaves smoothly.

Figure 11: (color online) Deformation energy $E(\beta )$ as a function of quadrupole deformation $\beta$ calculated for $^{32}$Mg with the SIII force and volume pairing interaction. Results of the PAV HFB+LN calculations (squares and triangles) are compared with the VAP PNP results (dots). The standard HFB result is shown by open triangles.

Following the example of $^{18}$O, also for $^{32}$Mg we performed PAV calculations wherein we took into account contributions from all the poles that contribute to the ground state configuration ($\beta=0$). The resulting PAV curve (solid triangles in Fig. 11) behaves smoothly but shows unphysical behavior at very large deformations. An explanation of this artifact follows from Fig. 12, where we plot energy contributions from the most important poles: (i) the neutron 1d$_{3/2,3/2}$ pole (solid circles; it leaves the unit-circle at $\beta\simeq0.22$ and we have to add its contribution beyond this point); (ii) the neutron 1d$_{7/2,1/2}$ pole (solid squares; it enters the contour at $\beta\simeq0.22$ and we have to subtract its contribution beyond this point); and (iii) the proton 1d$_{5/2,5/2}$ pole (solid triangles; it leaves the contour at $\beta\simeq0.1$ and we have to add its contribution beyond this point). Interestingly, pole contributions to the total projected energy oscillate with deformation. As expected, the residues vanish when the corresponding poles cross at the integration circle, cf. Fig. 8. However, oscillations between the crossing points can become quite large, as is the case for the proton 1d$_{5/2,5/2}$ pole; hence, the projected energy curve obtained by keeping contributions from the ground-state set of poles acquires strong unphysical oscillations.

Figure 12: (color online) Contributions to the PAV energy of $^{32}$Mg from the three selected poles as a function of deformation $\beta$.

In a search for the most sensible method of calculating the projected energy curve within the PAV approach, we can employ a prescription whereby the number of the lowest poles is kept fixed within the integration contour. This can be realized by keeping $r_0$ between the poles or at the pole crossing point. Since at the crossing points the poles vanish, at least in SIII calculations, this results in a smooth energy curve. Such an option is shown in Fig. 11 with solid squares. The resulting curve is indeed very smooth; however, at $\beta\simeq0.4$ there appears an unphysical bump, which makes this option as unacceptable as the other one.

We have also attempted calculating the energy curve within the VAP approach. In principle, the VAP approach could have generated problems related to the fact that the density-dependent terms of fractional order may lead to large negative contributions (see Sec. 3.6). In practice, this is never the case because the VAP method [23] is not implemented through an explicit minimization of the energy, but is carried out by solving variational equations that have been derived with the same incorrect treatment of cuts in the complex plane. In this context, it is worth emphasizing that the appearance of poles never leads to infinite total energies, but to discontinuities in the total energy. Therefore, there is no danger that the minimization procedure may attract a solution towards a pole.

The main problem in implementing the VAP method is related to the fact that unprojected quantities, e.g., particle $\rho_{nn'}$ or pairing $\tilde{\rho}_{nn'}$ densities, lose their usual physical meaning [23] in VAP. They depend on the internal normalization $N'={\mathrm Tr} \rho$ of the density $\rho_{nn'}$ that is not related to the particle number $N$ onto which the state is projected. Neither the total VAP energy nor other projected observables depend on the normalization $N'$. However, depending on the choice of the internal normalization $N'$, one obtains different canonical occupation probabilities; hence, the associated poles $z_n$ are not distributed in the same way as in the unprojected HFB case. Depending on the internal normalization $N'$, different poles $z_n$ enter the integration contour, and the convergence procedure cannot be easily controlled.

Additional problems arise when two poles are nearly degenerate. Although at the point of degeneracy the poles disappear, when the distance between the poles is small but nonzero and the integration contour is between them, one faces significant instabilities of the constrained VAP problem. During the iteration of VAP equations, one or both poles enter or leave the integration contour. The poles create jumps in the projected energy and, which is even more important, they create jumps in the deformation. The numerical algorithm enters a `ping-pong' regime, which cannot be overcome, and one cannot converge to any solution. Figure 11 offers a good illustration of this problem. The converged VAP energies for $^{32}$Mg are shown with solid circles. The converged solution can be found only in limited regions of deformation. The 1d and 1f neutron poles close to the contour spoil the convergence in the regions of $\beta \approx -0.2$, 0.25, and $0.5$. The same is true for the proton 1d states around $\beta \approx -0.2$ and 0.1. As a result, the VAP procedure could be solved only within small deformation intervals around $\beta \approx -0.2$, 0, and 0.4.