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 be the solution of the differential equation

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

where

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

(62) |

and

one obtains the recurrence relation (Numerov iterations),

In the present applications, is a 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 in an explicit form as

(65) |

whereby its inverse is given by

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

and | (66) |

only once in each -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, and , from the origin to a given point , 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.,

whereby we have

and

Next, another two linearly independent solutions, and , are found by the backward integration from the wall of the box to , again imposing that either the first or the second component of the quasiparticle wave function vanishes near for the Dirichlet boundary conditions, i.e.,

whereby we have

and

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 and depending on the four constants , , , and .

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

(73) |

we use the fact that in the vicinity of the origin is dominated by the centrifugal term, so that its diagonal matrix elements are close to . This leads to the result

(74) |

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

The eigenenergies are found by requiring that the full solution and its first derivative are continuous at (see for example [24]). This matching conditions read

The HFB quasiparticle wave functions having two components, the matching conditions amount to searching for zeros of the determinant of the matrix in function of the quasiparticle energy. For , the matching point 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 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 matrix out of and inverting it. Among the four possible choices one may have for we retain the one for which the product 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.