Density-dependent terms with fractional powers

Let us now analyze the terms in the DFT energy density that depend on fractional powers $\alpha$ of the local density. In many functionals related to the Skyrme interaction, and for the Gogny force, such terms are quite often postulated, both in the particle-hole and pairing channels (see Ref. [12] for a review). In particular, the familiar density-dependent term of the Skyrme force, which is proportional to $\rho^\gamma({\bf r})$, produces a contribution of the order of $\alpha=2+\gamma$ to the DFT energy density. Similarly, the density-dependent, zero-range term of the Gogny force yields a contribution to the DFT energy density that is of $\alpha=1+\gamma$ order, provided the particle-hole and pairing terms are consistently added. By taking into account the degeneracy factors $q$ discussed in Sec. 3.3, the resulting poles are of the order of $\alpha-p=2+\gamma-p$ and $\alpha-p=1+\gamma-p$, respectively. Since typical values of $\gamma$ are between 0 and 1, the DFT transition energy density always has poles at $z_n = \pm{i}\vert u_n/v_n\vert$ for non-degenerate states ($q=1$). But more importantly, the fractional powers lead to the multivalued DFT transition energy density on the complex-$z$ plane.

In the standard treatment of fractional powers $\alpha$ of a complex function, cuts along the negative real axis must be introduced. In order to apply this procedure to fractional powers of the local transition density (46), we must identify on the complex $z$ plane the lines along which $\rho_z({\bf r})$ is real negative.

Obviously, $\rho_z({\bf r})$ is real positive along the real $z$ axis and real along the imaginary $z$ axis. In order to simplify the discussion, let us assume that the sum in Eq. (46) is finite, which is always the case in any practical calculation. In such a case, $\rho_z({\bf r})$ has a finite number, say $M$, of different first-order poles along the positive imaginary axis, and the same number $M$ of poles along the negative imaginary axis. Moreover, since all coefficients in Eq. (46) are positive, $\rho_z({\bf r})$ must have a first-order zero between each pair of poles on the positive imaginary axis, and similarly on the negative imaginary axis. Since $\rho_z({\bf r})$ also has a second-order zero at $z=0$, we conclude that it has altogether $2M$ zeros on the imaginary axis. It is also obvious that $\rho_z({\bf r})$ is a rational function with a $2M$-order polynomial in the numerator, and thus we conclude that all the zeros of $\rho_z({\bf r})$ are located on the imaginary axis. Therefore, the cuts for possible fractional powers $\alpha$ must be located along the imaginary axis, and connect zeros of $\rho_z({\bf r})$ with its adjacent poles.

The above discussion is visualized in Fig. 2. The left portion shows schematically the transition density $\rho_{\mbox{\rm\scriptsize {Im}}[z]}({\bf r})$ along the imaginary axis $\mbox{\rm\scriptsize {Im}}[z]$ oriented vertically. The plot illustrates the transition density (46) in one selected point of space ${\bf r}$, i.e., values of wave functions at ${\bf r}$ enter only as numerical coefficients. There appear four poles and three zeros of $\rho_{\mbox{\rm\scriptsize {Im}}[z]}$ on the positive imaginary axis, the same number of poles and zeros on the negative imaginary axis, and the second-order zero at the origin. Sections of the imaginary axis where the density is negative are shaded. In the right portion of Fig. 2 we show poles (crossed circles) and zeros (full dots) of the transition density on the complex $z$ plane, along with the three integration contours $C0$, $C1'$, and $C2'$ discussed above. The cuts in the complex $z$-plane connecting zeros and poles, corresponding to real negative values of $\rho _z$, are indicated by vertical segments.

Figure 3: (color online) Ordinary ($z=1$) and transition ($z=i$) densities in $^{18}$O calculated in HFB+SLy4 as functions of $r$. The upper part of the Figure shows (in extended scale) the small region of radii where the transition density becomes negative.

While the location of the poles is independent of ${\bf r}$, the position of zeros of the transition density is $\bf r$-dependent. In order to visualize this, in Fig. 3 we plot the total ordinary ($z=1$) and transition ($z=i$) densities in $^{18}$O obtained within the SLy4 energy density functional. One can see that the transition density is positive almost everywhere; only in a very narrow region near $r\simeq5$fm does it become slightly negative, as shown in the upper part where the scale is expanded by a factor of 1000. Such a behavior of $\rho_z({\bf r})$ at $z^2=-1$ can be easily understood from Eq. (46). Indeed, strongly occupied states with $v_n^2\simeq 1$ always yield positive contributions, while negative contributions of states with $v_n^2<u_n^2$ can only appear in the surface region where the least bound canonical states dominate.

Figure 4: (color online) Bottom: ordinary ($z=1$) and transition ($z=i$) densities in $^{26}$O calculated in HFB+SLy4 as functions of $r$. Top: the real parts (solid lines) and imaginary parts (dashed lines) of the 1/6 powers of the ordinary ($z=1$) and transition ( $z=ie^{\pm {i}\varepsilon }$) densities.

By the same token, we can see that strong negative contributions may appear when the integration contour passes slightly below a pole located just above the unit radius, i.e., $v^2_n$ is slightly smaller then $u^2_n$. Such a situation is predicted in $^{26}$O, where the occupation probability of the canonical 2d$_{3/2}$ state equals 0.486. As shown in the bottom part of Fig. 4, the transition density at $z=i$ is dominated by this particular contribution and becomes strongly negative beyond $r\simeq2$fm.

