The Numerov algorithm for the HFB equation

The Numerov (Cowell) method [23] is the standard technique for a precise integration of second-order differential equations; here we only give some details pertaining to its application to the system of two HFB equations (34). Let function $y(r)$ be the solution of the differential equation

$$y''=F(y,r)\,.$$ (59)

The Numerov algorithm is based on the finite difference formula for three consecutive points on a mesh with step $h$,

$$y_{n+1}-2 y_n+y_{n-1}= \frac{h^2}{12}\left(y''_{n-1}+10 y''_n+y''_{n+1}\right)\,+\, ({h^6}),$$ (60)

where o$({h^6})$ means that the error in this relation is of the order $h^6$. Combining this formula with (60), one obtains the relation

$$\left(1-\frac{h^2}{12}F_{n+1}\right)y_{n+1} -2\left(1+\frac{5h^2}{12}F_n\right)y_n +\left(1-\frac{h^2}{12}F_{n-1}\right)y_{n-1}=0 ,$$ (61)

which is obviously also valid if $y$ is a vector and $F$ a square matrix. Introducing notations

$$A_n=\left(1-\frac{h^2}{12}F_n\right),$$ (62)

and

$$G_n=A_n y_n\,,$$ (63)

one obtains the recurrence relation (Numerov iterations),

$$\left\{\begin{matrix}G_{n+1} =& 12 y_n - 10 G_n - G_{n-1}\,, \hfill\\ [1.5mm] \hfill y_{n+1} =& A_{n+1}^{-1}G_{n+1}\,. \hfill\\ \end{matrix}\right.$$ (64)

In the present applications, $A_n$ is a $2\times2$ matrix (position and energy dependent) and the calculation of its inverse could be the bottleneck of the procedure if it is not computed cleverly. This can be done in the following way. Using the notations introduced in Sec. 2.5, one can write matrix $A_n$ in an explicit form as

$$A_n=\left(\begin{matrix}1 - \frac{h^2}{12}\frac{V_n}{M^*_n}+\frac... ...^2}{12}\frac{V_n}{M^*_n}-\frac{h^2}{12}\frac{E}{M^*_n} \\ \end{matrix}\right) ,$$ (65)

whereby its inverse is given by

