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

Shi Jin<EMAIL_ADDRESS>Department of Physics,University of Washington,
Seattle, Washington 98195–1560, USA  Aurel Bulgac<EMAIL_ADDRESS>Department of
Physics,University of Washington, Seattle, Washington 98195–1560, USA  Kenneth
Roche<EMAIL_ADDRESS>Department of Physics,University of Washington,
Seattle, Washington 98195–1560, USA  Pacific Northwest National Laboratory,
Richland, Washington 99352, USA  Gabriel Wlazłowski<EMAIL_ADDRESS>Faculty of Physics, Warsaw University of Technology, Ulica Koszykowa 75,
00–662 Warsaw, Poland  Department of Physics,University of Washington,
Seattle, Washington 98195–1560, USA

###### Abstract

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 wavefunctions 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 semi-magic) and axially deformed nuclei, the saddle-point
in the 240Pu 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.

###### pacs:

21.60.-n, 21.60.Jz

††preprint: NT@UW-16-09

##  I introduction

Density functional theory (DFT) and other self-consistent approaches [Hartree-
Fock (HF), Hartree-Fock-Bogoliubov (HFB) or Hartree-de Gennes] have played an
essential role in studying the properties of most nuclei across the nuclear
chart Bender  _et al._ (2003a). Present phenomenological nuclear energy
density functionals (EDF) allow for accurate descriptions of many bulk
properties of nuclei such as masses, radii and shapes, transition matrix
elements, potential energy surfaces and related inertial parameters, and even
non-equilibrium properties, when extended to time-dependent phenomena. The
time-dependent extension is straightforward Runge and Gross (1984); Marques
_et al._ (2006, 2012) and provides a unified approach for the study of both
structure and reaction dynamics. While pairing correlations are absent in
closed shell magic nuclei, there are a lot of nuclear problems where accurate
description of the pairing correlations is crucial. Examples of such problems
include large amplitude collective motion of open shell nuclei and dynamics of
vortices in neutron star crust Ring and Schuck (2004); Bertsch (1980);
Barranco  _et al._ (1990); Bertsch and Bulgac (1997); Bulgac  _et al._ (2016);
Anderson and Itoh (1975); Wlazłowski  _et al._ (2016).

The evaluation of the nuclear DFT is numerically demanding, particularly if
one considers large fermionic systems in three dimensions and large
deformations, without any symmetry constraints. Over the years many iterative
approaches for solving the HFB equations have been proposed, including
successive diagonalizations of the HFB or HF+BCS matrices, imaginary time
evolution Davies  _et al._ (1980); Ryssens  _et al._ (2015a), and gradient
methods Egido  _et al._ (1995); Robledo and Bertsch (2011), which typically
need significant memory requirements and operations of complexity $O(N^{3})$,
where $N$ is the dimension of the HFB matrix. For a review of modern
diagonalization software, see Ref. Golub and van der Vorst (2000).

The most standard approach, via series of direct diagonalizations of the HFB
Hamiltonian, can be divided into two main classes. In the first one, the HFB
problem is formulated in the configuration space by expanding the
quasiparticle states of HFB on a discrete basis of orthogonal functions,
usually provided by a (deformed) harmonic oscillator (HO) basis Stoitsov  _et
al._ (2013); Schunck  _et al._ (2012). Although typically very fast, this
approach suffers from truncation errors that typically lead to the incorrect
description of the asymptotic behavior of the system or of the large
deformations of the systems, e.g. when a nucleus fissions. An HO-basis does
not usually provide a very effective coverage of the relevant phase-space.
Sometimes this can be improved by introducing wavelets Pei  _et al._ (2014)
and complex energy Berggren states Michel  _et al._ (2002) to describe the
continuum spectrum. Nuclear systems are adequately described when the phase-
space in which the dynamics occurs is properly covered by the single-particle
basis wave functions. This space is characterized typically by a volume
$V=L^{3}$, where $L$ is a few times the nucleus radius and by a maximum
single-particle momentum $p_{\mathrm{cut}}=\hbar k_{\mathrm{cut}}$
proportional to the Fermi momentum $\hbar k_{F}$. The total number of quantum
states in such p hase space is

$\displaystyle{\mathcal{N}}_{PS}=4\frac{(2p_{\mathrm{cut}}L)^{3}}{(2\pi\hbar)^{3}}\propto
k_{F}^{3}R^{3},$ (1)

where the factor 4 arises from accounting for spin and isospin degrees of
freedom, see the discussion in Ref. Bulgac and Forbes (2013). The spatial
extension $L$ (that can be different in each Cartesian direction) is chosen
depending on the specific needs. These extensions in each spatial direction
are different for a nucleus with very extended density tails, for the
collision of two nuclei, for the fission of a heavy nucleus and its split into
two fragments or for a nucleus in the presence or absence of a vortex in the
neutron star crust. In an HO-basis, by increasing the size of the basis set to
cover the required spatial volume, one usually goes well beyond
$p_{\mathrm{cut}}$, leading to an inefficient coverage of the needed phase-
space. Pairing correlations typically lead to a filling of all the momentum
states up to $p_{\mathrm{cut}}$. The need to describe large deformations, the
tails of the density distributions, and particularly the large amplitude
dynamics of various nuclear systems requires large spatial simulation volumes
and a large momentum cutoff, thus resulting in a large number of phase-space
cells ${\mathcal{N}}_{PS}$.

Another approach is the direct HFB matrix diagonalization in the coordinate
space with a lattice spacing $dx=\pi/k_{\mathrm{cut}}$ chosen to ensure an
adequate coverage of the phase-space Bulgac and Forbes (2013); Ryssens  _et
al._ (2015a, b). Thus, one can obtain numerically accurate results for weakly
bound nuclei and large deformations Dobaczewski  _et al._ (1984, 1996). A
number of coordinate-space HFB solvers have been published over the years
Zhang  _et al._ (2013); Bennaceur and Dobaczewski (2005); Oberacker  _et al._
(2003); Terán  _et al._ (2003); Blazkiewicz  _et al._ (2005); Pei  _et al._
(2008, 2014), but solving the HFB equations in full three-dimensional (3D)
coordinate space is still a challenging problem because of the large dimension
of the HFB matrix discretized in a large box. To put this in perspective
consider calculations in medium size volume $32\times 32\times
48\;\textrm{fm}^{3}$ with lattice resolution $dx=1\;\textrm{fm}$ corresponding
to $p_{\mathrm{cut}}=\hbar\pi/dx\approx 600\;\textrm{MeV}/c$. Then, the
$N\times N$ HFB matrix has $N^{2}=(4\times 32\times 32\times 48)^{2}\approx
200~{}000^{2}$ matrix elements, and requires more than $0.5\,\text{TB}$ memory
just to store it. A typical diagonalization (which requires $O(N^{3})$
operations separately for protons and neutrons), takes about 40 minutes using
the high-performance linear algebra library ScaLAPACK sca  on the Edison
supercomputer at NERSC with 36 864 processor cores and a charged computational
cost of 49 152 CPU hours edi  (CPU, central processing unit). For self-
consistent convergence $\sim 100$ diagonalizations are typically required,
generating an enormous computational cost on the order of about 10 million CPU
hours per converged calculation. After each diagonalization, the
eigenfunctions are reduced to local densities needed to construct new
quasiparticle Hamiltonian matrix coefficients before the start of a next
iteration. The manner in which new matrix coefficients are constructed
reflects the fact that, for density functional theory, the many-body wave
function contains vastly more information than is needed in each iteration.

Here we present a new method for extracting densities directly from the HFB
Hamiltonian without calculating wavefunctions. The method is especially well
suited for large scale calculations that inevitably require an efficient use
of supercomputers. It is important to compare both the computational
(numerical) complexity of different algorithms that solve the same problem,
and the ease with which the methods can be partitioned into smaller problems
that can be effectively executed in parallel so that the entire process scales
well on today’s computers. An important indicator that characterizes this
property of an given algorithm is the strong scaling, which describes how the
time to compute a fixed problem depends on the aggregated scale of the
computing resources used to finish the problem. Ideal (linear) strong scaling
for an algorithm is achieved when the (wall) time to completion can be reduced
by a constant factor $k$ while increasing the aggregated machine scale by the
same factor $k$. The dense, direct eigenvalue decomposition based on data
decomposing the Hamiltonian over a set of processes does not exhibit perfect
strong scaling on modern parallel computers, particularly for very large
matrices, when communication between computational nodes starts dominating the
computational cost. The factorization requires frequent interruptions both in
communicating partial results and coordinating coarse phases of the algorithm
between computing units, and these synchronizations, coupled to increasingly
smaller work fractions for a fixed problem dominate the strong scaling
behavior. Eventually more time is spent managing the computation than
evaluating the algorithm. Moreover, a typical local energy density functional
leads to a sparse HFB matrix in the coordinate representation - a feature that
is not efficiently utilized by the eigenvalue decomposition with direct
methods. The method proposed here removes these weaknesses. The algorithm can
be constructed in a nearly communication-free manner and thus exhibits almost
prefect strong scaling over the number of points in the coordinate space
offering a near complete reduction of $O(N)$ in complexity when evaluating the
method in parallel. The only operation that involves the HFB matrix is a
matrix-vector (MV) multiplication, which can benefit easily from matrix
sparsity. In our implementation the MV multiplication is extremely efficiently
implemented using fast Fourier transform. The programming method is
straightforward in hybrid processing models that combine traditional CPUs with
hardware accelerators such as general purpose graphics processing units (GPUs)
cud . Presently, many leadership class computers are of this type.

In order to grasp the new method, consider a Hartree-Fock equation Ring and
Schuck (2004)

$\displaystyle
H\psi_{k}(\boldsymbol{r})=\varepsilon_{k}\psi_{k}(\boldsymbol{r}),$ (2)

where $\psi_{k}(\boldsymbol{r})$ is the wavefunction corresponding to the
energy level $\varepsilon_{k}$. Without loss of generality, we assume that the
eigenvalues $\varepsilon_{k}$ are positive. Our goal is to compute the
particle number density

$\rho(\boldsymbol{r})=\sum_{\varepsilon_{k}<\varepsilon_{F}}\psi_{k}(\boldsymbol{r})\psi_{k}^{*}(\boldsymbol{r}),$
(3)

where the summation includes only states up to a fixed Fermi energy
$\varepsilon_{F}$ for normal systems. The density can be obtained from the
Green’s function $G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})$, defined by the
linear equation (spin degrees freedom are suppressed here)

$\displaystyle(z-H)G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})=\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime}),$
(4)

where $z$ is a complex number. Notice that in this equation
$\boldsymbol{r}^{\prime}$ can be treated as parameter and thus for a fixed
value of $\boldsymbol{r}^{\prime}$ the Green’s function
$G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})$ is the solution of an
inhomogeneous Schrödinger equation. The formal solution of this equation is

$\displaystyle
G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})=\sum_{k}\frac{\psi_{k}(\boldsymbol{r})\psi_{k}^{*}(\boldsymbol{r}^{\prime})}{z-\varepsilon_{k}}.$
(5)

Once we obtain the Green’s function, forcing
$\boldsymbol{r}=\boldsymbol{r}^{\prime}$, the particle density can be
calculated via a contour integral

$\displaystyle\rho(\boldsymbol{r}^{\prime})=\frac{1}{2\pi
i}\oint_{C}dz\,\left.G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\right|_{\boldsymbol{r}=\boldsymbol{r}^{\prime}}.$
(6)

A contour can be chosen arbitrarily, but is required to enclose only the
eigenvalues in the the interval $[0,\varepsilon_{F}]$ on the real axis. The
problem of computing the density for a fixed point $\boldsymbol{r}^{\prime}$
is reduced to the calculation of the Green’s function
$G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})$ for all $z$’s on the contour,
which is equivalent to solving a set of independent linear equations (4). For
any given local Hamiltonian the density at any given point $\boldsymbol{r}$
can be extracted independently of any other point, enabling this step to be
executed in parallel, with no communications between processes. Thus by
construction, the present method exhibits perfect strong scaling up to the
number of coordinate points.

For fixed $z$ and $\boldsymbol{r}^{\prime}$ the algebraic problem (4) can be
solved using iterative methods. The conjugate orthogonal conjugate gradient
(COCG) method is a very efficient algorithm for solving a set of linear
equations $Ax=b$, assuming that matrix $A$ is symmetric and complex van der
Vorst and Melissen (1990). Recently, an extension of the method called the
shifted COCG has been implemented for electron systems Yamamoto  _et al._
(2008); Takayama  _et al._ (2006), nuclear shell-model Mizusaki  _et al._
(2010), the computation of the level density in nuclei Shimizu  _et al._
(2016), and other generalized eigenvalue problems Sakurai and Sugiura (2003);
Sakurai  _et al._ (2008); Ikegami  _et al._ (2010). The shifted variant solves
simultaneously a family of algebraic problems $(A-\sigma)x^{\sigma}=b$ for
many shifts $\sigma\in\mathbb{C}$ simultaneously, essentially with the same
speed as standard COCG for a single shift. In this way, the accuracy of
numerical estimation of the contour integral (6) can be refined to the desired
accuracy with almost no extra calculation cost. Taking as many computing units
as points in coordinate space, the computation time will be limited only by
the time needed to solve this single algebraic problem.

The purpose of the present work is to introduce the procedure in the context
of solving the HFB equation in 3D coordinate space. The Green’s function
$G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})$ of a HFB equation can be
obtained in a similar way as solving Eq. (4) using the shifted COCG method,
but for a generalized multi-component system with pairing and spin-orbit
coupling. The calculation of particle density in Eq. (6) will be generalized
to all types of local densities in the realistic nuclear energy density
functional (NEDF) J. Dobaczewski  _et al._ (2002); Chabanat  _et al._ (1998).
With the aim of taking advantage of existing and future computational
resources, we developed a highly efficient parallelized GPU code as the so-
called engine of the shifted COCG iteration to replace the direct parallel
diagonalization procedure in the code used in Refs. Stetcu  _et al._ (2011,
2015); Bulgac  _et al._ (2016). As a benchmark, full self-consistent HFB
calculations are performed in this work for the ground states of spherical
(magic or semi-magic) as well as axially deformed nuclei. The constrained HFB
(CHFB) calculation is also tested for the saddle point of
${}^{240}\mathrm{Pu}$ in the induced fission studied in Bulgac  _et al._
(2016). Finally, we apply the method to generate nontrivial states, relevant
for astrophysical applications, containing quantized vortices and nuclear
defects immersed in a superfluid neutron matter.

We emphasize that the concept of DFT is extensively used across many fields
dealing with fermionic systems, like quantum chemistry, solid state physics,
ultracold fermionic gases and many others. Here we focus only on the nuclear
case, as typically the nuclear EDF is very complicated in comparison to
functionals encountered in other fields. The method presented here is general
and can be applied to other fermionic systems as well.

##  II Theoretical Framework

###  II.1 Nuclear DFT within local density approximation

Local density approximation is one of the most successful concepts introduced
to DFT approaches. It was extended to superfluid systems, namely superfluid
local density approximation (SLDA) Bulgac and Yu (2002); Yu and Bulgac (2003),
and assumes that energy density functional (EDF) $\mathcal{E}$ depends on
various local densities. In nuclear systems, a generic EDF is represented as a
sum of the kinetic $\mathcal{E}_{\mathrm{kin}}$, the nuclear
$\mathcal{E}_{\mathrm{nuclear}}$, the Coulomb $\mathcal{E}_{\mathrm{Coul}}$,
and the pairing $\mathcal{E}_{\mathrm{pair}}$ contributions

$\displaystyle\mathcal{E}=\mathcal{E}_{\mathrm{kin}}+\mathcal{E}_{\mathrm{nuclear}}+\mathcal{E}_{\mathrm{Coul}}+\mathcal{E}_{\mathrm{pair}}.$
(7)

The kinetic part depends on the kinetic densities and the effective nucleon
masses

$\displaystyle\mathcal{E}_{\mathrm{kin}}(\boldsymbol{r})=\sum_{q=n,p}\frac{\hbar^{2}\tau_{q}(\boldsymbol{r})}{2m_{q}(\boldsymbol{r})}.$
(8)

