For several decades, the ground state of nuclei have been studied within the self-consistent mean-field approximation. Such a description of the atomic nucleus, has ability to properly account for the bulk properties of nuclei such as masses, radii or shape. The Hartree-Fock method provides a good approximation of closed shell magic nuclei, however, the pairing correlations constitute an essential ingredient for the description of open shell nuclei. These effects are usually described by the Hartree-Fock plus BCS (HFBCS) or Hartree-Fock-Bogolyubov (HFB) methods.
Two main classes of numerical methods have been so far used for the implementation of the HFBCS or HFB methods. In the first one, one employs the expansion of the quasiparticle states on a discrete basis of orthogonal functions, usually provided by the harmonic oscillator potential. In this case the non-linear HFB(CS) equations are formulated in a matrix form and can be solved by using an iterative procedure (hereafter Bogolyubov iterations). In the second class, one uses the direct integration of the equations in the coordinate representation. This is usually done within the box boundary conditions to discretize the spectrum of quasiparticle states. The first approach is more elegant and efficient for finite-range interactions, however, it has a disadvantage that the use of harmonic oscillator wave functions may perturb the correct description of the asymptotics of the system. On the other hand, the asymptotic part of the quasiparticle wave functions does not suffer from this possible weakness when the problem is solved in the coordinate basis, but here one has to solve a set of integro-differential equations and as a consequence this method is usually only applied for the zero-range Skyrme force, for which the problem reduces to a set of differential equations.
The question of the asymptotic properties of the atomic nuclei (neutron skins, halos) is an urgent challenge in theoretical nuclear structure, and is particularly important due to the availability of beams of exotic nuclei1, which allow for relevant experiments. The HFB problem in coordinate basis, along with the use of the effective Skyrme force is indeed a powerful tool for studying ground-state properties of nuclei, and such studies have been and will be pursued especially in the regions of nuclear chart where the experiments are being performed.
Two computer codes that solve self-consistent equations on the basis have recently been published [1,2] for non-spherical shapes. We refer the reader to these publications for a review of approaches and methods that can be and have been used in such cases. Within the assumed spherical symmetry, which is the subject of the present study, the first self-consistent solutions were obtained in the 1970's [3,4], and these old codes were later used for decades by different groups, however, they have never been published.
The first description and implementation of the HFB problem within the coordinate-space spherical symmetry was presented in Ref. . The code written during this work was later updated many times since its early version, and distributed widely, however, neither this code was ever the object of an official publication. Indeed, the absence of accompanying manual makes it hard to use, and its older versions are sometimes used despite the fact that the newer ones have been made available. For coordinate-space spherical symmetry, the code solving the HFBCS equations for the Skyrme interaction was published in Ref.  and that for the relativistic mean field Hartree-Bogolyubov method in Ref. . A similar unpublished code also exists to treat pairing correlations for the finite-range Gogny interaction , although the code that would solve the full HFB problem within these conditions is not yet available.
The aim of this work is to provide a modern code (named HFBRAD) for the spherical Skyrme-HFB problem in the coordinate representation. The code constitutes a completely rewritten version of that constructed in Ref. , but it also features implementation of the pairing renormalization . Although such code has to be used with caution in open-shell nuclei, because there the deformation effects not included here are essential, it nevertheless provides a fast and easy tool for a quick estimate of nuclear properties, which can precede much more time consuming calculations beyond the spherical symmetry. In this paper we also provide a clear description of the different Skyrme forces implemented in the code, terms neglected in the energy functional, and parameters in the pairing channel.
This paper is organized as follows. In section 2 we briefly show how the Skyrme-HFB equations are derived and discuss their main properties. Notably, a significant part is devoted to the problem of the divergence of the energy due to the zero range of the force in the pairing channel. In section 3 we present several parametrizations of the Skyrme force implemented in the code, with a particular attention paid to the pairing channel. In section 4 we give definitions and meanings of the various quantities and observables which are the result of calculations. Some aspects of the numerical treatment of the problem are presented in section 5, and finally, sections 6 and 7 give a precise description of the input and output data files.