We are now ready to discuss contour integration of terms depending on fractional powers of the transition density. Contour $C0$ shown in Fig. 2 crosses the imaginary axis in sections where there is no cut, and thus it always stays on the same Riemann sheet. On the other hand, contours $C1$ and $C2$ of Fig. 1, or contours $C1'$ and $C2'$ of Fig. 2, cross the imaginary axis by passing through cuts onto another Riemann sheet. Since the transition density (46) is an even function of $z$, the phase of the fractional power $\alpha$ of the transition density increases or decreases by $2\pi\alpha$ when going across each of the two cuts. Therefore, after returning to $z=1$, $\rho_z({\bf r})$ is multiplied by $\exp(\pm4\pi\alpha)$, and thus it is not a continuous function at $z=1$, unless $\alpha=k/2$. This is quite unacceptable as the presence of the phase creates serious problems in interpreting the projected DFT energies (34) (see, e.g., the sum rule condition discussed in Sec. 3.5).

Formally, by using powers of square roots in density-dependent terms, i.e., $\alpha=k/2$, one can guarantee that the integration contours return onto the original Riemann sheet, and that the transition energy density is a continuous function of $z$. However, even in such a case, one important property of the DFT transition energy density (33) is lost, namely, the density-dependent term in the energy density becomes an odd function of $z$, and the corresponding term in the integrand (35) becomes an even function of $z$. This is so, because the square root has opposite signs on the two Riemann sheets in question. Consequently, contour integrals of such terms would vanish and the density-dependent terms would yield zero contribution to the PN-projected energy. This is a rather disastrous result. Hence, we are forced to conclude that the use of continuous contours for fractional powers is not a viable prescription for constructing the projected DFT energies.

Let us now discuss the way of evaluating contour integrals in all practical PNP calculations up to now. Unfortunately, such calculations have always disregarded the analytic structure of the underlying integrands. In fact, the fractional powers of transition density,

$$\rho^\alpha_z({\bf r})= \vert\rho_z({\bf r})\vert \exp\left\{i\arg\left[\rho_z({\bf r})\right]/\alpha\right\} ,$$ (52)

are practically determined by computer compilers. In Eq. (52), the so-called argument $\arg\left[\rho_z({\bf r})\right]$ of $\rho_z({\bf r})$ is defined as the phase of the complex variable $\rho_z({\bf r})$; hence, it is contained in the interval of $-\pi\div\pi$. This usual prescription corresponds to stepping over the cut whenever the contour approaches the imaginary axis for $\arg\left[\rho_z({\bf r})\right]=\pm\pi$, i.e., for real negative transition densities $\rho_z({\bf r})<0$. In this way, the integrand is always calculated on the same Riemann sheet, but the integration contour is not closed.

Figure 5: Modified closed contour $C1''$ that encircles the zero of the transition density. The poles (zeros) of $\rho _z$ are marked by crossed circles) (dots). The equivalent circular contour $C1'''$ lying between the zero of $\rho _z$ and the previous pole is also indicated. See text for more details.

The contour can be closed by adding a piece that goes around the zero of the transition density $\rho_z({\bf r})$. This is illustrated in Fig. 5 that shows a modification of contour $C1'$ of Fig. 2 near the positive imaginary $z$-axis. [An analogous mirror-like detour is made near the negative imaginary axis.] The resulting contour $C1''$ always stays on the same Riemann sheet and, therefore, the integration result does not depend on the radius. On the other hand, the contribution due to the additional path surrounding the zero is affected by the discontinuity of the integrand along the cut. Such a discontinuity in the transition density is shown in the top panel of Fig. 4 for $\alpha$=1/6 (52). Since the ordinary density ($z=1$) is real and positive, its fractional power is also real and positive. On the other hand, transition densities (46) corresponding to $z_\pm=ie^{\pm{i}\varepsilon}$ with $\varepsilon=0.5^\circ$, i.e., near the positive imaginary axis, are complex. While their real parts are practically identical on both sides of the cut, their imaginary parts have opposite signs; hence, a discontinuity is encountered. [Of course, since the directions of integration are opposite, the contributions to the PNP energy from both segments $z_\pm$ of the additional path are identical.]

When transforming the contour integration in Eq. (34) into the integral along the imaginary axis $y$, one must do a change of variables from $z$ to $iy$. This introduces the additional factor $i^{-N}$ in the integrand. For that reason, for even $N$, the discontinuity in the imaginary part of the density-dependent term of fractional order contributes to the real part of the projected DFT energy. As the discussion in Sec. 3.4 proves, the same holds for odd values of $N$.

Figure 5 also shows the circular contour $C1'''$ lying between the zero of $\rho _z$ and the previous pole $\vert z_{n-1}\vert$. This contour is formally equivalent to the deformed contour $C1''$ but it is easier to handle in practical applications. The radius of $C1'''$ must be slightly greater than $\vert z_{n-1}\vert$ and smaller than the lowest zero of $\rho_z({\bf r})$, minimized over the whole space $\bf r$, associated with the branching point corresponding to $z_n$. The use of contour $C1'''$ guarantees that the integration of fractional-order terms is done properly.

Altogether, blind application of prescription (52) can lead to spurious and entirely uncontrolled contributions to the projected DFT energies. Excepting Ref. [23], this fact has been entirely overlooked in all practical applications of the PNP method to date, and casts serious doubts on the reliability of the obtained results. Largest contributions are, of course, obtained when the integration contour passes slightly below a pole of the DFT transition energy density. For the Skyrme functionals in Sec. 4.2, we present specific examples of such situations.

The appearance of spurious contributions is, in fact, independent of the order of divergence at the pole. Therefore, it also shows up for ``integrable" poles, diverging with powers of $\alpha-p=1+\gamma-p<1$, discussed for the Gogny force in Ref. [37].