In calculations we included the simplest center of mass correction by
replacing the bare nucleon mass $m$ with $m/(1-1/A)$. The total particle
number is $A=N+Z$ where $N$ and $Z$ are respectively the neutron ($n$) and
proton ($p$) numbers. The Coulomb contribution is composed of the direct part
and the exchange part

$\displaystyle\begin{split}&\mathcal{E}_{\mathrm{Coul}}(\boldsymbol{r})=\mathcal{E}^{d}_{\mathrm{Coul}}(\boldsymbol{r})+\mathcal{E}^{e}_{\mathrm{Coul}}(\boldsymbol{r})\\\
&=\frac{e^{2}}{2}\int\frac{\rho_{p}(\boldsymbol{r})\rho_{p}(\boldsymbol{r}^{\prime})}{\lvert\boldsymbol{r}-\boldsymbol{r}^{\prime}\rvert}d^{3}\boldsymbol{r}^{\prime}-\frac{3e^{2}}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho_{p}^{4/3}(\boldsymbol{r}).\end{split}$
(9)

The pairing energy in SLDA depends on the local anomalous density:

$\displaystyle\mathcal{E}_{\mathrm{pair}}(\boldsymbol{r})=\sum_{q=n,p}g_{\mathrm{eff}}(\boldsymbol{r})\lvert\nu_{q}(\boldsymbol{r})\rvert^{2}$
(10)

and the effective pairing coupling strength $g_{\mathrm{eff}}(\boldsymbol{r})$
is obtained via a renormalization Bulgac and Yu (2002); Yu and Bulgac (2003);
Borycki  _et al._ (2006) of the bare pairing strength, typically parametrized
as

$\displaystyle
g_{0}(\boldsymbol{r})=g_{0}\left(1-\alpha\frac{\rho(\boldsymbol{r})}{\rho_{0}}\right),$
(11)

where $\rho_{0}=0.16~{}\mathrm{fm}^{-3}$ is the saturation density. The
parameter $\alpha=0,1,1/2$ corresponds to volume, surface, and mixed pairing
respectively J. Dobaczewski  _et al._ (2002); Bertsch  _et al._ (2009).

The nuclear part is the most complicated. Over the years, many forms of the
nuclear functional have been proposed, see Refs. Bender  _et al._ (2003a);
Bulgac  _et al._ (2015); Dutra  _et al._ (2012) for review. Typically, the
functional depends on various proton and neutron densities, including normal
$\rho(\boldsymbol{r})$, kinetic $\tau(\boldsymbol{r})$, spin
$\boldsymbol{s}(\boldsymbol{r})$ or spin kinetic energy densities
$\boldsymbol{T}(\boldsymbol{r})$. In a high accuracy nuclear EDF (NEDF)
various currents are present as well, such as the normal current
$\boldsymbol{j}(\boldsymbol{r})$, and the spin current densities
$\boldsymbol{J}(\boldsymbol{r})$. In our previous works Stetcu  _et al._
(2011, 2015); Bulgac  _et al._ (2016) we used the popular parametrization SLy4
Chabanat  _et al._ (1998); Bender  _et al._ (2003b) of the Skyrme NEDF, that
has a rather generic form,

$\displaystyle\begin{split}\mathcal{E}_{\mathrm{Skyrme}}&=\mathcal{E}_{\rho^{2}}+\mathcal{E}_{\rho^{\gamma}}+\mathcal{E}_{\rho\Delta\rho}+\mathcal{E}_{\rho\tau}+\mathcal{E}_{\rho\nabla
J}\\\
&=\sum_{t=0,1}\left(C_{t}^{\rho}\rho_{t}^{2}+C_{t}^{\gamma}\rho_{t}^{2}\rho_{0}^{\gamma}+C_{t}^{\rho\Delta\rho}\rho_{t}\Delta\rho_{t}\right.\\\
&\quad\quad\quad\left.+C_{t}^{\tau}\rho_{t}\tau_{t}+C_{t}^{\nabla
J}\rho_{t}\boldsymbol{\nabla}\cdot\boldsymbol{J}_{t}\right),\end{split}$ (12)

where $\rho_{0}=\rho_{n}+\rho_{p}$ and $\rho_{1}=\rho_{n}-\rho_{p}$ (and
similar for $\tau_{0,1}$ and $\boldsymbol{J}_{0,1}$) are isoscalar and
isovector densities respectively, and $C$’s are coupling constants.

Starting from the NEDF defined above, the derived HFB equation is a
4-component eigenvalue problem:

$\displaystyle\begin{gathered}H\begin{pmatrix}u_{k\uparrow}\\\
u_{k\downarrow}\\\ v_{k\uparrow}\\\
v_{k\downarrow}\end{pmatrix}=E_{k}\begin{pmatrix}u_{k\uparrow}\\\
u_{k\downarrow}\\\ v_{k\uparrow}\\\ v_{k\downarrow}\end{pmatrix}\\\
H=\begin{pmatrix}h_{\uparrow\uparrow}-\mu&h_{\uparrow\downarrow}&0&\Delta\\\
h_{\downarrow\uparrow}&h_{\downarrow\downarrow}-\mu&-\Delta&0\\\
0&-\Delta^{*}&-h^{*}_{\uparrow\uparrow}+\mu&-h^{*}_{\uparrow\downarrow}\\\
\Delta^{*}&0&-h^{*}_{\downarrow\uparrow}&-h^{*}_{\downarrow\downarrow}+\mu\end{pmatrix}\end{gathered}$
(15)

where we have suppressed the spatial coordinate $\boldsymbol{r}$ and $k$ is
the label of each quasiparticle wavefunction
$[u_{k\sigma}(\boldsymbol{r}),v_{k\sigma}(\boldsymbol{r})]$ where
$\sigma=\uparrow,\downarrow$. The local particle-hole Hamiltonian $h$ is
obtained by taking the appropriate functional derivatives of the energy
density functional. For the Skyrme functional (12) it takes the form Ring and
Schuck (2004):

$\displaystyle\begin{split}h_{\sigma,\sigma^{\prime}}(\boldsymbol{r})=\left(-\boldsymbol{\nabla}\cdot\frac{\hbar^{2}}{2m^{*}}\boldsymbol{\nabla}+U\right)\delta_{\sigma,\sigma^{\prime}}\\\
-i\boldsymbol{W}\cdot(\boldsymbol{\nabla}\times\boldsymbol{\sigma})_{\sigma,\sigma^{\prime}},\end{split}$
(16)

where $m^{*}(\boldsymbol{r})$ is the effective mass, $U(\boldsymbol{r})$ is
the central-part of the mean-field potential, and
$\boldsymbol{W}(\boldsymbol{r})$ is the spin-orbit potential (for their
explicit forms see Ref. Chabanat  _et al._ (1998)). The local pairing field
$\Delta(\boldsymbol{r})$ is defined as a function of the anomalous density

$\displaystyle\Delta(\boldsymbol{r})=-g_{\mathrm{eff}}(\boldsymbol{r})\nu(\boldsymbol{r}).$
(17)

The HFB Hamiltonian is a functional of local densities, which are determined
by the quasiparticle wavefunctions
$[u_{k\sigma}(\boldsymbol{r}),v_{k\sigma}(\boldsymbol{r})]$. The explicit
expressions for the most important ones are

$\displaystyle\rho(\mathbf{r})$
$\displaystyle=\sum_{k,\sigma}v^{*}_{k\sigma}(\mathbf{r})v_{k\sigma}(\mathbf{r}),$
(18) $\displaystyle\nu(\boldsymbol{r})$
$\displaystyle=\sum_{k}v^{*}_{k\uparrow}(\boldsymbol{r})u_{k\downarrow}(\boldsymbol{r}),$
(19) $\displaystyle\tau(\boldsymbol{r})$
$\displaystyle=\sum_{k,\sigma}\boldsymbol{\nabla}v^{*}_{k\sigma}(\mathbf{r})\cdot\boldsymbol{\nabla}v_{k\sigma}(\mathbf{r}),$
(20) $\displaystyle\boldsymbol{J}(\boldsymbol{r})$
$\displaystyle=\left.\frac{1}{2i}(\boldsymbol{\nabla}-\boldsymbol{\nabla}^{\prime})\times\boldsymbol{s}(\boldsymbol{r},\boldsymbol{r}^{\prime})\right|_{\boldsymbol{r}=\boldsymbol{r}^{\prime}},$
(21)

where

$\displaystyle\begin{split}s_{x}(\boldsymbol{r},\boldsymbol{r}^{\prime})&=\sum_{k}\left(v_{k\uparrow}^{*}(\boldsymbol{r})v_{k\downarrow}(\boldsymbol{r}^{\prime})+v_{k\downarrow}^{*}(\boldsymbol{r})v_{k\uparrow}(\boldsymbol{r}^{\prime})\right),\\\
s_{y}(\boldsymbol{r},\boldsymbol{r}^{\prime})&=i\sum_{k}\left(v_{k\uparrow}^{*}(\boldsymbol{r})v_{k\downarrow}(\boldsymbol{r}^{\prime})-v_{k\downarrow}^{*}(\boldsymbol{r})v_{k\uparrow}(\boldsymbol{r}^{\prime})\right),\\\
s_{z}(\boldsymbol{r},\boldsymbol{r}^{\prime})&=\sum_{k}\left(v_{k\uparrow}^{*}(\boldsymbol{r})v_{k\uparrow}(\boldsymbol{r}^{\prime})-v_{k\downarrow}^{*}(\boldsymbol{r})v_{k\downarrow}(\boldsymbol{r}^{\prime})\right).\\\
\end{split}$ (22)

The summations over $k$ should be performed for quasiparticle states with
quasiparticle energies $E_{k}$ that satisfy $0<E_{k}<E_{\mathrm{cut}}$, where
$E_{\mathrm{cut}}$ is the energy cutoff related to the momentum cutoff
$E_{\mathrm{cut}}=\hbar^{2}k^{2}_{\mathrm{cut}}/(2m)$ ($m$ stands for the mass
of nucleon). It should be chosen sufficiently large to ensure the convergence
of all observables. This convention is applied to all summations over $k$
throughout this paper. In the next section, we will show how to extract these
densities directly from the Green’s function without explicit diagonalization
of the HFB matrix. The method can be extended to other (not listed) densities.

###  II.2 Green’s function and local densities

Denoting the $4\times 4$ HFB matrix in Eq. (15) as $H$, the Green’s function
$G(\boldsymbol{r},\boldsymbol{r}^{\prime},z)$ is the solution of a matrix
equation,

$\displaystyle(zI_{4}-H)G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})=\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})I_{4},$
(23)

where $I_{4}$ stands for a $4\times 4$ unit matrix. Now the Green’s function
is a $4\times 4$ blocked matrix in the form:

$\displaystyle
G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})=\sum_{k}\frac{1}{z-E_{k}}\begin{pmatrix}u_{k\uparrow}(\boldsymbol{r})\\\
u_{k\downarrow}(\boldsymbol{r})\\\ v_{k\uparrow}(\boldsymbol{r})\\\
v_{k\downarrow}(\boldsymbol{r})\end{pmatrix}\cdot\begin{pmatrix}u_{k\uparrow}^{*}(\boldsymbol{r}^{\prime})\\\
u_{k\downarrow}^{*}(\boldsymbol{r}^{\prime})\\\
v_{k\uparrow}^{*}(\boldsymbol{r}^{\prime})\\\
v_{k\downarrow}^{*}(\boldsymbol{r}^{\prime})\end{pmatrix}^{T}.$ (24)

In the above equation the summation over $k$ includes all eigenstates. As we
discussed in Sec. II.1, we need to calculate the normal, anomalous, kinetic
and spin-orbit densities and to construct the HFB matrix in the self-
consistent iterations. A closer look at the explicit expressions for the
densities (18)-(22) reveals that we need to extract only 5 of 16 entries,
containing
$\\{v_{k\uparrow}(\boldsymbol{r}),v_{k\downarrow}(\boldsymbol{r})\\}\otimes\\{v_{k\uparrow}^{*}(\boldsymbol{r}^{\prime}),v_{k\downarrow}^{*}(\boldsymbol{r}^{\prime})\\}$
and
$v^{*}_{k\uparrow}(\boldsymbol{r})u_{k\downarrow}(\boldsymbol{r}^{\prime})$
products. In the next subsections, we provide prescriptions for extraction of
the local densities assuming that Green’s function can be efficiently
computed. For simplification of the formulas we introduce the notation

$\displaystyle
G_{\phi\rho,\psi\sigma}(z,\boldsymbol{r};\boldsymbol{r}^{\prime})=\sum_{k}\frac{\phi_{k\rho}(\boldsymbol{r})\psi^{*}_{k\sigma}(\boldsymbol{r}^{\prime})}{z-E_{k}}$
(25)

for submatrices of the Green’s function, where $\phi,\psi=\\{u,v\\}$ are
wavefunction coordinates and $\rho,\sigma=\\{\uparrow,\downarrow\\}$ are spin
coordinates.

####  II.2.1 Normal and anomalous density

The simplest products in the Green’s function are the normal density
$\rho(\boldsymbol{r})$ in Eq. (18) and the anomalous density
$\nu(\boldsymbol{r})$ in Eq. (19). For $\rho(\boldsymbol{r})$, one needs to
extract the $G_{v\uparrow,v\uparrow}$ and $G_{v\downarrow,v\downarrow}$
components from the $G$ matrix and add them. Next, by performing the contour
integral, like in Eq. (6), we obtain the normal density for a selected
$\boldsymbol{r}^{\prime}$ point

$\displaystyle\rho(\boldsymbol{r}^{\prime})=\frac{1}{2\pi
i}\sum_{\sigma}\oint_{C}dz\,\left.G_{v\sigma,v\sigma}(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\right|_{\boldsymbol{r}=\boldsymbol{r}^{\prime}},$
(26)

where the contour integral encompasses the interval $[0,E_{\mathrm{cut}}]$ of
the real axis. The expression for the anomalous density requires only the
$G_{u\downarrow,v\uparrow}$ component and by analogy reads

$\displaystyle\nu(\boldsymbol{r}^{\prime})=\frac{1}{2\pi
i}\oint_{C}dz\,\left.G_{u\downarrow,v\uparrow}(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\right|_{\boldsymbol{r}=\boldsymbol{r}^{\prime}}.$
(27)

The extension to spin densities (22) is straightforward.

####  II.2.2 Kinetic density

Without losing generality, we suppress the spin degrees of freedom in the
equations and write the (normal) Eq. (18) and the kinetic energy densities Eq.
(20) in a simplified form:

$\displaystyle\begin{split}\rho(\boldsymbol{r})&=\sum_{k}\psi^{*}_{k}(\boldsymbol{r})\psi_{k}(\boldsymbol{r}),\\\
\quad\tau(\boldsymbol{r})&=\sum_{k}\boldsymbol{\nabla}\psi^{*}_{k}(\boldsymbol{r})\cdot\boldsymbol{\nabla}\psi_{k}(\boldsymbol{r}).\end{split}$
(28)

In order to evaluate $\tau(\boldsymbol{r})$, we first calculate the Laplacian
of the number density:

$\displaystyle\Delta\rho(\boldsymbol{r})=2{\mathrm{R}e}[\sum_{k}\psi^{*}_{k}(\boldsymbol{r})\Delta\psi_{k}(\boldsymbol{r})+\boldsymbol{\nabla}\psi_{k}^{*}(\boldsymbol{r})\cdot\boldsymbol{\nabla}\psi_{k}(\boldsymbol{r})].$
(29)

Then $\tau(\boldsymbol{r})$ can be obtained via

$\displaystyle\tau(\boldsymbol{r})=\frac{1}{2}\Delta\rho(\boldsymbol{r})-\mathrm{Re}\left(\sum_{k}\psi^{*}_{k}(\boldsymbol{r})\Delta\psi_{k}(\boldsymbol{r})\right).$
(30)

Now the problem is reduced to the calculation of the quantity
$\sum_{k}\psi^{*}_{k}(\boldsymbol{r})\Delta\psi_{k}(\boldsymbol{r})$. Recall
that

$\displaystyle\begin{split}\sum_{k}\psi_{k}(\boldsymbol{r})\psi_{k}^{*}(\boldsymbol{r}^{\prime})=\frac{1}{2\pi
i}\oint_{C}dzG(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\end{split}$ (31)

and if we apply the Laplacian $\Delta$ on $\boldsymbol{r}$ on both sides of
the equation and set $\boldsymbol{r}=\boldsymbol{r}^{\prime}$ we obtain

$\displaystyle\begin{split}&\sum_{k}\Delta\psi_{k}(\boldsymbol{r})\left.\psi_{k}^{*}(\boldsymbol{r}^{\prime})\right|_{\boldsymbol{r}=\boldsymbol{r}^{\prime}}\\\
&=\frac{1}{2\pi i}\oint_{c}dz\left.\Delta
G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\right|_{\boldsymbol{r}=\boldsymbol{r}^{\prime}}.\end{split}$
(32)

Computation of the derivatives of the Green’s functions introduces significant
numerical costs. We avoid the calculation of the Laplacian by taking advantage
of an integration by parts

$\displaystyle\begin{split}\left.\Delta
G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\right|_{\boldsymbol{r}=\boldsymbol{r}^{\prime}}&=\int
d\boldsymbol{r}\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})\Delta
G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\\\ &=\int
d\boldsymbol{r}\left(\Delta\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})\right)G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\end{split}$
(33)

