Pairing correlations and divergence of the energy

The use of a zero range interaction in the particle-particle channel leads to the energy divergence problem [9,19] related to the fact that the short-range behaviour of the abnormal density is

$$\tilde\rho(r)=\lim_{r'\rightarrow r}\tilde\rho(r,r') \propto\frac{1}{\vert r-r'\vert} .$$ (37)

The other abnormal densities diverge as well, but along that the pairing correlations also induce the divergence of the kinetic density $\tau(r)$. One simple way to overcome the divergence of the energy is to introduce a cut-off energy $E_{\mathrm{cut}}$ in the summations of all the densities

$$\sideset{}{_i}\sum_{E_i>0} \ \rightarrow\ \sideset{}{_i}\sum_{0<E_i<E_{\mathrm{cut}}}\,.$$ (38)

This cut-off somehow simulates the finite range of the interaction in the pairing channel. It therefore constitutes an additional parameter of the interaction and must be adjusted along with the other quantities which define the pairing interaction, see e.g. Ref. [18].

A more elegant way to prevent the divergence was recently proposed in Refs. [9,19]. It is based on the subtractions of the divergent part of the abnormal density. In its present implementation in the code HFBRAD, this method can only be applied if the pairing force does not contains non-local terms (i.e. gradient terms). In that case, the pairing energy depends only on $\tilde\rho$ (and possibly on $\rho$ but not on $\tilde\tau$ or $\tilde J$), i.e.,

$$\mathcal{E}_{\mathrm{pair}}=\sum_q \left[\frac{t_0'}{4}(1-x_0')+\... ...{24}(1-x_3')\rho^{\gamma'}\right] \tilde\rho_q^2=\sum_q g[\rho]\tilde\rho_q^2 ,$$ (39)

and the pairing potential depends linearly on $\tilde\rho$,

$$\tilde U_q=2 g[\rho]\tilde\rho_q\,.$$ (40)

The regularized pairing field can then be written as

$$\tilde U^{\mathrm{reg}}_q=2 g[\rho]\tilde\rho_q^{\mathrm{reg}} \equiv2 g_{\mathrm{eff}}[\rho]\tilde\rho_q\,.$$ (41)

Calculation of the various densities is performed by summing up the quasiparticle contributions up to a maximum energy $E_{\mathrm{max}}$, but because of the regularization, this maximum energy has not the meaning of a cut-off energy $E_{\mathrm{cut}}$.

It can be shown [19] that when we evaluate the total energy of the system, the pairing energy must be defined as

$$\mathcal{E}_{\mathrm{pair}}=\sum_q g_{\mathrm{eff}}[\rho]\tilde\rho_q^2\,.$$ (42)

With this definition, the kinetic and pairing energies both diverge but the sum, and indeed the total energy, does not [19]. If the pairing energy density involves kinetic or gradient terms (like it was implemented for the Skyrme interaction SkP [5]), the regularization would require the subtraction of higher order terms. This possibility in not implemented in the present version of the code. Finally, let us note that relations (21)-(26) in Ref. [19] imply that the variation of quantities $k_F$ and $k_c$ with respect to the densities should be included in the mean fields. This contribution is supposed to be small and is not taken into account in the code. On the other hand, for density-dependent pairing interaction, the variation with respect to the explicit density dependence is taken into account in the code, and gives a contribution to the energy of the order of a few keV.