Coulomb Interaction

In the case of proton states, one has to add to the central potential the direct Coulomb field

$$V^C_d({\bf r})=e^2 \int d^3 {\bf r^\prime} \frac{\rho_p({\bf r^\prime})}{\vert{\bf r}-{\bf r^\prime}\vert},$$ (50)

as well as the exchange Coulomb field, which in the present implementation is treated within the Slater approximation:

$$V^C_{ex}({\bf r})= -e^2\left({\textstyle{\frac{3}{\pi}}}\right)^{1/3} \rho_p^{1/3}({\bf r}) .$$ (51)

The integrand in the direct term (50) has a logarithmic singularity at the point ${\bf r}$= ${\bf r^\prime}$. A way to bypass this difficulty is to use the Vautherin prescription [30], i.e., to employ the identity

$$\triangle_{{\bf r}^\prime}\vert{\bf r}-{\bf r^\prime}\vert= 2/\vert{\bf r}-{\bf r^\prime}\vert,$$ (52)

and then integrate by parts the integral in Eq. (50). As a results, one obtains a singularity-free expression

$$V^C_d({\bf r})=\frac{e^2}{2} \int d^3 {\bf r^\prime}~ \vert{... ...prime}\vert~\triangle_{{\bf r}^\prime} \rho_p({\bf r^\prime}).$$ (53)

In cylindrical coordinates, after integrating over the azimuthal angle $\varphi$, one finds

$$\displaystyle V^C_d(r^\prime,z^\prime)=2 e^2 \int_0^\infty r ... ...\left(\frac{4rr^\prime}{d(r,z)}\right) ~\triangle \rho_p(r,z),$$ (54)

where $d(r,z)=\left[(z-z^\prime)^2+(r+r^\prime)^2\right]$ and $E(x)$ is the complete elliptic integral of the second kind that can be approximated by a standard polynomial formula [31].

Equivalently, one can use the prescription developed originally for calculations with the finite-range (Gogny) force [3]. It consists of expressing the Coulomb force as a sum of Gaussians:

$$\displaystyle \frac{ 1 }{ \vert {\bf r} - {\bf r}' \vert } =... ... \mu }{ \mu^{2}}~ e^{-\frac{({\bf r} - {\bf r}')^2 }{ \mu^2}},$$ (55)

which gives

$$\displaystyle V^C_d({\bf r})=e^2 \frac{ 2 }{ \sqrt{\pi}} \in... ...rac{\mbox{\scriptsize {d}} \mu }{ \mu^{2}}~ I_{\mu} ({\bf r}),$$ (56)

where the integral

$$\displaystyle I_{\mu} ({\bf r}) = \int\mbox{\scriptsize {d}}^... ...{-\frac{({\bf r} - {\bf r}')^2 }{ \mu^2}} \; \rho ({\bf r}')$$ (57)

can be easily calculated in cylindrical coordinates. After integrating over the azimuthal angle $\varphi$, one finds

$$\displaystyle I_{\mu}(r^\prime,z^\prime)= 2 \pi \int_0^\infty... ...^2}} I_0\left(\frac{2 r r^\prime}{\mu^2}\right)~ \rho_p(r,z),$$ (58)

where $I_0(x)$ is the Bessel function that can also be approximated by a standard polynomial formula [31].

In order to perform the remaining one-dimensional integration in Eq. (56), the variable $\mu$ is changed to

$$\xi = b / \sqrt{b^2 + \mu^{2}} ,$$ (59)

where $b$ is the largest of the two HO lengths $b_z$ and $b_\bot$. This change of variable is very convenient, since then the range of integration becomes [0, 1]. The integral (56) is accurately computed by using a 30-point Gauss-Legendre quadrature with respect to $\xi$.

We have tested the precision of both prescriptions, Eqs. (53) and (56), and checked that the second one gives better results within the adopted numbers of Gauss-Hermite and Gauss-Laguerre points that are used for calculating proton densities. Therefore, in the code HFBTHO (v1.66p) this second prescription is used, while the first one remains in the code, but is inactive.