where the second equality uses the boundary condition
$\lim_{r\to\infty}G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})=0$. Finally

$\displaystyle\begin{split}\sum_{k}\psi^{*}_{k}(\boldsymbol{r}^{\prime})&\Delta\psi_{k}(\boldsymbol{r}^{\prime})=\\\
=&\frac{1}{2\pi i}\oint_{C}dz\int
d\boldsymbol{r}\left(\Delta\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})\right)G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\end{split}$
(34)

which is a contour integral of a convolution. Therefore, to obtain the total
kinetic energy density Eq. (20), one just needs to extract
$G_{v\uparrow,v\uparrow}$ and $G_{v\downarrow,v\downarrow}$ components from
the $G$ matrix, and calculate the kinetic energy density
$\tau=\tau_{\uparrow}+\tau_{\downarrow}$ based on Eqs. (30) and (34). When
implemented numerically $\Delta\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})$
is the numerical implementation of the Laplacian applied to the
$\delta$-function on the lattice.

####  II.2.3 Spin-current density

The spin-current density Eq. (21) can be written in the explicit form:

$\displaystyle\begin{split}J_{x}(\boldsymbol{r})&=\mathrm{Im}\sum_{k}\left(v^{*}_{k\downarrow}(\boldsymbol{r})\frac{\partial}{\partial
y}v_{k\downarrow}(\boldsymbol{r})-v^{*}_{k\uparrow}(\boldsymbol{r})\frac{\partial}{\partial
y}v_{k\uparrow}(\boldsymbol{r})\right)+\mathrm{Re}\sum_{k}\left(v^{*}_{k\uparrow}(\boldsymbol{r})\frac{\partial}{\partial
z}v_{k\downarrow}(\boldsymbol{r})-v^{*}_{k\downarrow}(\boldsymbol{r})\frac{\partial}{\partial
z}v_{k\uparrow}(\boldsymbol{r})\right),\\\
J_{y}(\boldsymbol{r})&=\mathrm{Im}\sum_{k}\left(v^{*}_{k\uparrow}(\boldsymbol{r})\frac{\partial}{\partial
x}v_{k\uparrow}(\boldsymbol{r})+v^{*}_{k\downarrow}(\boldsymbol{r})\frac{\partial}{\partial
x}v_{k\downarrow}(\boldsymbol{r})-v^{*}_{k\downarrow}(\boldsymbol{r})\frac{\partial}{\partial
z}v_{k\uparrow}(\boldsymbol{r})-v^{*}_{k\uparrow}(\boldsymbol{r})\frac{\partial}{\partial
z}v_{k\downarrow}(\boldsymbol{r})\right),\\\
J_{z}(\boldsymbol{r})&=\mathrm{Im}\sum_{k}\left(v^{*}_{k\downarrow}(\boldsymbol{r})\frac{\partial}{\partial
y}v_{k\uparrow}(\boldsymbol{r})+v^{*}_{k\uparrow}(\boldsymbol{r})\frac{\partial}{\partial
y}v_{k\downarrow}(\boldsymbol{r})\right)+\mathrm{Re}\sum_{k}\left(v^{*}_{k\downarrow}(\boldsymbol{r})\frac{\partial}{\partial
x}v_{k\uparrow}(\boldsymbol{r})-v^{*}_{k\uparrow}(\boldsymbol{r})\frac{\partial}{\partial
x}v_{k\downarrow}(\boldsymbol{r})\right),\end{split}$ (35)

which are combinations of the quantities
$\sum_{k}\psi^{*}_{k}\boldsymbol{\nabla}\psi_{k}$ for different spin
combinations. The evaluation procedure in this case is similar to the one we
used for the density $\sum_{k}\psi^{*}_{k}\Delta\psi_{k}$ and we obtain

$\displaystyle\begin{split}&\sum_{k}\psi^{*}_{k}(\boldsymbol{r}^{\prime})\boldsymbol{\nabla}\psi_{k}(\boldsymbol{r}^{\prime})=\\\
&=-\frac{1}{2\pi i}\oint_{C}dz\int
d\boldsymbol{r}\left(\boldsymbol{\nabla}\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})\right)G(z,\boldsymbol{r};\boldsymbol{r}^{\prime}),\end{split}$
(36)

where $\boldsymbol{\nabla}$ represents the gradient operator on
$\boldsymbol{r}$. One needs to extract all spin components $G_{v\rho,v\sigma}$
from the $G$ matrix and form the appropriate combinations.

##  III Numerical implementation

###  III.1 Shifted COCG method

Throughout our theoretical framework, the core problem is to solve Eq. (23)
for all $z$’s on the contour. We solve these equations separately for each
coordinate point $\boldsymbol{r}^{\prime}$. When discretized, this problem is
reduced to the linear equations for a given set of contour points $z_{m}$
$(m=0,1,\cdots,m_{\text{max}})$. These linear equations are called “shifted”
linear equations or shifted linear systems because the matrices $z_{m}-H$ are
connected by scalar shifts. We solve this problem by using the shifted COCG
method Yamamoto  _et al._ (2008); Takayama  _et al._ (2006), which is an
iterative method for solving large-scale shifted linear systems with symmetric
matrices. The details of this algorithm are available in Ref. Yamamoto  _et
al._ (2008); Takayama  _et al._ (2006) and here we give only a brief review
and illustrate how to apply it to the HFB matrix in Eq. (15).

For a given symmetric matrix $A$, we want to solve the linear equation:

$\displaystyle A\boldsymbol{x}=\boldsymbol{b},$ (37)

which is called the reference system, and its shifted equations

$\displaystyle(A+\sigma I)\boldsymbol{x}^{\sigma}=\boldsymbol{b},$ (38)

where $\sigma$ is a scalar complex shift, and $\boldsymbol{x}^{\sigma}$ is the
solution of the corresponding shifted system. The reference system is solved
by the COCG method van der Vorst and Melissen (1990). We define
$\boldsymbol{x}_{n}$ as the approximate solution in the $n$-th iteration,
$\boldsymbol{r}_{n}$ as the corresponding residual vector
$\boldsymbol{r}_{n}=\boldsymbol{b}-A\boldsymbol{x}_{n}$, search direction
vector $\boldsymbol{p}_{n}$, and other coefficients $\alpha_{n},\beta_{n}$.
With the initial condition
$\boldsymbol{x}_{0}=\boldsymbol{0},\boldsymbol{r}_{0}=\boldsymbol{b},\boldsymbol{p}_{0}=\boldsymbol{b},\alpha_{-1}=1,\beta_{-1}=0$
we have to perform the following iterations:

$\displaystyle\alpha_{n-1}$
$\displaystyle=\frac{\boldsymbol{r}_{n-1}^{T}\cdot\boldsymbol{r}_{n-1}}{\boldsymbol{p}_{n-1}^{T}\cdot
A\boldsymbol{p}_{n-1}},$ (39) $\displaystyle\boldsymbol{x}_{n}$
$\displaystyle=\boldsymbol{x}_{n-1}+\alpha_{n-1}\boldsymbol{p}_{n-1},$ (40)
$\displaystyle\boldsymbol{r}_{n}$
$\displaystyle=\boldsymbol{r}_{n-1}-\alpha_{n-1}A\boldsymbol{p}_{n-1},$ (41)
$\displaystyle\beta_{n-1}$
$\displaystyle=\frac{\boldsymbol{r}_{n}^{T}\cdot\boldsymbol{r}_{n}}{\boldsymbol{r}_{n-1}^{T}\cdot\boldsymbol{r}_{n-1}},$
(42) $\displaystyle\boldsymbol{p}_{n}$
$\displaystyle=\boldsymbol{r}_{n}+\beta_{n-1}\boldsymbol{p}_{n-1}$ (43)

where the $T$ represents transpose only, without complex conjugation and
$\cdot$ stands for the implied scalar product between left and right vectors.
Note, the evaluation of $A\boldsymbol{p}_{n-1}$ is the most computationally
consuming part of an iteration, as it requires a matrix-vector multiplication.
The residual vector measures the accuracy of solution in the $n$-th iteration
and it is used as a breaking condition for iterations.

Shifted systems Eq. (38) can also be solved by COCG algorithm. For each
shifted equation we introduce the corresponding vectors
$\boldsymbol{x}_{n}^{\sigma}$, $\boldsymbol{r}_{n}^{\sigma}$, and
$\boldsymbol{p}_{n}^{\sigma}$ and coefficients $\alpha_{n}^{\sigma}$ and
$\beta_{n}^{\sigma}$, and initialize them with the same initial conditions as
the reference system. However, the iterations (39)-(43) for the shifted
systems can benefit from their collinearity, meaning that the residuals
vectors for the reference and shifted systems are connected Yamamoto  _et al._
(2008)

$\displaystyle\boldsymbol{r}_{n}^{\sigma}=\frac{1}{\pi^{\sigma}_{n}}\boldsymbol{r}_{n},$
(44)

where the proportionality constant in each iteration is given by
($\pi_{0}^{\sigma}=\pi_{-1}^{\sigma}=1$)

$\displaystyle\begin{split}\pi_{n}^{\sigma}=\left(1+\alpha_{n-1}\sigma+\alpha_{n-1}\frac{\beta_{n-2}}{\alpha_{n-2}}\right)\pi_{n-1}^{\sigma}\\\
-\alpha_{n-1}\frac{\beta_{n-2}}{\alpha_{n-2}}\pi_{n-2}^{\sigma}.\end{split}$
(45)

Thus, for the shifted systems instead of evaluating the time consuming Eq.
(41) one can evaluate the simpler Eq. (44). Moreover, the coefficients
$\alpha^{\sigma}$ and $\beta^{\sigma}$ are also connected with corresponding
coefficients of the reference system

$\displaystyle\alpha_{n-1}^{\sigma}=\frac{\pi_{n-1}^{\sigma}}{\pi^{\sigma}_{n}}\alpha_{n-1},$
(46)
$\displaystyle\beta_{n-1}^{\sigma}=\left(\frac{\pi_{n-1}^{\sigma}}{\pi_{n}^{\sigma}}\right)^{2}\beta_{n-1}.$
(47)

These equations replace Eqs. (39) and (42). Often the reference system is
called the seed system as it “seeds” data for shifted systems. The iterations
are executed simultaneously for all systems and they end when the desired
accuracy is achieved for all equations. Since when evaluating the density we
set $\boldsymbol{r}=\boldsymbol{r}^{\prime}$ in the Green’s function in the
vector Eq. (44), one needs to evaluate only one component, and not the entire
vector; see also Sec. III.5.

Assume for a moment that the HFB matrix is symmetric. This is the case when
the spin-orbit term in the nuclear density functional is ignored and the
pairing potential is real. Then we can straightforwardly let $A=z_{0}I-H$,
where $z_{0}$ is selected point from contour. Shifts are given by
$\sigma_{m}=z_{m}-z_{0}$. We solve the problem on a Cartesian mesh grid with a
lattice size $N_{x}\times N_{y}\times N_{z}$ and lattice spacing $dx=dy=dz$.
After discretization, the HFB matrix has dimensions $N\times N$ with
$N=4N_{x}N_{y}N_{z}$, but in this case spin-up and spin-down components are
decoupled and the matrix size is effectively $2N_{x}N_{y}N_{z}$. For a fixed
spatial point $\boldsymbol{r}^{\prime}$ function
$\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})$ becomes vector, where all $N$
elements are zeros except one [$=1/(dxdydz)$], corresponding to the selected
position. We denote this vector by
$\delta_{\boldsymbol{r}^{\prime}}(\boldsymbol{r})$. Since we do not need all
elements of $G$ matrix, it is sufficient to solve problems
$A\boldsymbol{x}_{1}=\boldsymbol{b}_{1}$ and
$A\boldsymbol{x}_{2}=\boldsymbol{b}_{2}$ (together with shifted counterparts)
where $\boldsymbol{b}_{1}=(0,0,\delta_{\boldsymbol{r}^{\prime}},0)^{T}$ and
$\boldsymbol{b}_{2}=(0,0,0,\delta_{\boldsymbol{r}^{\prime}})^{T}$. Then
solutions $\boldsymbol{x}_{1}$ and $\boldsymbol{x}_{2}$ are third and fourth
columns of Eq. (24).

In general, the HFB matrix is a hermitian but not a symmetric matrix, which
cannot be solved using the shifted COCG method directly. However, one can use
the COCG method designed for symmetric matrices in the case of Hermitian
matrices by performing a simple matrix transformation. Assume that $H$ is an
arbitrary Hermitian matrix which can be divided into real and imaginary parts:
$H=H_{x}+iH_{y}$, with obvious symmetry properties $H_{x}=H_{x}^{T}$ and
$H_{y}=-H_{y}^{T}$. If we also divide the eigenvectors of $H$ into real and
imaginary parts as $\psi=x+iy$. The eigenvalue problem $H\psi=\lambda\psi$
will be converted to:

$\displaystyle(H_{x}+iH_{y})(x+iy)=\lambda(x+iy).$ (48)

After collecting the real and imaginary terms it will become an equivalent
eigenvalue problem:

$\displaystyle\begin{pmatrix}H_{x}&-H_{y}\\\
H_{y}&H_{x}\end{pmatrix}\begin{pmatrix}x\\\
y\end{pmatrix}=\lambda\begin{pmatrix}x\\\ y\end{pmatrix},$ (49)

but now the matrix on the left hand side is real symmetric. The Green’s
function $G=G_{x}+iG_{y}$ of a general HFB matrix $H$ is then the solution of
the equation [equivalent to Eq. (23)]

$\displaystyle\left(zI-\right.H^{\prime}\left.\right)\begin{pmatrix}G_{x}\\\
G_{y}\end{pmatrix}=\begin{pmatrix}\delta^{(4)}\\\ 0\end{pmatrix},$ (50)

where

$\displaystyle H^{\prime}=\begin{pmatrix}H_{x}&-H_{y}\\\
H_{y}&H_{y}\end{pmatrix},\quad\delta^{(4)}=\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})\otimes
I_{4}.$ (51)

As before we need to solve two linear equations
$A^{\prime}\boldsymbol{x}^{\prime}_{1}=\boldsymbol{b}^{\prime}_{1}$ and
$A^{\prime}\boldsymbol{x}_{2}=\boldsymbol{b}^{\prime}_{2}$ (and their shifted
equation) with $A^{\prime}=zI-H^{\prime}$ in doubled space. Now, vectors of
length $2N$ need to be set for
$\boldsymbol{b}^{\prime}_{1}=(0,0,\delta_{\boldsymbol{r}^{\prime}},0,0,0,0,0)^{T}$
and
$\boldsymbol{b}^{\prime}_{2}=(0,0,0,\delta_{\boldsymbol{r}^{\prime}},0,0,0,0)^{T}$.

###  III.2 The form of the integration contour

Before a further discussion of the shifted COCG method, we should provide a
clear definition of the contour, which determines the shifted systems and is
fundamental for the convergence behavior and error control of the whole
algorithm. We need to calculate an integral of the form $\frac{1}{2\pi
i}\oint_{C}dzf(z)$ where the contour encloses exactly a segment of the real
axis $[0,E_{\mathrm{cut}}]$. It is natural to choose the contour to be
symmetric with respect to the real axis since the poles (the HFB spectrum) are
real. We can parametrize the contour as $z(\varphi)$ and a simple choice of
its form is an ellipse:

