Coordinate-space solver for superfluid many-fermion systems with the shifted conjugate-orthogonal conjugate-gradient method

Self-consistent approaches to superfluid many-fermion systems in three dimensions (and their subsequent use in time-dependent studies) require a large number of diagonalizations of very large dimension Hermitian matrices, which results in enormous computational costs. We present an approach based on the shifted conjugate-orthogonal conjugate-gradient (COCG) Krylov method for the evaluation of the Green's function, from which we subsequently extract various densities (particle number, spin, current, kinetic energy, anomalous, etc.) of a nuclear system. The approach eschews the determination of the quasiparticle wave functions and their corresponding quasiparticle energies, which never explicitly appear in the construction of a single-particle Hamiltonian or needed for the calculation of various static nuclear properties, which depend only on densities. As benchmarks we present calculations for nuclei with axial symmetry, including the ground state of spherical (magic or semimagic) and axially deformed nuclei, the saddle point in the ${}^{240}\mathrm{Pu}$ constrained fission path, and a vortex in the neutron star crust, and demonstrate the superior efficiency of the shifted COCG Krylov method over traditional approaches.

Figure 1

The ellipse contour $C:z\left(\phi \right)$ in Eq. (50) adopted to perform the contour integral. The thick solid line represents the positive energy continuum quasiparticle states with $E>-\mu$, while the crosses represent the positive energy discrete bound quasiparticle states $0<E_i<-\mu$. Density of integration points is larger for parts of the contour close to the real axis (depicted by blue and red colors) and decreases as we go far away from the real axis (depicted by green color).

Figure 2

The convergence behavior of the shifted COCG method for the HFB equation with W-S potential. The top and bottom rows show results for spherically symmetric and deformed W-S potentials respectively. The spherically symmetric problem was solved on lattice ${16}^3$, while the deformed problem used a lattice ${16}^2\times 20$. In both cases, the lattice spacing is $dx=1.25\mathrm{fm}$. Panels (a) and (d) show the distribution of the norms of the residual vectors $\parallel r_n^{\sigma }\parallel$ on the semiellipse contour $z\left(\phi \right)$ for iteration number $n=4000$. Three colors (red, green, blue) correspond to different parts of the contour, as depicted on Fig. 1. Panels (b) and (e) show the convergence behavior of $\parallel r_n^{\sigma }\parallel$ for three points $z\left(0\right)$ (red), $z(\pi /2)$ (green), and $z\left(\pi \right)$ (blue) as a function of iteration number. Panels (c) and (f) show the difference between the value obtained from COCG method in the $n\mathrm{th}$ iteration and the exact value for the normal (red) and anomalous (blue) densities.

Figure 3

Comparisons were made between the local densities calculated by shifted COCG iteration and the direct diagonalization approaches along $z$ axis for fixed $x=0$ and $y=0$ coordinates with three different lattice constants $dx$. In the case with spherical symmetric potential (top row) the system size is $L_x=L_y=L_z=20$ fm, and in the deformed case (bottom row) the system size is $L_x=L_y=20$ and $L_z=25$ fm. The left column shows differences for the normal density $\delta \rho =|{\rho }_{\mathrm{cocg}}-{\rho }_{\mathrm{diag}.}|$, and the right column shows them for the anomalous density $\delta \nu =|{\nu }_{\mathrm{cocg}}-{\nu }_{\mathrm{diag}.}|$.

Figure 4

Total density distribution $\rho \left(\mathit{r}\right)={\rho }_n\left(\mathit{r}\right)+{\rho }_p\left(\mathit{r}\right)$ for the saddle-point configuration of ${}^{240}\mathrm{Pu}$ in the $y=0$ plane.

Figure 5

The lowest energy states generated by shifted COCG method for neutron background density $n=0.014{\mathrm{fm}}^{-3}$ [panels (a) and (b)] and $0.031{\mathrm{fm}}^{-3}$ [panels (c) and (d)]. In each box of size $75\times 60{\mathrm{fm}}^2$ we show the absolute value of the neutron pairing potential $\mathrm{\Delta }\left(\mathit{r}\right)$ (upper half) and the total density distribution $\rho \left(\mathit{r}\right)$ (lower half). The black lines separating blue and red regions correspond to a value $1\mathrm{MeV}$ for the paring and of $0.07{\mathrm{fm}}^{-3}$ for the density respectively.