$$\begin{multline} \displaystyle A^{-1}_n=\frac{1}% {\left(1 - \frac{h^2 V_n}{12 M... ...}{M^*_n}+\frac{h^2}{12}\frac{E}{M^*_n} \\ \end{matrix}\right)\,. \end{multline}$$

It appears clearly that a lot of time can be saved by calculating and storing the energy independent terms,

$$1 - \frac{h^2 V_n}{12 M^*_n}\,, \ \frac{h^2}{12 M^*_n}\,, \ \frac{h^2 W_n}{12 M^*_n} \ \ \ \left(1 - \frac{h^2 V_n}{12 M^*_n}\right)^2+ \left(\frac{h^2 W_n}{12 M^*_n}\right)^2 ,$$ (66)

only once in each $(\ell j q)$-block; then the energy-dependent terms can be in each iteration calculated very rapidly.

The eigen energies are found by integrating two linearly independent solutions, $y^{(1+)}(r)$ and $y^{(2+)}(r)$, from the origin to a given point $r_m$, called the matching point. These two solutions are selected by imposing that either the first or the second component of the quasiparticle wave function vanishes near the origin, i.e.,

$$y^{(1+)}(r)\sim\left(\begin{matrix}r^{\ell+1}\\ 0\\ \end{matrix}\... ...(2+)}(r)\sim\left(\begin{matrix}0 \\ r^{\ell+1}\\ \end{matrix}\right)\,, \ \ \ \ r\rightarrow 0\,,$$ (67)

whereby we have

$$y^{(1+)}_0 = \left(\begin{matrix}0\\ 0\\ \end{matrix}\right)\,, \... ...mment_mark>181 y^{(2+)}_0 = \left(\begin{matrix}0 \\ 0\\ \end{matrix}\right)\,,$$ (68)

and

$$y^{(1+)}_1 = \left(\begin{matrix}Ah^{\ell+1}\\ 0\\ \end{matrix}\r... ...>182 y^{(2+)}_1 = \left(\begin{matrix}0 \\ Bh^{\ell+1}\\ \end{matrix}\right)\,.$$ (69)

Next, another two linearly independent solutions, $y^{(1-)}(r)$ and $y^{(2-)}(r)$, are found by the backward integration from the wall of the box $R_{\mathrm{box}}\equiv r_N$ to $r_m$, again imposing that either the first or the second component of the quasiparticle wave function vanishes near $r_N$ for the Dirichlet boundary conditions, i.e.,

$$y^{(1-)} \sim \left(\begin{matrix}r-r_N\\ 0\\ \end{matrix}\right)... ...85 y^{(2-)} \sim \left(\begin{matrix}0 \\ r-r_N\\ \end{matrix}\right)\,, \ \ \ \ r\rightarrow r_N\,,$$ (70)

whereby we have

$$y^{(1-)}_N = \left(\begin{matrix}0\\ 0\\ \end{matrix}\right)\,, \... ...mment_mark>186 y^{(2-)}_N = \left(\begin{matrix}0 \\ 0\\ \end{matrix}\right)\,,$$ (71)

and

$$y^{(1-)}_{N-1} = \left(\begin{matrix}Ch\\ 0\\ \end{matrix}\right)... ..._mark>187 y^{(2-)}_{N-1} = \left(\begin{matrix}0 \\ Dh\\ \end{matrix}\right)\,.$$ (72)

For the Neumann boundary conditions, at the left-hand-sides of Eqs. (71)-(73) one replaces functions by derivatives. Finally, one obtains four solutions valid everywhere except at the matching point $r_m$ and depending on the four constants $A$, $B$, $C$, and $D$.

The boundary conditions (69)-(70) and (72)-(73) determine the initial values $(G_0, G_1)$ and $(G_N, G_{N-1})$, respectively, needed to start the Numerov iterations of (65). Obviously, a correct limit should be taken when calculating values of $G_0$ from Eq. (64). It turns out that one obtains $G_0$=0 for all partial waves except $\ell=1$, which is the case that deserves a little more attention. For example, focusing our attention on one of the solutions,

$$y^{(1+)}(r)=\left(\begin{matrix}Ar^2\\ 0\\ \end{matrix}\right)\,,$$ (73)

we use the fact that in the vicinity of the origin $A(r)$ is dominated by the centrifugal term, so that its diagonal matrix elements are close to $-\frac{h^2}{12}\frac{2}{r^2}$. This leads to the result

$$G^{(1+)}_0=\lim_{r\rightarrow 0}A(r)y(r)= - \frac{Ah^2}{6} \left(\begin{matrix}1 \\ 0 \\ \end{matrix}\right).$$ (74)

An analogous result can be obtained for the other solution at the origin.

The eigenenergies are found by requiring that the full solution $y=A y^{(+)}_1 + B y^{(+)}_2 = C y^{(-)}_1 + D y^{(-)}_2$ and its first derivative are continuous at $r=r_m$ (see for example [24]). This matching conditions read

$$\left(\begin{matrix}y_m^{(1+)} & y_m^{(2+)} & - y_m^{(1-)} & - y_... ...ix}\right)=M_4 \left(\begin{matrix}A \\ B \\ C \\ D \\ \end{matrix}\right)=0\,.$$ (75)

The HFB quasiparticle wave functions having two components, the matching conditions amount to searching for zeros of the determinant of the $4\times4$ matrix $M_4$ in function of the quasiparticle energy. For $\ell=0$, the matching point $r_m$ is chosen to be near the point where the central part of the Hartree-Fock potential is half of its maximum, while for increasing values of $\ell$ it is gradually shifted to the exterior. A systematic search allows to find the intervals where the zeros of this determinant are located, and then the Newton method is used to find their exact locations to an arbitrary precision.

For vanishing determinant, Eq. (76) is used to determine three of the multiplicative constants in terms of the fourth one. This is done by extracting a $3\times3$ matrix $M_3$ out of $M_4$ and inverting it. Among the four possible choices one may have for $M_3$ we retain the one for which the product $M_3\times M_3^{-1}$ is closest to unity. Finally, the fourth multiplicative constant is fixed by the normalization condition of the wave function.

When approaching the self-consistent solution, the eigenenergies do not change very much between two consecutive Bogolyubov iterations. We take advantage of this fact by keeping in memory the intervals where the solutions are located and using them as initial guesses for the next iteration. This results in finding solutions for eigenenergies faster and faster as we get closer to the converged self-consistent state.