$\displaystyle\begin{split}z(\varphi)&=\frac{E_{\mathrm{cut}}}{2}+\frac{E_{\mathrm{cut}}}{2}\cos\varphi+ih\sin\varphi\\\
&=E_{\mathrm{cut}}\cos^{2}\frac{\varphi}{2}+ih\sin\varphi\end{split}$ (52)

where $0\leq\varphi\leq 2\pi$ and the height of ellipse $h$ is the parameter
that should be chosen carefully, see Fig. 1.

Figure 1:  The ellipse contour $C:z(\varphi)$ in Eq. (52) 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).

The integrand function behaves smoothly, only for points that are far away
from real axis. As the contour approaches real axis, the function
$f(z)\sim\sum_{k}\frac{1}{z-E_{k}}$ as a function of $z$ starts to exhibit
fast oscillations, making accurate numerical integration difficult. For this
reason, we choose non-uniformly distributed integration points along the
contour. The density of points distribution depends on distance from the real
axis with more points closer to the real axis, and less far away from the real
axis. The height $h$ of the eclipse is set to be significantly larger than the
expected average separation between poles, and we compute the integrand
function for angles $\varphi$ given by a distribution function ($u$ and the
$l$ subscript corresponds to upper and lower parts of the contour with respect
to the real axis)

$\displaystyle\begin{split}\varphi_{u}&=\frac{\pi}{2}\left\\{1+\tanh\left[\alpha\tan\left(\phi-\frac{\pi}{2}\right)\right]\right\\},\\\
\varphi_{l}&=\varphi_{u}+\pi,\quad 0\leq\phi\leq\pi,\end{split}$ (53)

where $\alpha$ is a parameter (usually $\alpha\sim$1-5) and the integral over
$\phi$ is discretized with a set of evenly distributed points $\phi_{n}$ in
the interval $[0,\pi]$ with step size $\Delta\phi$. Thus the contour integral
can be converted into an weighted sum

$\displaystyle\frac{1}{2\pi
i}\oint_{C}dzf(z)\cong\sum_{\sigma=u,l;n}f\\{z[\varphi_{\sigma}(\phi_{n})]\\}w_{\sigma}(\phi_{n})\Delta\phi.$
(54)

Typically, the number of integration points does not exceed $10^{4}$, and
further increases do not improve the numerical result. The computation of the
integrand along the contour is very fast and represents a very small fraction
of the entire calculation time, which is dominated by the calculation of the
Green’s function at the reference point.

###  III.3 Convergence behavior

(a)  (b)  (c)

(d)  (e)  (f)

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\;\textrm{fm}$. Panels (a) and (d) show the distribution of
the norms of the residual vectors $\lVert r_{n}^{\sigma}\rVert$ on the semi-
ellipse contour $z(\varphi)$ 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 $\lVert
r_{n}^{\sigma}\rVert$ for three points $z(0)$ (red), $z(\pi/2)$ (green) and
$z(\pi)$ (blue) as a function of iteration number. Panels (c) and (f) show the
difference between the value obtained from COCG method in the $n$-th iteration
and the exact value for the normal (red) and anomalous (blue) densities.

The norm of the residual vectors $\boldsymbol{r}_{n}$ at the reference and of
the shifted systems $\boldsymbol{r}_{n}^{\sigma}$ quantify the error of the
approximate solution $\boldsymbol{x}_{n}$ and $\boldsymbol{x}_{n}^{\sigma}$ in
each iteration. Since we need to calculate the Green’s functions
$G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})$ for all $z$’s on the contour
defined in Sec. III.2, the convergence behavior on the contour needs to be
studied.

As a simple but non-trivial test case, we consider only neutrons, and choose a
phenomenological Wood-Saxon (W-S) model Ring and Schuck (2004); Bohr and
Mottelson (1998). The central part of the single-particle Hamiltonian Eq. (16)
is given by $U(\boldsymbol{r})=V_{0}f(\boldsymbol{r})$, where
$f(\boldsymbol{r})$ has the symmetrized W-S form:

$\displaystyle f(\boldsymbol{r})$
$\displaystyle=\frac{1-\exp\left(\frac{-2R_{0}}{a}\right)}{\left[1+\exp\left(\frac{r-R_{0}}{a}\right)\right]\left[1+\exp\left(\frac{-r-R_{0}}{a}\right)\right]}$
(55)

and $V_{0}=-50~{}\mathrm{MeV}$, $r=\sqrt{x^{2}+y^{2}+(z/\beta)^{2}}$,
$R_{0}=1.12A_{0}^{1/3}-0.86A_{0}^{-1/3}~{}\mathrm{fm}$, $A_{0}=100$,
$a=0.54~{}\mathrm{fm}$. By varying the parameter $\beta$, we calculate the
densities in a W-S potential in both spherical and axially deformed cases. The
spin-orbit potential $\boldsymbol{W}(\boldsymbol{r})$ is chosen as
$\boldsymbol{W}(\boldsymbol{r})=\lambda\boldsymbol{\nabla}U(\boldsymbol{r})$
where $\lambda=-0.5~{}\mathrm{fm}^{2}$, and the effective mass is set to be
the bare mass of neutron $mc^{2}=939.565~{}\mathrm{MeV}$. The self-consistent
pairing field $\Delta(\boldsymbol{r})$ in Eq. (17) should be determined by the
anomalous density $\nu(\boldsymbol{r})$, but here we still use a
phenomenological W-S potential shape
$\Delta(\boldsymbol{r})=12/A_{0}^{1/2}f(\boldsymbol{r})$. In the test, the
chemical potential is fixed at $\mu=-7~{}\mathrm{MeV}$. The numerical tests
are performed in a cubic box of size $L_{x}=L_{y}=L_{z}=20~{}\mathrm{fm}$ for
the spherical W-S model ($\beta=1$), and a rectangular box of size
$L_{x}=L_{y}=20,L_{z}=25~{}\mathrm{fm}$ for the deformed W-S model
($\beta=1.5$). In both cases, the lattice constant
$dx=dy=dz=1.25~{}\mathrm{fm}$. The number of lattice points is
$N_{x}=L_{x}/dx$ and respectively for $y$ and $z$ direction. The energy cutoff
is set to $E_{\mathrm{cut}}=100~{}\mathrm{MeV}$.

We find that the convergence properties of the algorithm are strongly
correlated with the density of particles for a given point. For points with
higher particle densities the method needs more iterations to converge. This
behavior can be easily understood if we recall that the starting point for
iterations is $\boldsymbol{x}_{0}=0$; see Sec. III.1. For the points where the
density is higher the convergence is noticeably slower. In the following
tests, we consider the convergence properties of the method for the points
requiring the highest number of iterations. These points are located close to
the center of the simulation box.

Results for the spherically symmetric case ($\beta=1$) are presented in the
top row of Fig. 2. The Figure 2(a) shows the distribution of
$\lVert{\boldsymbol{r}}_{n}^{\sigma}\rVert$ on the semi-ellipse
$z(\varphi),0\leq\varphi\leq\pi$ on the upper half complex plane, for
iteration number $n=4000$, and fixed position $\boldsymbol{r}^{\prime}$
located close to center of the box. In the lower half plane (
$\pi\leq\varphi\leq 2\pi$ ) the behavior is identical due to the reflection
symmetry of the integration contour. The shifted system converges most quickly
in the middle part of the contour (represented by color green), more slowly
near the left end of contour $z(\pi)=0$ (represented by color blue), and very
slowly near the right end of the contour $z(0)=E_{\mathrm{cut}}$ (represented
by the color red). A closer look at the convergence behavior at the
representative points in these three parts, respectively $z(\pi/2)$, $z(\pi)$,
and $z(0)$, is shown on Figure 2 (b). Besides the rapid convergence of the
middle point $z(\pi/2)$, the residual at the origin $z(\pi)=0$ also has a
stable decrease with iterations. However, the iterations for the right end
$z(0)=E_{\mathrm{cut}}$ fail to converge, and $\lVert r_{n}^{\sigma}\rVert$
keeps oscillating around $10^{-2}$. Consider the convergence analysis of the
conjugate gradient (CG) method Trefethen and Bau (1997) in which case the
convergence ratio depends on the 2-norm condition number $\kappa$ of the
matrix $A$:

$\displaystyle\frac{\lVert r_{n}\rVert}{\lVert r_{0}\rVert}\leq
2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{n},$ (56)

where the condition number
$\kappa(A)=\lvert\lambda(A)\rvert_{\mathrm{max}}/\lvert\lambda(A)\rvert_{\mathrm{min}}$
is the ratio of maximum and minimal absolute value of the eigenvalues of
matrix $A$. From Eq. (56) it is clear that a larger $\kappa$ leads to a slower
convergence. The order of magnitude of the condition number for matrix
$A=zI-H$ can be estimated as $\frac{\max|z-E_{k}|}{\min|z-E_{k}|}$, where
$E_{k}$ are eigenvalues of the HFB matrix from interval
$[0,E_{\mathrm{cut}}]$. Definitely, the condition number decreases as the
imaginary part of $z$ increases, thus confirming the finding that the
convergence is fastest in middle part of the semi-ellipse. The largest values
of $\kappa$ occur at points close to the real axis. We know that the spectrum
of a HFB matrix is discrete if $\lvert E_{k}\rvert<\lvert\mu\rvert$ and
continuous otherwise Bulgac ; Dobaczewski  _et al._ (1984). Around $z=0$ there
is a gap in the spectrum and therefore for $z\rightarrow 0$
$\kappa=O(E_{\mathrm{cut}}/\Delta)$, where $\Delta$ is average value of the
pairing gap in the system. When $z\rightarrow E_{\mathrm{cut}}$,
$\min|z-E_{k}|$ can be arbitrarily small (only set by lattice resolution), and
in the limit of zero lattice spacing $dx$ the condition number diverges. This
“quasi-singular” matrix equation is hard to solve within a reasonable
iteration number. However, these continuous states near the energy cutoff have
very small occupation probabilities ($<10^{-5}$) and make negligible
contributions to the local densities. While it is possible that
preconditioning the COCG method can help Gu  _et al._ (2014); Li  _et al._
(2011, 2008); Wang  _et al._ (2009) our efforts in this direction did not
yield any advantages.

Figure 2(c) tracks the error of the particle density $\rho$ and anomalous
density $\nu$ obtained by the calculated Green’s function in each iteration.
The exact value of these two densities $\rho_{0}$ ($\sim 10^{-1}$) and
$\nu_{0}$ ($\sim 10^{-3}$) are calculated by direct diagonalization of the
same HFB matrix via ScaLAPACK. It shows that the particle density $\rho$ can
converge to an accuracy of $10^{-8}$ within 1000 iterations, and this accuracy
will not improve as the iteration goes past that. The anomalous density $\nu$
will converge to an accuracy of $10^{-6}$ within the same iteration count,
which is quite acceptable, since the contribution to the energetics of a
nucleus from the pairing field is much smaller in absolute terms. The pairing
field is almost two orders of magnitude smaller than the single-particle
potential. This behavior of the anomalous density is due to the fact that the
high-energy quasiparticle continuum states contribute with significant weights
as the anomalous density diverges when the energy cutoff is increased Bulgac ;
Bulgac and Yu (2002); Borycki  _et al._ (2006). The same tests for deformed
W-S model are presented in Figure. 2(d)–2(f). They show a similar pattern of
the residual distribution on the contour. The local densities can reach the
same accuracy but with more iterations, because of the slower convergence of
points near $E=0$. In both the spherical and deformed cases, the convergence
behavior at $E=0$ can represent the convergence behavior of the final local
densities, and thus is chosen to be the breaking condition for the iterations.

(a)  (b)

(c)  (d)

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_{\textrm{cocg}}-\rho_{\textrm{diag.}}|$, and the right
column shows them for the anomalous density
$\delta\nu=|\nu_{\textrm{cocg}}-\nu_{\textrm{diag.}}|$.

Until now our discussion of the convergence behavior has been for a fixed
spatial point in the lattice system. To have a comprehensive view of the
accuracy of this method, we also compare the difference between the local
densities $\rho$ and $\nu$ calculated with the shifted COCG method and with
the direct diagonalization method on all lattice spatial points. Figure. 3
displays the differences for spherical and deformed W-S model along the $z$
axis for fixed $x=0$ and $y=0$ coordinates. Because the local densities on
this spatial line have the largest magnitude throughout the box, we will also
get the largest errors on this line. To study the influence of the lattice
constants on the accuracy, besides the previous $dx=1.25$, we also performed
calculations with smaller lattice constants $dx=1.0$ and $dx=0.8333$. From the
figures one can see the global maximum error of $\rho$ and $\nu$ are
respectively $10^{-8}$ and $10^{-6}$, which matches the convergence behavior
of Fig. 2. The kinetic and spin-current densities are calculated with similar
accuracy as achieved for the calculation of normal density. Moreover, the
accuracy is little affected by the value of the lattice constant, a fact which
can be expected as argued in Sec. I; see also Refs. Bulgac and Forbes (2013);
Ryssens  _et al._ (2015b).

###  III.4 Exploiting symmetries

Densities on each spatial point in the system are calculated independently.
The number of points to be processed is equal to the number of lattice points
$N_{x}N_{y}N_{z}$. This number can be easily reduced if the system exhibits
symmetries, like reflection symmetries, axial symmetry, or spherical symmetry.
For axially or spherically symmetric systems, one can argue that it is more
profitable to exploit this symmetry directly on the level of HFB equations,
i.e., assume the correct symmetry for the wavefunctions and solve the
constrained HFB problem. However, if the solution of HFB problem will be used
as initial point for an unconstrained 3D simulation, then this approach is not
accurate enough. Typically, the solutions obtained from a solver which
explicitly uses spherical or axial symmetries once discretized on a 3D spatial
lattice are no longer orthogonal to each other and do not represent self-
consistent eigenstates (with the necessary numerical accuracy for a stable
numerical integration in time). Self-consistent iterations in full 3D space
are still required to get a high quality state. We demonstrate the utilization
of axial symmetry for reducing the computational cost, while the underlying
HFB matrix is defined in the full 3D coordinate space.

Consider the scalar local densities represented in cylindrical coordinates
$(r,\phi,z)$:

$\displaystyle\rho(\boldsymbol{r})$ $\displaystyle=\rho(r,z),$
$\displaystyle\nu(\boldsymbol{r})$ $\displaystyle=\nu(r,z),$
$\displaystyle\tau(\boldsymbol{r})$ $\displaystyle=\tau(r,z).$

We only need to calculate the densities on the points with different values of
$(r,z)$’s in the system. This reduces the number of points to be calculated to
approximately $1/8$ of the total number of points in the full 3D lattice
system, as only the points with $0\leq x\leq y\leq L/2~{}(L=L_{x}=L_{y})$ need
to be explicitly considered. Thus the performance will improve by a factor of
8. A reflection symmetry along the $z$-axis can add another factor of 2, thus
an overall speedup of 16.

We re-expand the vector density $\boldsymbol{J}(\boldsymbol{r})$ along the
three unit directions of the cylindrical coordinate system as $J_{r}$,
$J_{\phi}$, $J_{z}$, which can be related to the Cartesian components $J_{x}$,
$J_{y}$, $J_{z}$ by the transformation

$\displaystyle J_{r}(r,z)$ $\displaystyle=J_{x}\cos\phi+J_{y}\sin\phi,$
$\displaystyle J_{\phi}(r,z)$ $\displaystyle=J_{y}\cos\phi-J_{x}\sin\phi,$
(57) $\displaystyle J_{z}(r,z)$ $\displaystyle=J_{z},$

and reversely

$\displaystyle J_{x}$ $\displaystyle=J_{r}\cos\phi-J_{\phi}\sin\phi,$
$\displaystyle J_{y}$ $\displaystyle=J_{r}\sin\phi+J_{\phi}\cos\phi,$ (58)
$\displaystyle J_{z}$ $\displaystyle=J_{z}.$

The components $J_{r}$, $J_{\phi}$, $J_{z}$ have axial symmetry, i.e., all of
them depend only on $r$ and $z$, while the Cartesian components $J_{x}$,
$J_{y}$, and $J_{z}$ depend on all three spatial coordinates. In our approach,
we first calculate the Cartesian components $J_{x}$, $J_{y}$, $J_{z}$ for the
set of points with different $(r,z)$, we transform them into the cylindrical
components $J_{r}$, $J_{\phi}$, $J_{z}$ via Eqs. (III.4), and using the axial
symmetry we cover all the equivalent points in Cartesian space via reverse
transformation, using Eqs. (III.4).

Numerically, we found that the axial symmetry is not strictly realized due to
the finite size of the lattice constant. In the tests discussed above for a
spherical or deformed W-S model, we start from a HFB Hamiltonian with strict
axial symmetry, and calculate the local densities on all points of the lattice
system without enforcing any symmetry constraints. In theory, for example, the
scalar densities ($\rho$, $\nu$, $\tau$) on the point of Cartesian coordinates
(in fm): $(0,5,0)$ should have the same value ad on the point $(3,4,0)$, but
numerically they are different due to the cubic symmetry of the simulation
box. The relative difference between the scalar densities on these pairs of
pseudo-equivalent points varies from $10^{-2}$ to $10^{-4}$ as the lattice
constants decreases from $dx=1.25~{}\mathrm{fm}$ to $dx=0.8333~{}\mathrm{fm}$.
In order to fix this issue, we include all points in these pairs into the
calculation, with a small increase of the total number of lattice points
needed to be calculated.

This lattice effect becomes more pronounced when one calculates the
cylindrical components of the vector density $\boldsymbol{J}(\boldsymbol{r})$.
The components $J_{r}$ and $J_{z}$ can be treated in the same manner as the
scalar densities, but $J_{\phi}$, which is expected theoretically to vanish,
suffers from a numerical noise, which oscillates around zero. The average
magnitude of this noise decreases from $10^{-5}$ to $10^{-7}$ as $dx$
decreases from $1.25~{}\mathrm{fm}$ to $0.8333~{}\mathrm{fm}$ and for
simplicity we just force them to be zero, which will bring some discrepancy in
the calculation of $\boldsymbol{J}(\boldsymbol{r})$ when comparing the COCG
and the direct diagonalization approaches. We will discuss this discrepancy in
realistic calculations for finite nuclei in Sec. IV and show that it vanishes,
as naturally expected, as the lattice constant goes to 0.

###  III.5 Computational cost

The computational cost of the shifted COCG iteration is set by the cost of
solving the reference system, as it involves matrix-vector (MV) multiplication
between the matrix $A=z_{0}-H$ and the vector $\boldsymbol{p}_{n}$ of size
$2N$, where $N=4N_{x}N_{y}N_{z}$. The non-locality in the HFB matrix (see the
Appendix for explicit form) comes from the Laplacian or gradient operators.
When discretized, the gradient or the Laplacian of a function (vector) can be
obtained either via a finite difference formula or through the Fourier
transform. In both cases the operation is represented by a sparse matrix.
Exploiting the sparsity accelerates the MV operations significantly. We
compute derivatives using Fourier transforms, as it is a more accurate method
than using finite difference formulas Ryssens  _et al._ (2015b). Moreover, we
advocate the direct use of fast Fourier transforms (FFT) due to a manageable
complexity $O(N\log N)$. The computation of derivatives with FFT is expected
to be faster than accurate multipoint finite difference algorithms Aichinger
and Krotscheck (2005). For shifted systems, the collinear theorem insures that
there is no need to evaluate MV products. The vector-scalar arithmetic in Eqs.
(40), (43) and (44) for $\sigma\neq 0$ makes a considerable contribution to
the total computational cost if all $2Nm$ elements of
$\boldsymbol{x}^{\sigma}_{n}$, $\boldsymbol{p}^{\sigma}_{n}$,
$\boldsymbol{r}^{\sigma}_{n}$ are calculated, where $m$ is the number of
points on the contour. However, in order to obtain the local densities on a
fixed spatial point $\boldsymbol{r}^{\prime}$, we do not need to know the
Green’s function $G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})$ for all spatial
points $\boldsymbol{r}$ in the system. This is because integrand function
$f(z)$ requires either
$\left.G(z,\boldsymbol{r};\boldsymbol{r}^{\prime})\right|_{\boldsymbol{r}=\boldsymbol{r}^{\prime}}$
or overlap of Green’s function with
$\boldsymbol{\nabla}\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})$ or
$\Delta\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})$. Derivatives of
$\delta(\boldsymbol{r}-\boldsymbol{r}^{\prime})$ have only
$l=N_{x}+N_{y}+N_{z}-2$ non-zero elements (which are located on three lines in
$x,y,z$ directions that cross at $\boldsymbol{r}^{\prime}$). The vector-scalar
arithmetic can be executed only for elements required by the integrand
functions. Finally, the size of $\boldsymbol{x}^{\sigma}_{n}$,
$\boldsymbol{p}^{\sigma}_{n}$, $\boldsymbol{r}^{\sigma}_{n}$ is reduced to the
order $2lm$, which is a small number when compared with $N$ for the reference
system. The cost related to solving the shifted systems turns out to be
negligible. Therefore, the theoretical upper bound of the computational cost
within one full self-consistent iteration of shifted COCG method at $N$
coordinates is determined by the FFT operation [$O(N\log N)$] and the maximum
number of COCG Krylov iterations $M$ required to converge each point in the
reference system. The number of Krylov iterations to converge is problem
dependent and theoretically it scales as $O(N)$, but in all of our realistic
calculation we find that $M\approx 4,000$ is sufficient to get a high accuracy
solution for nuclear problems. Final numerical complexity scales as
$O(MN^{2}\log N)$ where in practice $M\ll N$, but in the worst case $M$ can be
of the same order as $N$. On the other hand, the direct diagonalization scales
like $O(N^{3})$. It is not a priori clear that the COCG approach is preferred.
It is at this point that leveraging the strong scalability feature of the COCG
approach wins. Let us use a distinct parallel process for each point in the
reference system reducing the calculation to the cost $O(MN\log N)$. Thus,
within the errors of the converged results we present, the COCG approach gains
a clear advantage for larger dimensional systems.

We have also investigated other algorithms that can be applied to non-
symmetric shifted linear systems, the shifted BiCG-Stab($l$) Frommer (2003)
and the shifted GMRES Frommer and Glässner (1998). In these algorithms, we do
not need to play the tricks described in Sec. III.1 for the Hermitian matrix
$H$ and the size of matrix $A$ will not be increased by a factor of two. But,
these algorithms require two sequential matrix-vector (MV) products in each
iteration; thus in total there is no profit per iteration. Among the tested
methods, the COCG method exhibits so far the best convergence properties.

Finally, we emphasize that the presented method can efficiently utilize
heterogeneous computers. In our implementation we perform calculations both
using CPUs as well as highly efficient multi-threaded GPUs. In our experience,
the GPU implementation of shifted COCG is more than 50$\times$ faster than its
CPU counterpart. Table 1 compares the charged core hours between our GPU code
implementing the shifted COCG method and the CPU code using the direct
diagonalization method (ScaLAPACK) for two problems with different lattice
dimensions. On supercomputers, this is the most relevant quantity to be
compared as it essentially represents the cost that a user has to pay for
calculations. For shifted COCG method, the axial symmetry described in Sec.
III.4 is implemented. These timing tests were performed at OLCF Titan tit  and
NERSC Edison edi  supercomputers. In both problems, shifted COCG is faster
than diagonalization, and this advantage becomes more pronounced in the
problem with larger dimension (no.2).

No. | $N$ |  $t_{\mathrm{cocg}}$ (CPU hrs) |  $t_{\mathrm{diag}}$ (CPU hrs)
---|---|---|---
1 | $20\times 20\times 40\times 4$ | 649.23 (T) | 1,547.71 (T)
2 | $50\times 50\times 40\times 2$ | 9,318.4 (T) | 46,694.4 (E)
Table 1:  Timing comparison between the shifted COCG method and
diagonalization method using ScaLAPACK for solving two problems with different
dimensions $N$s in each self-consistent iteration. The first case is for
finite nuclei 240Pu in the lattice of dimension $20\times 20\times 40$ (the
extra factor of 4 is included for spin-orbit and pairing). The second case is
for the study of a nucleus immersed in a neutron superfluid sea of dimension
$50\times 50\times 40$, where spin-orbit is ignored (the extra factor of 2 is
for pairing only). In both problems, $t_{\mathrm{cocg}}$ and
$t_{\mathrm{diag}}$ denote the charged core hours of the GPU code using
shifted COCG method and the CPU code using ScaLAPACK in one self-consistent
iteration. The capital letter in parentheses denotes the computing facility
where the timing test is performed: “T” represents Titan and “E” represents
Edison tit ; edi  and the CPU hours were determined in units of CPU hours
according to the corresponding policy for these supercomputers. On Titan an
hour using a single computing node, which has 16 CPUs and one GPU, is charged
as 30 CPU h. On Edison an hour using a CPU is charged as 2 CPU h. For problem
no.2 the estimated charge in case of using Titan for direct diagonalization is
about 87 500 CPU hrs.

##  IV Benchmark examples

In this section, the Green’s function and the shifted COCG method, which we
denote XCOCG, is benchmarked by solving the self-consistent HFB equation in 3D
coordinate space and compared with the codes used in Refs. Stetcu  _et al._
(2011, 2015); Bulgac  _et al._ (2016), which we denote XDIAG. The XDIAG code
extracts wave functions of HFB Hamiltonian (15) in the discrete variable
representation (DVR) basis Bulgac and Forbes (2013) via a direct
diagonalization. 111Direct methods yield the exact solution if the precision
is unlimited whereas indirect methods may be even more accurate and are often
more efficient depending on system details. Also, iterative methods more tend
to naturally damp out roundoff errors that accumulate and become difficult
with large $N$ problems evaluated directly. The diagonalization procedure is
executed parallel using ScaLAPACK library sca . Next, densities are formed
from the wave-functions using formulas (18-22). These operations form a single
self-consistent iteration. We mix the intermediate solutions during the
iteration process using a linear or a Broyden mixing algorithm Baran  _et al._
(2008). The XCOCG code is a modified version of the XDIAG code where parallel
diagonalization procedure and computation of densities is replaced by the COCG
method. Densities are extracted according formulas provided in Sec. II.2. The
COCG part is ported to GPUs. Both codes provide the same results up to the
accuracy specified in Sec. III.3.

We performed tests on the nuclei with axial symmetry, including spherical and
axially deformed nuclei. In all of these realistic tests, the Skyrme NEDF SLy4
Chabanat  _et al._ (1998) is used in p-h channel, and the SLDA treatment
Bulgac and Yu (2002); Yu and Bulgac (2003) for the pairing interaction is used
in p-p channel with bare coupling strength
$g_{0}(\boldsymbol{r})=g_{0}=-233~{}\mathrm{MeV}$ (volume pairing), and energy
cutoff $E_{\mathrm{cut}}=100~{}\mathrm{MeV}$. These are typical values used in
realistic calculations. We conclude this section by presenting states
generated by the shifted COCG method used for studies of interaction between
quantized vortices and nuclei in neutron start crust Wlazłowski  _et al._
(2016). These calculations are not feasible for XDIAG code in a reasonable
time.

| $dx=1.25~{}\mathrm{fm}$ | $dx=1.0~{}\mathrm{fm}$ | $dx=0.8~{}\mathrm{fm}$ | HFBRAD
---|---|---|---|---
$E_{\mathrm{kin}}$ | 3868.984 | 3866.098 | 3866.187 | 3866.163
$E_{\rho^{2}}$ | -22401.965 | -22383.542 | -22384.385 | -22384.462
$E_{\rho^{\gamma}}$ | 14548.826 | 14536.219 | 14536.840 | 14536.890
$E_{\rho\Delta\rho}$ | 315.563 | 315.236 | 315.286 | 315.288
$E_{\rho\tau}$ | 1332.134 | 1330.147 | 1330.208 | 1330.216
$E_{\rho\nabla J}$ | -96.592 | -96.451 | -96.446 | -96.446
$E_{\mathrm{Coul}}$ | 796.848 | 796.600 | 796.607 | 796.645
$E_{\mathrm{tot}}$ | -1636.202 | -1635.693 | -1635.703 | -1635.707
Table 2:  Results of HFB + SLy4 calculations for 208Pb using XCOCG in a cubic box of size $32^{3}~{}\mathrm{fm}^{3}$ with lattice constants $dx=1.25,1.0,0.8~{}\mathrm{fm}$. Results obtained with the spherical 1D code HFBRAD are presented for comparison. All energies are in MeV.  | $dx=1.25~{}\mathrm{fm}$ | $dx=1.0~{}\mathrm{fm}$ | $dx=0.8333~{}\mathrm{fm}$
---|---|---|---
| XDIAG | XCOCG | $\Delta E$ | XDIAG | XCOCG | $\Delta E$ | XDIAG | XCOCG | $\Delta E$
$E_{\mathrm{kin}}$ | 623.983 | 624.115 | 0.132 | 624.755 | 624.781 | 0.026 | 624.836 | 624.837 | <0.001
$E_{\rho^{2}}$ | -3714.569 | -3714.702 | -0.133 | -3722.988 | -3722.977 | 0.011 | -3723.483 | -3723.483 | <0.001
$E_{\rho^{\gamma}}$ | 2396.941 | 2396.967 | 0.026 | 2403.198 | 2403.190 | -0.008 | 2403.555 | 2403.554 | <0.001
$E_{\rho\Delta\rho}$ | 106.284 | 106.368 | 0.084 | 106.957 | 106.957 | <0.001 | 106.992 | 106.992 | <0.001
$E_{\rho\tau}$ | 173.227 | 173.122 | -0.105 | 173.725 | 173.712 | -0.013 | 173.746 | 173.746 | <0.001
$E_{\rho\nabla J}$ | -1.233 | -1.267 | -0.033 | -1.279 | -1.287 | -0.008 | -1.282 | -1.282 | <0.001
$E_{\mathrm{Coul}}$ | 71.483 | 71.484 | <0.001 | 71.532 | 71.532 | <0.001 | 71.535 | 71.535 | <0.001
$E_{\mathrm{tot}}$ | -343.885 | -343.914 | -0.029 | -344.100 | -344.093 | 0.007 | -344.101 | -344.101 | <0.001
Table 3:  Results of spherical HFB + SLy4 calculations for 40Ca using the XCOCG and XDIAG approaches with a cubic 3D lattice of size $L_{x}=20~{}\mathrm{fm}$ with different mesh size $dx$. All energies are in MeV. In the column of $\Delta E$, “$<0.001$” means $\lvert\Delta E\lvert$ is less than 1 keV and can be negligible.  | $dx=1.25~{}\mathrm{fm}$ | $dx=1.0~{}\mathrm{fm}$ | $dx=0.8333~{}\mathrm{fm}$
---|---|---|---
| XDIAG | XCOCG | $\Delta E$ | XDIAG | XCOCG | $\Delta E$ | XDIAG | XCOCG | $\Delta E$
$E_{\mathrm{kin}}$ | 1092.413 | 1092.737 | 0.324 | 1090.235 | 1090.293 | 0.058 | 1090.107 | 1090.122 | 0.015
$E_{\rho^{2}}$ | -6335.819 | -6336.367 | -0.547 | -6319.916 | -6320.133 | -0.218 | -6318.595 | -6318.650 | -0.055
$E_{\rho^{\gamma}}$ | 4134.013 | 4134.153 | 0.140 | 4122.721 | 4122.864 | 0.143 | 4121.738 | 4121.775 | 0.037
$E_{\rho\Delta\rho}$ | 150.018 | 150.256 | 0.238 | 150.031 | 150.055 | 0.024 | 149.967 | 149.971 | 0.004
$E_{\rho\tau}$ | 341.670 | 341.492 | -0.178 | 339.456 | 339.460 | 0.005 | 339.326 | 339.329 | 0.003
$E_{\rho\nabla J}$ | -52.747 | -52.819 | -0.072 | -52.371 | -52.380 | -0.010 | -52.367 | -52.365 | 0.001
$E_{\mathrm{Coul}}$ | 129.928 | 129.917 | -0.011 | 129.761 | 129.762 | <0.001 | 129.751 | 129.751 | <0.001
$E_{\mathrm{pair}}$ | -3.834 | -3.904 | -0.070 | -3.989 | -4.000 | -0.011 | -3.986 | -3.994 | -0.009
$E_{\mathrm{tot}}$ | -544.358 | -544.534 | -0.175 | -544.073 | -544.080 | -0.007 | -544.058 | -544.061 | -0.003
Table 4:  Same as Table 3, but for 62Ni.  | $dx=1.25~{}\mathrm{fm}$ | $dx=0.9375~{}\mathrm{fm}$
---|---|---
| XDIAG | XCOCG | $\Delta E(\Delta Q)$ | XDIAG | XCOCG | $\Delta E(\Delta Q)$
$E_{\mathrm{kin}}$ | 1841.034 | 1841.444 | 0.410 | 1838.435 | 1838.438 | 0.003
$E_{\rho^{2}}$ | -10381.846 | -10383.607 | -1.760 | -10364.579 | -10364.463 | 0.115
$E_{\rho^{\gamma}}$ | 6715.556 | 6716.620 | 1.064 | 6702.965 | 6702.865 | -0.100
$E_{\rho\Delta\rho}$ | 209.650 | 210.055 | 0.405 | 208.432 | 208.427 | -0.005
$E_{\rho\tau}$ | 592.475 | 592.366 | -0.109 | 590.686 | 590.669 | -0.017
$E_{\rho\nabla J}$ | -63.446 | -63.716 | -0.270 | -62.429 | -62.422 | 0.008
$E_{\mathrm{Coul}}$ | 230.311 | 230.288 | -0.023 | 230.359 | 230.358 | -0.001
$E_{\mathrm{pair}}$ | -3.067 | -3.083 | -0.016 | -3.231 | -3.233 | -0.002
$E_{\mathrm{tot}}$ | -859.333 | -859.632 | -0.299 | -859.361 | -859.360 | 0.001
$Q_{20}$ | 1077.44 | 1080.62 | 3.18 | 1047.81 | 1047.65 | -0.16
Table 5:  Results of HFB + SLy4 calculations for 102Zr using XCOCG and XDIAG
in a rectangular box of size $22.5\times 22.5\times 30~{}\mathrm{fm}^{3}$ with
lattice constants $dx=1.25$ and $0.9375~{}\mathrm{fm}$ . All energies are in
MeV and the quadrupole moment $Q_{20}$ is in $\mathrm{fm}^{2}$.

###  IV.1 Spherical limit: Doubly magic nuclei 208Pb, 40Ca and semi magic
nucleus 62Ni

A common approach to testing the accuracy of a 3D coordinate solver is to
compare the 3D results for a spherical nucleus with the results obtained by a
1D spherical code that also represents the single-particle wave functions in
coordinate space. We can choose extremely fine lattice constants for the 1D
solver and its results can be considered to be very accurate. Instead of a
benchmark with XDIAG, we compute the double magic nucleus 208Pb using our
XCOCG code in a cubic box of size $32^{3}~{}\mathrm{fm}^{3}$ with lattice
constants $dx=1.25,1.0,0.8~{}\mathrm{fm}$, and compare them with the results
obtained with the HFBRAD code Bennaceur and Dobaczewski (2005) code (the
lattice constant $dx=0.05~{}\mathrm{fm}$). In Table 2 we compare various
contributions to the total energy $E_{\mathrm{tot}}$, computed as volume
integrals from corresponding terms in NEDF for 208Pb in these situations. From
$dx=1.25$ to $0.8~{}\mathrm{fm}$, the difference of total energy between XCOCG
and HFBRAD decreases from 0.5 MeV to 4 keV. In particular, from $dx=1.0$ to
$0.8~{}\mathrm{fm}$, the values of the energy terms have a steady convergence
to those solved by HFBRAD with the maximum difference $\leq 100~{}$keV. This
numerical accuracy and convergence pattern is similar to the results in Ref.
Ryssens  _et al._ (2015b), which uses Lagrange-mesh representation Baye (2015)
in the calculation of spatial derivatives. The Lagrange-mesh method is
equivalent to the DVR method using FFT on a 3D spatial lattice.

The benchmarking with XDIAG starts with the case of spherical nuclei. We solve
the HFB problem for the doubly magic nucleus ${}^{40}\mathrm{Ca}$ and semi-
magic nucleus ${}^{62}\mathrm{Ni}$ to demonstrate the accuracy of our solver
in the cases with and without pairing. In both cases, we use a cubic lattice
of size $20^{3}~{}\mathrm{fm}^{3}$ with different lattice constants
$\mathrm{d}x=1.25,1.0,0.8333~{}\mathrm{fm}$. In Table 3 we compare various
contributions to the total energy of 40Ca calculated by XCOCG and XDIAG
methods respectively. The main source of differences between the results of
the two solvers is due to the neglect of the azimuthal component of the spin-
orbit density $\boldsymbol{J}_{\phi}(\boldsymbol{r})$ when we utilize the
axial symmetry of the system. This lattice effect vanishes as the lattice
constant decreases. In particular, when $dx\leq 1~{}\mathrm{fm}$, the
difference between the total energies is less than 10 keV. Notice also that
the XCOCG energies are always lower than the XDIAG energies and we attribute
this to the fact that $\boldsymbol{J}_{\phi}(\boldsymbol{r})\equiv 0$ in
XCOCG. The same kind of calculations are performed for the semi-magic nucleus
62Ni with non-zero neutron pairing; see Table 4. Compared to 40Ca, the
difference of total energy is larger in $dx=1.25$ case ($0.175$ MeV). This is
due to the much larger magnitude of the spin-orbit contribution, which brings
larger error in $\boldsymbol{J}$ caused by the finite lattice effects.
Similarly to 40Ca, as the lattice constant becomes finer, this lattice effect
vanishes and the difference of total energy drops to values below $10$ keV
when $dx\leq 1~{}\mathrm{fm}$.

###  IV.2 Axially deformed nucleus: 102Zr

The advantage of solving the HFB equation in a coordinate space basis is that
it can correctly describe the asymptotic behavior of quasiparticle
wavefunctions of nuclei with large deformations and weak binding energies
Dobaczewski  _et al._ (1996). Following Refs. Pei  _et al._ (2008); Terán  _et
al._ (2003); Blazkiewicz  _et al._ (2005), we choose the neutron rich Zr
isotope 102Zr, which has a large prolate deformation, as the testing ground
and we choose a rectangular box of size $22.5\times 22.5\times
30~{}\mathrm{fm}^{3}$ to fit its large deformation. As in Sec. IV.1, we
compare the XCOCG and XDIAG results with different lattice sizes $dx=1.25$ and
$0.9375~{}\mathrm{fm}$. These are shown in Table. 5. The quadrupole moment
$Q_{20}$ of the nucleus is also listed in each case, where

$\displaystyle
Q_{20}=\langle\hat{Q}\rangle=\int(2z^{2}-x^{2}-y^{2})\rho(\boldsymbol{r})\,d^{3}\boldsymbol{r}$
(59)

and $\rho(\boldsymbol{r})=\rho_{n}(\boldsymbol{r})+\rho_{p}(\boldsymbol{r})$.

The difference in the total energy is 0.3 MeV between XCOCG and XDIAG for
$dx=1.25$, which almost equals to the difference in the spin-orbit energy 0.27
MeV. This error in the spin-orbit energy contribution will affect the position
of the HFB minimum, leading to a difference in the quadrupole moment $Q_{20}$
in the ground state of $3~{}\mathrm{fm}^{2}$ and in the energy terms
$E_{\rho^{2}},E_{\rho^{\gamma}}$ of 1–2 MeV between XCOCG and XDIAG. For
$dx=0.9375\;\textrm{fm}$, the errors are significantly smaller, $\leq
1~{}\mathrm{keV}$ for the total energy and $0.16~{}\mathrm{fm}^{2}$ for the
quadrupole moment.

###  IV.3 Constrained HFB: saddle-point of 240Pu fission

Induced fission of 240Pu is a frequent benchmark for many implementations
based on DFT methods Younes and Gogny (2009); Schunck  _et al._ (2014); Bulgac
_et al._ (2016); Regnier  _et al._ (2016). In a constrained HFB calculation in
Ref. Bulgac  _et al._ (2016) the nucleus is brought to a shape and an energy
near the outer saddle point of the fission barrier (at zero temperature), used
as the initial state of a time-dependent SLDA (TDSLDA) simulation. In this
work, we use our XCOCG code to reproduce the configuration on this saddle
point and compare with that obtained via the XDIAG code.

| XDIAG | XCOCG | $\Delta E$
---|---|---|---
$E_{\mathrm{kin}}$ | 4418.132 | 4419.125 | 0.993
$E_{\rho^{2}}$ | -25104.582 | -25109.048 | -4.467
$E_{\rho^{\gamma}}$ | 16238.995 | 16241.781 | 2.786
$E_{\rho\Delta\rho}$ | 408.182 | 409.032 | 0.850
$E_{\rho\tau}$ | 1465.062 | 1464.900 | -0.161
$E_{\rho\nabla J}$ | -111.516 | -112.088 | -0.572
$E_{\mathrm{Coul}}$ | 901.083 | 901.062 | -0.021
$E_{\mathrm{pair}}$ | -8.801 | -8.818 | -0.017
$E_{\mathrm{tot}}$ | -1793.439 | -1794.054 | -0.615
Table 6:  Results of HFB + SLy4 calculations for the saddle point of 240Pu
using XCOCG and XDIAG in a rectangular box of size $25\times 25\times
50~{}\mathrm{fm}^{3}$ with mesh size $dx=dy=dz=1.25~{}\mathrm{fm}$. All
energies are in MeV.

To obtain a nucleus with a given quadrupole moment
$\langle\hat{Q}\rangle=Q_{0}$, we need to minimize the Routhian

$\displaystyle E^{\prime}=E_{\mathrm{total}}+c(\langle\hat{Q}\rangle-
Q_{0})^{2},$ (60)

which is equivalent to adding a Lagrange multiplier in the single-particle
Hamiltonian: $h^{\prime}=h+2c(\langle\hat{Q}\rangle-Q_{0})\hat{Q}$. In the
self-consistent calculation, the constraint strength $c$ will be updated in
each iteration using the augmented Lagrangian method Staszczak  _et al._
(2010). This saddle-point configuration of 240Pu is prepared in a rectangular
box of size $25\times 25\times 50~{}\mathrm{fm}^{3}$, with lattice spacing
$dx=1.25~{}\mathrm{fm}$. Following paper Bulgac  _et al._ (2016), the
quadrupole moment constraint is set to $Q_{0}=16500~{}\mathrm{fm}^{2}$ and an
additional auxiliary external field is turned on for the formation of the neck
in this configuration.

In Fig. 4 we show the density profile for the converged solution.

Figure 4:  Total density distribution
$\rho(\boldsymbol{r})=\rho_{n}(\boldsymbol{r})+\rho_{p}(\boldsymbol{r})$ for
the saddle-point configuration of 240Pu in the $y=0$ plane.

Table 6 compares the XCOCG and XDIAG results for the saddle-point
configuration of 240Pu. As we discussed earlier, the lattice effects in the
calculation of $\boldsymbol{J}$ cause the change of the HFB minimum, which
results in large differences in the energy terms $E_{\mathrm{kin}}$,
$E_{\rho^{2}}$, $E_{\rho^{\gamma}}$, and $E_{\rho\Delta\rho}$. In this case,
because of the strong pairing, and extremely large deformation of the heavy
nucleus, any minute change of HFB minimum will result in larger difference in
various energy terms than in the cases of the nuclei we tested earlier. On the
other hand, the larger amplitude of $\boldsymbol{J}$, which is reflected by
the absolute value of $E_{\rho\nabla J}$, also makes the lattice effects more
pronounced. In principle, these differences will reduce and eventually vanish
as $dx\to 0$ as we demonstrated in the earlier sections. But the large
dimension of this system and the considerable number of self-consistent
iterations in the constrained HFB problem will make the calculation on a finer
lattice grid more expensive with the current computational resources.

###  IV.4 A nucleus immersed in a neutron superfluid sea

We used the shifted COCG method to generate initial states for studies of
quantum vortex dynamics in a neutron star crust Wlazłowski  _et al._ (2016).
The stationary state is a superfluid neutron medium containing a quantum
vortex and an immersed nucleus located in the vicinity of the topological
defect or vortex. The calculations were performed with the FaNDF0 nuclear
density functional constructed by Fayans et al. Fayans (1998); Fayans  _et
al._ (2000), which is particularly well suited for these type of studies. The
bare pairing coupling constant in Eq. (10) was chosen to reproduce the
rescaled BCS ${}^{1}S_{0}$ pairing gap in neutron matter; for more details,
see Ref. Wlazłowski  _et al._ (2016).

For these studies we used a simulation volume of size $75\times 75\times
60\;\textrm{fm}^{3}$ with a lattice spacing $dx=1.5\;\textrm{fm}$. The energy
cut-off was chosen to be $75\;\textrm{MeV}$. In the box, we place a tube
(simulated by a flat-bottomed external potential) that we fill with superfluid
neutrons of density $n=0.014\;\textrm{fm}^{-3}$ or $0.031\;\textrm{fm}^{-3}$.
The problem has been simplified by dropping the spin-orbit term, which is not
expected to play a major role in vortex pinning. The simplification results in
an HFB matrix (15) with a simpler block structure
($h_{\uparrow\downarrow}=h_{\downarrow\uparrow}=0$) and a smaller dimension
$2N_{x}N_{y}N_{z}$:

$\displaystyle\begin{pmatrix}h_{\uparrow\uparrow}-\mu&\Delta\\\
\Delta^{*}&-h^{*}_{\downarrow\downarrow}+\mu\end{pmatrix}\begin{pmatrix}u_{k\uparrow}\\\
v_{k\downarrow}\end{pmatrix}=E_{k}\begin{pmatrix}u_{k\uparrow}\\\
v_{k\downarrow}\end{pmatrix}.$ (61)

After this reduction, the HFB matrix has the size $200,000\times 200,000$ and
it still represents a very demanding problem for the traditional approaches to
determine the stationary states, see Section III.5 and Table 1. We solved this
problem successfully with moderate computational costs to achieve self-
consistency using GPUs and the COCG approach described here on the Titan
supercomputer tit . Only after the final iteration, do we use a
diagonalization to generate the wave functions on the Edison supercomputer edi
.

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

In Fig. 5 we present stationary configurations (with constraints) for two
background neutron densities: $n=0.014\;\textrm{fm}^{-3}$ and
$0.031\;\textrm{fm}^{-3}$. The “nuclear defect” consists of $Z=50$ protons.
Two mutual configurations were considered: (1) a quantum vortex attached to
the nucleus (pinned case); and (2) a nucleus outside the vortex core (unpinned
case). The position of nucleus in the box was fixed by adding a constraint to
the density functional for the center of mass of the protons, in a similar
fashion as it was done for the quadrupole moment Eq. (60) and for the for
saddle point in 240Pu test case. The vortex was generated by imprinting the
correct phase pattern for the pairing field $\Delta$ in the neutron channel
$\Delta(\rho,z,\phi)=|\Delta(\rho,z)|\exp(i\phi)$, where
$\rho=\sqrt{x^{2}+y^{2}}$ is the distance from the center of the tube and
$\phi=\tan^{-1}\tfrac{y}{x}$. In the pinned case, we took advantage of the
axial symmetry of the problem and of the reflection symmetry with respect to
one plane in the unpinned configuration.

From the energetics of these systems we determine that for both densities the
configuration with the nucleus located outside the vortex core has a lower
energy per particle than the pinned configuration by about $6$ and
$4\;\textrm{keV}$ respectively for densities $n=0.014\;\textrm{fm}^{-3}$ and
$0.031\;\textrm{fm}^{-3}$. Thus one can expect that the effective vortex-
nucleus interaction is repulsive in nature, which was further confirmed by
studying the motion of the vortex in dynamical simulations Wlazłowski  _et
al._ (2016).

##  V Further extension

###  V.1 Linear response

The free linear response of a many-fermion system can be evaluated using the
same COCG approach. For simplicity we will illustrate this procedure here only
for a normal system, as the extension to the superfluid case is
straightforward. The free polarization operator is defined as Shlomo and
Bertsch (1975); Bertsch and Tsai (1975)

$\displaystyle\begin{split}\Pi_{0}(\omega,\boldsymbol{r},\boldsymbol{r}^{\prime})=\sum_{k}\psi_{k}(\boldsymbol{r})\psi^{*}_{k}(\boldsymbol{r}^{\prime})\\\
\times\left[G(\varepsilon_{k}+\omega,\boldsymbol{r},\boldsymbol{r}^{\prime})+G(\varepsilon_{k}-\omega,\boldsymbol{r},\boldsymbol{r}^{\prime})\right],\end{split}$
(62)

where $\psi_{k}(\boldsymbol{r})$ and
$G(\boldsymbol{r},\boldsymbol{r}^{\prime},z)$ were defined in Eqs. (2) and
(4), and the summation is over occupied levels. In practice the polarization
operator is evaluated for a complex energy $\omega+i\gamma$, since the
spreading width $\Gamma^{\downarrow}$ is either not accounted for in the
random phase approximation, or in order to imitate to some extent the
imaginary part of the optical potential. The spreading width accounts for the
fragmentation of the particle-hole transition strength due to the coupling to
more complex states. If the complex integration contour in Eq. (6) is chosen
so that $|\mathrm{Im}\;z|<\gamma$ one can easily show that

$\displaystyle\begin{split}\Pi_{0}(\omega+i\gamma,\boldsymbol{r},\boldsymbol{r}^{\prime})=\frac{1}{2\pi
i}\oint_{C}G(z,\boldsymbol{r},\boldsymbol{r}^{\prime})\\\
\times\left[G(z+\omega+i\gamma,\boldsymbol{r},\boldsymbol{r}^{\prime})+G(z-\omega-i\gamma,\boldsymbol{r},\boldsymbol{r}^{\prime})\right].\end{split}$
(63)

It remains to be established if this method for the evaluation of the free
polarization operator
$\Pi_{0}(\omega+i\gamma,\boldsymbol{r},\boldsymbol{r}^{\prime})$ and the
subsequent determination of the full response is competitive and under what
conditions with the finite-amplitude method Nakatsukasa  _et al._ (2007);
Avogadro and Nakatsukasa (2011); Nakatsukasa (2014); Hinohara  _et al._
(2013); Oishi  _et al._ (2016); Mustonen  _et al._ (2014) and/or the time-
dependent approach Stetcu  _et al._ (2011, 2015) for deform open-shell nuclei.

###  V.2 Extraction of eigenvalues and eigenvectors using the shifted COCG
Krylov method

The shifted COCG Krylov method can be used also to determine the eigenvalues
and the eigenvectors using the approach described in Refs. Sakurai and Sugiura
(2003); Sakurai  _et al._ (2008); Ikegami  _et al._ (2010). In this approach
one has at first to evaluate the moments

$\displaystyle\mu_{k}=\frac{1}{2\pi
i}\oint_{\Gamma}dz\sum_{n}\frac{(z-\epsilon)^{k}}{z-\varepsilon_{n}},\;\;\;k=0,1,...,2N-1.$
(64)

where $\epsilon$ is located inside the contour $\Gamma$ enclosing a segment on
the real axis with a known number $N$ of the eigenvalues, and where
$\varepsilon_{n}$ are eigenvalues. Once these moments are computed the
eigenvalues are obtained by solving a generalized eigenvalue problem for two
matrices of size $N\times N$. The number of eigenvalues in a given interval is
not known a priori and some eigenvalues could also be degenerate. In the
presence of degeneracies one has to disentangle the corresponding
eigenvectors. If the degeneracy is due to spherical or axial symmetry one can
introduce slightly different lattice constants $dx$, $dy$, and $dz$ or a very
weak external field and lift the degeneracies at a level that has no
noticeable consequence on the physics studied. Subsequently, the solution of
the Schrödinger equation for the corresponding eigenvector $\psi_{n}$ in the
case of a known non-degenerate eigenvalue $(H-\varepsilon_{n})\psi_{n}=0$ is a
trivial linear algebra problem. In the case of Kramers degeneracies, one can
easily separate the two degenerate eigenvectors. The unknown number $N$ of
eigenvalues in a given energy interval can be determined by evaluating the
trace (integral over all coordinates and summation over all four components)
of the Green’s function of the Hamiltonian (15)

$\displaystyle N=\mathrm{Tr}\left(\frac{1}{2\pi
i}\oint_{\Gamma}dz\frac{1}{z-H}\right).$ (65)

Thus the need to use diagonalization of very large matrices can be completely
eschewed.

##  VI Conclusions

In this paper, we describe a new approach for solving the HFB type of
equations in a coordinate representation that is different from the
traditional approaches based typically on direct diagonalizations. In the
present approach, there is no need to evaluate the individual single-particle
wave functions and their energies, but instead we calculate the Green’s
function of the HFB equations, from which we extract various densities after
evaluating a contour integral. The Green’s functions are obtained by solving a
set of linear equations with scalar shifts using the iterative shifted COCG
Krylov method. We demonstrate the high accuracy of the iterative shifted COCG
approach by solving typical nuclear problems with and without complicated
constraints, such as the fission outer saddle-point of 240Pu and a quantum
vortex state in a neutron star crust. A notable advantage of this algorithm is
its suitability for efficient parallelization and effective utilization of
heterogeneous computing platforms. The method becomes computationally superior
for large spatial lattice sizes that are otherwise computationally very
expensive for standard approaches, such as a direct diagonalization.

###### Acknowledgements.

We thank G.F. Bertsch for discussions and I. Stetcu and P. Magierski for many
helpful discussions and code sharing. This work was supported in part by U.S.
DOE Office of Science Grant DE-FG02-97ER41014. This work was also supported in
part by the Polish National Science Center (NCN) under Contract No.
UMO-2014/13/D/ST3/01940. Calculations have been performed: at the OLCF Titan -
resources of the Oak Ridge Leadership Computing Facility, which is a DOE
Office of Science User Facility supported under Contract DE-AC05-00OR22725; at
NERSC Edison - resources of the National Energy Research Scientific computing
Center, which is supported by the Office of Science of the U.S. Department of
Energy under Contract No. DE-AC02-05CH11231; and at Moonlight, resources of
the Institutional Computing Program at Los Alamos National Laboratory.

## APPENDIX: Equation of Green’s function for HFB equation

For the HFB Hamiltonian $H$ in Eq. (15), the single-particle Hamiltonian in
Eq. (16) reads

$\displaystyle
h=-\boldsymbol{\nabla}\cdot\frac{\hbar^{2}}{2m^{*}(\boldsymbol{r})}\boldsymbol{\nabla}+U(\boldsymbol{r})-i\boldsymbol{W}(\boldsymbol{r})\cdot\left(\boldsymbol{\nabla}\times\boldsymbol{\sigma}\right),$
(66)

where the kinetic energy term is represented by a real symmetric operator in
numerical implementations $T$. When discretized, besides the local term
$U(\boldsymbol{r})_{mn}=U_{n}\delta_{nm}$, the nonlocal term $T$ and spin-
orbit terms require more attention. For $T$, since

$\displaystyle\begin{split}&-\boldsymbol{\nabla}\cdot\frac{\hbar^{2}}{2m^{*}(\boldsymbol{r})}\boldsymbol{\nabla}v(\boldsymbol{r})=-\frac{1}{2}\left[\frac{\hbar^{2}}{2m^{*}(\boldsymbol{r})}\nabla^{2}v(\boldsymbol{r})\right.\\\
&\left.+\nabla^{2}\left(\frac{\hbar^{2}}{2m^{*}(\boldsymbol{r})}v(\boldsymbol{r})\right)-\left(\nabla^{2}\frac{\hbar^{2}}{2m^{*}(\boldsymbol{r})}\right)v(\boldsymbol{r})\right],\end{split}$
(67)

we have

$\displaystyle\begin{split}&T_{nm}=\left(-\boldsymbol{\nabla}\cdot\frac{\hbar^{2}}{2m^{*}(\boldsymbol{r})}\boldsymbol{\nabla}\right)_{nm}=\\\
&-\frac{1}{2}(\nabla^{2})_{nm}\left(\frac{\hbar^{2}}{2m^{*}_{n}}+\frac{\hbar^{2}}{2m^{*}_{m}}\right)+\frac{1}{2}\left(\nabla^{2}\frac{\hbar^{2}}{2m^{*}}\right)_{n}\delta_{nm},\end{split}$
(68)

where the Laplacian operator $\nabla^{2}$ is a symmetric matrix in the
discrete variable representation (DVR) basis Bulgac and Forbes (2013) and
$m^{*}_{n/m}=m^{*}(\boldsymbol{r}_{n/m})$. A similar symmetrization is also
performed for the spin-orbit term

$\displaystyle\begin{split}&\boldsymbol{W}(\boldsymbol{r})\cdot\left(\boldsymbol{\nabla}\times\boldsymbol{\sigma}\right)\\\
&=\frac{1}{2}\left[\boldsymbol{W}(\boldsymbol{r})\cdot\left(\boldsymbol{\nabla}\times\boldsymbol{\sigma}\right)+\boldsymbol{\nabla}\cdot\left(\boldsymbol{\sigma}\times\boldsymbol{W}(\boldsymbol{r})\right)\right],\end{split}$
(69)

where the gradient operator
$\boldsymbol{\nabla}=\\{\partial_{x},\partial_{y},\partial_{z}\\}$ is an anti-
symmetric operator in each spatial direction in the DVR representation. In the
p-p channel, the complex pairing field $\Delta(\boldsymbol{r})$ is diagonal
when discretized, $\Delta(\boldsymbol{r})_{mn}=\Delta_{n}\delta_{nm}$.
Finally, after separating the real and imaginary parts of the HFB Hamiltonian
$H=A+iB$, the equation of Green’s function $G=G_{x}+iG_{y}$ for a HFB equation
is

$\displaystyle(zI-H^{\prime})\begin{pmatrix}G_{x}\\\
G_{y}\end{pmatrix}=\begin{pmatrix}\delta^{(4)}\\\ 0\end{pmatrix},\quad
H^{\prime}=\begin{pmatrix}A&-B\\\ B&A\end{pmatrix},$ (70)

and where

$\displaystyle A$
$\displaystyle=\begin{pmatrix}T+\tilde{U}&\widetilde{W_{x}\partial_{z}}-\widetilde{W_{z}\partial_{x}}&0&\mathrm{Re}~{}\Delta\\\
\widetilde{W_{z}\partial_{x}}-\widetilde{W_{x}\partial_{z}}&T+\tilde{U}&-\mathrm{Re}~{}\Delta&0\\\
0&-\mathrm{Re}~{}\Delta&-T-\tilde{U}&\widetilde{W_{z}\partial_{x}}-\widetilde{W_{x}\partial_{z}}\\\
\mathrm{Re}~{}\Delta&0&\widetilde{W_{x}\partial_{z}}-\widetilde{W_{z}\partial_{x}}&-T-\widetilde{U}\end{pmatrix},$
(71) $\displaystyle B$
$\displaystyle=\begin{pmatrix}\widetilde{W_{y}\partial_{x}}-\widetilde{W_{x}\partial_{y}}&\widetilde{W_{z}\partial_{y}}-\widetilde{W_{y}\partial_{z}}&0&\mathrm{Im}~{}\Delta\\\
\widetilde{W_{z}\partial_{y}}-\widetilde{W_{y}\partial_{z}}&\widetilde{W_{x}\partial_{y}}-\widetilde{W_{y}\partial_{x}}&-\mathrm{Im}~{}\Delta&0\\\
0&\mathrm{Im}~{}\Delta&\widetilde{W_{y}\partial_{x}}-\widetilde{W_{x}\partial_{y}}&\widetilde{W_{z}\partial_{y}}-\widetilde{W_{y}\partial_{z}}\\\
-\mathrm{Im}~{}\Delta&0&\widetilde{W_{z}\partial_{y}}-\widetilde{W_{y}\partial_{z}}&\widetilde{W_{x}\partial_{y}}-\widetilde{W_{y}\partial_{x}}\end{pmatrix},$
(72)

with $\tilde{U}=U-\mu$ and $\widetilde{W_{i}\partial_{j}}$ the anti-commutator
of $W_{i}$ and $\partial_{j}$:

$\displaystyle\widetilde{W_{i}\partial_{j}}=\frac{1}{2}\left(W_{i}\partial_{j}+\partial_{j}W_{i}\right),\quad
i,j=x,y,z,i\neq j,$ (73)

One can show that $A=A^{T}$ and $B=-B^{T}$. In the shifted COCG method, we do
not need to construct these $4N_{xyz}\times 4N_{xyz}$ matrices explicitly
because we only need the MV product between these matrices and the vectors
$\boldsymbol{x}_{n}$, $\boldsymbol{p}_{n}$ (see Sec. III.1). Among the MV
product operation, the product between the local part of $H$ and the vectors
can be regarded as a vector-vector product. The product between the vectors
and the non-local part of $H$, due to the Laplacian and gradient operators,
can be performed using the fast Fourier transform.

## References

  * Bender  _et al._ (2003a) M. Bender, P.-H. Heenen, and P.-G. Reinhard, “Self-consistent mean-field models for nuclear structure,” Rev. Mod. Phys. 75, 121–180 (2003a).
  * Runge and Gross (1984) E. Runge and E. K. U. Gross, “Density-functional theory for time-dependent systems,” Phys. Rev. Lett. 52, 997–1000 (1984).
  * Marques  _et al._ (2006) M.A.L. Marques, C. A. Ullrich, F. Nogueira, A. Rubio, K. Burke, and E. U. K. Gross, eds.,  _Time-Dependent Density Functional Theory_, Lecture Notes in Physics, Vol. 706 (Springer, Berlin, 2006).
  * Marques  _et al._ (2012) A.L. Miguel Marques, T. Neepa Maitra, M.S. Fernando Nogueira, E.K.U. Gross, and Angel Rubio, eds.,  _Fundamentals of Time-Dependent Density Functional Theory_, Lecture Notes in Physics, Vol. 837 (Springer Berlin Heidelberg, 2012).
  * Ring and Schuck (2004) P. Ring and P. Schuck,  _The Nuclear Many-Body Problem_ (Springer, Berlin, 2004).
  * Bertsch (1980) G. Bertsch, “The nuclear density of states in the space of nuclear shapes,” Physics Letters B 95, 157 – 159 (1980).
  * Barranco  _et al._ (1990) F. Barranco, G.F. Bertsch, R.A. Broglia, and E. Vigezzi, “Large-amplitude motion in superfluid Fermi droplets,” Nuclear Physics A 512, 253 – 274 (1990).
  * Bertsch and Bulgac (1997) G. F. Bertsch and A. Bulgac, “Comment on “spontaneous fission: A kinetic approach”,” Phys. Rev. Lett. 79, 3539–3539 (1997).
  * Bulgac  _et al._ (2016) A. Bulgac, P. Magierski, K. J. Roche, and I. Stetcu, “Induced fission of ${}^{240}\mathrm{Pu}$ within a real-time microscopic framework,” Phys. Rev. Lett. 116, 122504 (2016).
  * Anderson and Itoh (1975) P. W. Anderson and N. Itoh, “Pulsar glitches and restlessness as a hard superfluid phenomenon,” Nature (London) 256, 25 (1975).
  * Wlazłowski  _et al._ (2016) G. Wlazłowski, K. Sekizawa, P. Magierski, A. Bulgac, and M. M. Forbes, “Vortex pinning and dynamics in the neutron star crust,” Phys. Rev. Lett. 117, 232701 (2016).
  * Davies  _et al._ (1980) K.T.R. Davies, H. Flocard, S. Krieger, and M.S. Weiss, “Application of the imaginary time step method to the solution of the static Hartree-Fock problem,” Nucl. Phys. A 342, 111 – 123 (1980).
  * Ryssens  _et al._ (2015a) W. Ryssens, V. Hellemans, M. Bender, and P.-H. Heenen, “Solution of the Skyrme–HF+BCS equation on a 3D mesh, II: A new version of the Ev8 code,” Computer Physics Communications 187, 175 – 194 (2015a).
  * Egido  _et al._ (1995) J.L. Egido, J. Lessing, V. Martin, and L.M. Robledo, “On the solution of the Hartree-Fock-Bogoliubov equations by the conjugate gradient method,” Nucl. Phys. A 594, 70 – 86 (1995).
  * Robledo and Bertsch (2011) L. M. Robledo and G. F. Bertsch, “Application of the gradient method to Hartree-Fock-Bogoliubov theory,” Phys. Rev. C 84, 014312 (2011).
  * Golub and van der Vorst (2000) G. H. Golub and H. A. van der Vorst, “Eigenvalue computation in the 20th century,” Journal of Computational and Applied Mathematics 123, 35 – 65 (2000).
  * Stoitsov  _et al._ (2013) M.V. Stoitsov, N. Schunck, M. Kortelainen, N. Michel, H. Nam, E. Olsen, J. Sarich, and S. Wild, “Axially deformed solution of the Skyrme-Hartree–Fock–Bogoliubov equations using the transformed harmonic oscillator basis (II) HFBTHO v2.00d: A new version of the program,” Computer Physics Communications 184, 1592 – 1604 (2013).
  * Schunck  _et al._ (2012) N. Schunck, J. Dobaczewski, J. McDonnell, W. Satuła, J.A. Sheikh, A. Staszczak, M. Stoitsov, and P. Toivanen, “Solution of the Skyrme–Hartree–Fock–bogolyubov equations in the cartesian deformed harmonic-oscillator basis.: (VII) HFODD (v2.49t): A new version of the program,” Computer Physics Communications 183, 166 – 192 (2012).
  * Pei  _et al._ (2014) J. C. Pei, G. I. Fann, R. J. Harrison, W. Nazarewicz, Yue Shi, and S. Thornton, “Adaptive multi-resolution 3d Hartree-Fock-Bogoliubov solver for nuclear structure,” Phys. Rev. C 90, 024317 (2014).
  * Michel  _et al._ (2002) N. Michel, W. Nazarewicz, M. Płoszajczak, and K. Bennaceur, “Gamow shell model description of neutron-rich nuclei,” Phys. Rev. Lett. 89, 042502 (2002).
  * Bulgac and Forbes (2013) Aurel Bulgac and Michael McNeil Forbes, “Use of the discrete variable representation basis in nuclear physics,” Phys. Rev. C 87, 051301 (2013).
  * Ryssens  _et al._ (2015b) W. Ryssens, P.-H. Heenen, and M. Bender, “Numerical accuracy of mean-field calculations in coordinate space,” Phys. Rev. C 92, 064318 (2015b).
  * Dobaczewski  _et al._ (1984) J. Dobaczewski, H. Flocard, and J. Treiner, “Hartree-Fock-Bogolyubov description of nuclei near the neutron-drip line,” Nuclear Physics A 422, 103 – 139 (1984).
  * Dobaczewski  _et al._ (1996) J. Dobaczewski, W. Nazarewicz, T. R. Werner, J. F. Berger, C. R. Chinn, and J. Dechargé, “Mean-field description of ground-state properties of drip-line nuclei: Pairing and continuum effects,” Phys. Rev. C 53, 2809–2840 (1996).
  * Zhang  _et al._ (2013) Y. N. Zhang, J. C. Pei, and F. R. Xu, “Hartree-Fock-Bogoliubov descriptions of deformed weakly bound nuclei in large coordinate spaces,” Phys. Rev. C 88, 054305 (2013).
  * Bennaceur and Dobaczewski (2005) K. Bennaceur and J. Dobaczewski, “Coordinate-space solution of the Skyrme–Hartree–Fock–Bogolyubov equations within spherical symmetry. The program HFBRAD (v1.00),” Computer Physics Communications 168, 96 – 122 (2005).
  * Oberacker  _et al._ (2003) V. E. Oberacker, A. S. Umar, E. Terán, and A. Blazkiewicz, “Hartree-Fock-Bogoliubov calculations in coordinate space: Neutron-rich sulfur, zirconium, cerium, and samarium isotopes,” Phys. Rev. C 68, 064302 (2003).
  * Terán  _et al._ (2003) E. Terán, V. E. Oberacker, and A. S. Umar, “Axially symmetric Hartree-Fock-Bogoliubov calculations for nuclei near the drip lines,” Phys. Rev. C 67, 064314 (2003).
  * Blazkiewicz  _et al._ (2005) A. Blazkiewicz, V. E. Oberacker, A. S. Umar, and M. Stoitsov, “Coordinate space Hartree-Fock-Bogoliubov calculations for the zirconium isotope chain up to the two-neutron drip line,” Phys. Rev. C 71, 054321 (2005).
  * Pei  _et al._ (2008) J. C. Pei, M. V. Stoitsov, G. I. Fann, W. Nazarewicz, N. Schunck, and F. R. Xu, “Deformed coordinate-space Hartree-Fock-Bogoliubov approach to weakly bound nuclei and large deformations,” Phys. Rev. C 78, 064306 (2008).
  * (31) http://www.netlib.org/scalapack/ .
  * (32) http://www.nersc.gov/users/computational-systems/edison/ .
  * (33) http://www.nvidia.com/object/cuda_home_new.html .
  * van der Vorst and Melissen (1990) H. A. van der Vorst and J. B. M. Melissen, “A petrov-galerkin type method for solving $Axk=b$, where $A$ is symmetric complex,” IEEE Transactions on Magnetics 26, 706–708 (1990).
  * Yamamoto  _et al._ (2008) Susumu Yamamoto, Tomohiro Sogabe, Takeo Hoshi, Shao-Liang Zhang, and Takeo Fujiwara, “Shifted conjugate-orthogonal–conjugate-gradient method and its application to double orbital extended hubbard model,” Journal of the Physical Society of Japan 77, 114713 (2008).
  * Takayama  _et al._ (2006) R. Takayama, T. Hoshi, T. Sogabe, S.-L. Zhang, and T. Fujiwara, “Linear algebraic calculation of the Green’s function for large-scale electronic structure theory,” Phys. Rev. B 73, 165108 (2006).
  * Mizusaki  _et al._ (2010) Takahiro Mizusaki, Kazunari Kaneko, Michio Honma, and Tetsuya Sakurai, “Filter diagonalization of shell-model calculations,” Phys. Rev. C 82, 024310 (2010).
  * Shimizu  _et al._ (2016) Noritaka Shimizu, Yutaka Utsuno, Yasunori Futamura, Tetsuya Sakurai, Takahiro Mizusaki, and Takaharu Otsuka, “Stochastic estimation of nuclear level density in the nuclear shell model: An application to parity-dependent level density in 58Ni,” Physics Letters B 753, 13 – 17 (2016).
  * Sakurai and Sugiura (2003) Tetsuya Sakurai and Hiroshi Sugiura, “A projection method for generalized eigenvalue problems using numerical integration,” Journal of Computational and Applied Mathematics 159, 119 – 128 (2003).
  * Sakurai  _et al._ (2008) Tetsuya Sakurai, Yoshihisa Kodaki, Hiroto Tadano, Daisuke Takahashi, Mitsuhisa Sato, and Umpei Nagashima, “A parallel method for large sparse generalized eigenvalue problems using a gridrpc system,” Future Generation Computer Systems 24, 613 – 619 (2008).
  * Ikegami  _et al._ (2010) Tsutomu Ikegami, Tetsuya Sakurai, and Umpei Nagashima, “A filter diagonalization for generalized eigenvalue problems based on the sakurai–sugiura projection method,” Journal of Computational and Applied Mathematics 233, 1927 – 1936 (2010).
  * J. Dobaczewski  _et al._ (2002) J. Dobaczewski, W. Nazarewicz, and M.V. Stoitsov, “Nuclear ground-state properties from mean-field calculations,” Eur. Phys. J. A 15, 21–26 (2002).
  * Chabanat  _et al._ (1998) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, “A Skyrme parametrization from subnuclear to neutron star densities Part II. Nuclei far from stabilities,” Nuclear Physics A 635, 231 – 256 (1998).
  * Stetcu  _et al._ (2011) I. Stetcu, A. Bulgac, P. Magierski, and K. J. Roche, “Isovector giant dipole resonance from the 3D time-dependent density functional theory for superfluid nuclei,” Phys. Rev. C 84, 051309 (2011).
  * Stetcu  _et al._ (2015) I. Stetcu, C. A. Bertulani, A. Bulgac, P. Magierski, and K. J. Roche, “Relativistic Coulomb Excitation within the Time Dependent Superfluid Local Density Approximation,” Phys. Rev. Lett. 114, 012701 (2015).
  * Bulgac and Yu (2002) Aurel Bulgac and Yongle Yu, “Renormalization of the Hartree-Fock-Bogoliubov equations in the case of a zero range pairing interaction,” Phys. Rev. Lett. 88, 042504 (2002).
  * Yu and Bulgac (2003) Yongle Yu and Aurel Bulgac, “Energy density functional approach to superfluid nuclei,” Phys. Rev. Lett. 90, 222501 (2003).
  * Borycki  _et al._ (2006) P. J. Borycki, J. Dobaczewski, W. Nazarewicz, and M. V. Stoitsov, “Pairing renormalization and regularization within the local density approximation,” Phys. Rev. C 73, 044319 (2006).
  * Bertsch  _et al._ (2009) G. F. Bertsch, C. A. Bertulani, W. Nazarewicz, N. Schunck, and M. V. Stoitsov, “Odd-even mass differences from self-consistent mean field theory,” Phys. Rev. C 79, 034306 (2009).
  * Bulgac  _et al._ (2015) Aurel Bulgac, Michael McNeil Forbes, and Shi Jin, “Nuclear energy density functionals: What do we really know?” arXiv:1506.09195 (2015).
  * Dutra  _et al._ (2012) M. Dutra, O. Lourenço, J. S. Sá Martins, A. Delfino, J. R. Stone, and P. D. Stevenson, “Skyrme interaction and nuclear matter constraints,” Phys. Rev. C 85, 035201 (2012).
  * Bender  _et al._ (2003b) Michael Bender, Paul-Henri Heenen, and Paul-Gerhard Reinhard, “Self-consistent mean-field models for nuclear structure,” Rev. Mod. Phys. 75, 121–180 (2003b).
  * Bohr and Mottelson (1998) A. Bohr and B.R. Mottelson,  _Nuclear Structure_, Vol. 1 (World Scientific, Singapore, 1998).
  * Trefethen and Bau (1997) L.N. Trefethen and D. Bau,  _Numerical Linear Algebra_ (Society for Industrial and Applied Mathematics, Philadelphia, 1997).
  * (55) Aurel Bulgac, “Hartree-Fock-Bogoliubov approximation for finite systems,” CIP, Bucharest, Report No. FT-194-1980; nucl-th/9907088 .
  * Gu  _et al._ (2014) X. M. Gu, T. Z. Huang, L. Li, H. B. Li, T. Sogabe, and M. Clemens, “Quasi-minimal residual variants of the COCG and COCR methods for complex symmetric linear systems in electromagnetic simulations,” IEEE Transactions on Microwave Theory and Techniques 62, 2859–2867 (2014).
  * Li  _et al._ (2011) Yue Hui Li, Zai Ping Nie, Min Meng, Xiang Qian Zhang, and Xiang Yang Sun, “An efficient MAINV preconditioned COCG method for FEM analysis of millimeter wave filters,” Journal of Infrared, Millimeter, and Terahertz Waves 32, 216–224 (2011).
  * Li  _et al._ (2008) L. Li, T.-Z. Huang, and Z.-G. Ren, “A preconditioned COCG method for solving complex symmetric linear systems arising from scattering problems,” Journal of Electromagnetic Waves and Applications 22, 2023–2034 (2008).
  * Wang  _et al._ (2009) Yingxi Wang, Kui Du, and Weiwei Sun, “Preconditioning iterative algorithm for the electromagnetic scattering from a large cavity,” Numerical Linear Algebra with Applications 16, 345–363 (2009).
  * Aichinger and Krotscheck (2005) M. Aichinger and E. Krotscheck, “A fast configuration space method for solving local Kohn–Sham equations,” Computational Materials Science 34, 188 – 212 (2005).
  * Frommer (2003) A. Frommer, “Bicgstab(ℓ) for families of shifted linear systems,” Computing 70, 87–109 (2003).
  * Frommer and Glässner (1998) Andreas Frommer and Uwe Glässner, “Restarted GMRES for shifted linear systems,” SIAM Journal on Scientific Computing 19, 15–26 (1998).
  * (63) https://www.olcf.ornl.gov/computing-resources/titan-cray-xk7/ .
  * Baran  _et al._ (2008) Andrzej Baran, Aurel Bulgac, Michael McNeil Forbes, Gaute Hagen, Witold Nazarewicz, Nicolas Schunck, and Mario V. Stoitsov, “Broyden’s method in nuclear structure calculations,” Phys. Rev. C 78, 014318 (2008).
  * Baye (2015) Daniel Baye, “The lagrange-mesh method,” Physics Reports 565, 1 – 107 (2015).
  * Younes and Gogny (2009) W. Younes and D. Gogny, “Microscopic calculation of ${}^{240}\mathrm{Pu}$ scission with a finite-range effective force,” Phys. Rev. C 80, 054313 (2009).
  * Schunck  _et al._ (2014) N. Schunck, D. Duke, H. Carr, and A. Knoll, “Description of induced nuclear fission with Skyrme energy functionals: Static potential energy surfaces and fission fragment properties,” Phys. Rev. C 90, 054305 (2014).
  * Regnier  _et al._ (2016) D. Regnier, N. Dubray, N. Schunck, and M. Verrière, “Fission fragment charge and mass distributions in ${}^{239}\mathrm{Pu}(n,f)$ in the adiabatic nuclear energy density functional theory,” Phys. Rev. C 93, 054611 (2016).
  * Staszczak  _et al._ (2010) A. Staszczak, M. Stoitsov, A. Baran, and W. Nazarewicz, “Augmented lagrangian method for constrained nuclear density functional theory,” Eur. Phys. J. A 46, 85–90 (2010).
  * Fayans (1998) S. A. Fayans, “Towards a universal nuclear density functional,” Journal of Experimental and Theoretical Physics Letters 68, 169–174 (1998).
  * Fayans  _et al._ (2000) S.A. Fayans, S.V. Tolokonnikov, E.L. Trykov, and D. Zawischa, “Nuclear isotope shifts within the local energy-density functional approach,” Nuclear Physics A 676, 49 – 119 (2000).
  * Shlomo and Bertsch (1975) S. Shlomo and G. Bertsch, “Nuclear response in the continuum,” Nuclear Physics A 243, 507 – 518 (1975).
  * Bertsch and Tsai (1975) G.F. Bertsch and S.F. Tsai, “A study of the nuclear response function,” Physics Reports 18, 125 – 158 (1975).
  * Nakatsukasa  _et al._ (2007) Takashi Nakatsukasa, Tsunenori Inakura, and Kazuhiro Yabana, “Finite amplitude method for the solution of the random-phase approximation,” Phys. Rev. C 76, 024318 (2007).
  * Avogadro and Nakatsukasa (2011) Paolo Avogadro and Takashi Nakatsukasa, “Finite amplitude method for the quasiparticle random-phase approximation,” Phys. Rev. C 84, 014314 (2011).
  * Nakatsukasa (2014) Takashi Nakatsukasa, “Finite amplitude method in linear response TDDFT calculations,” Journal of Physics: Conference Series 533, 012054 (2014).
  * Hinohara  _et al._ (2013) Nobuo Hinohara, Markus Kortelainen, and Witold Nazarewicz, “Low-energy collective modes of deformed superfluid nuclei within the finite-amplitude method,” Phys. Rev. C 87, 064309 (2013).
  * Oishi  _et al._ (2016) Tomohiro Oishi, Markus Kortelainen, and Nobuo Hinohara, “Finite amplitude method applied to the giant dipole resonance in heavy rare-earth nuclei,” Phys. Rev. C 93, 034329 (2016).
  * Mustonen  _et al._ (2014) M. T. Mustonen, T. Shafer, Z. Zenginerler, and J. Engel, “Finite-amplitude method for charge-changing transitions in axially deformed nuclei,” Phys. Rev. C 90, 024308 (2014).