# Solving the three-dimensional Skyrme Hartree-Fock-Bogoliubov problem using
the mixed-basis method

Yue Shi (石跃) corresponding author<EMAIL_ADDRESS>Department of Physics,
Harbin Institute of Technology, Harbin 150001, People’s Republic of China
Nobuo Hinohara  Center for Computational Sciences, University of Tsukuba,
Tsukuba 305-8577, Japan  Faculty of Pure and Applied Sciences, University of
Tsukuba, Tsukuba 305-8571, Japan

###### Abstract

Background

The symmetry-unrestricted Hartree-Fock-Bogoliubov (HFB) simulation is
important for describing various quantum many-body systems. However, the HFB
problem in Cartesian coordinate space is numerically challenging.

Purpose

For describing ground states without imposing axial symmetry and looking ahead
to future extension for dynamics with full time dependence, we present a
numerically efficient implementation of the three-dimensional (3D) HFB code.

Methods

We develop a 3D Skyrme HFB code based on the mixed-basis representation
(HFBmix) which consists of two harmonic-oscillator (HO) bases in the $x$\- and
$y$-directions, and finite-difference (FD) basis in the $z$-direction in
solving the nuclear 3D HFB problem.

Results

The results show very well agreement among all the three codes (HFBmix, HO3D,
and hfodd). Especially for the HF calculations, the differences in total
energies are on the order of a few keV for the lightest O and Mg nuclei. The
HFBmix is applied to spherical, prolate, and triaxial systems, and gives the
same quadrupole moments for the deformed nuclei as those of the HO-based
calculations. Feasibility of the HFBmix is demonstrated in the fission isomer
and barrier calculations of 240Pu.

Conclusions

The HFBmix is useful for solving the nuclear 3D HFB problem for its numerical
efficiency. Future work will include the analysis of deformed drip-line
systems and systematic potential-energy surface calculation for fission-path
analysis as well as the time-dependent extension of the HFBmix code for
dynamics calculations.

##  I Introduction

For various quantum many-body simulations in quantum chemistry, condensed-
matter physics, atomic physics, and nuclear physics, one deals with problems
where constituent particles extend far spatially [1]. For example, inside a
crystal lattice, the wave functions of the electrons spread extensively within
the Coulomb potential exerted by the atomic nuclei which are distributed
periodically. Similarly, a composite molecule often constitutes many atoms
situated at a significant distance. While the electrons of such a molecule may
be localized around a certain atom in the ground state, they may be extended
over the whole molecule in the excited states. For some drip-line nuclei where
the Fermi energy is close to zero, the dilute nuclear material densities may
extend spatially when constituent nucleons occupy orbits with low orbital
angular momenta (halo nuclei). In these situations, a microscopic description
for the many-body wave functions of the quantum systems requires the solution
in a three-dimensional (3D) coordinate/momentum space.

In theoretical nuclear physics, a multitude of independent-particle (IP)
models play important roles. The conventional IP models use the well-tuned
modified harmonic oscillator (HO) or the Woods-Saxon potential plus the spin-
orbit (SO) potential to reproduce the experimental magic numbers [2]. Such
mean-field models together with the Strutinski shell correction method [3]
achieve quantitative success in describing a variety of observables in low-
energy nuclear physics. These models are known as microscopic-macroscopic
models.

Starting from an effective nucleon-nucleon interaction or Lagrangian together
with a variational Hartree-Fock (HF) procedure, one determines the one-body
potential in a self-consistent manner [4]. These IP models are known as the
self-consistent mean-field theories which are based on the zero-range Skyrme
[5, 6] or the finite-range Gogny [7, 8, 9] forces, and the relativistic mean-
field theory [10, 11, 12, 13]. With the advances in synthesizing the exotic
nuclei with extreme $N/Z$ ratio, modern self-consistent mean-field theories
are more capable of accounting for new features which were not present in
nuclei near the valley of stability, such as halo nuclei mentioned above. This
requires the models to be flexible enough to work in a 3D coordinate space,
yet efficient enough for the beyond-mean-field extensions and to cover a large
parameter space [14, 4].

However, the requirements are, in general, numerically formidable with the
Hartree-Fock-Bogoliubov (HFB) treatment of the pairing correlation. Indeed, in
a typical nuclear 3D coordinate-space HFB simulation, one needs to repeatedly
diagonalize the Hamiltonian matrix of the dimension of a few hundreds of
thousands [15]. Recent application of the shifted Krylov method [16, 17] and
expansion in the three-dimensional wave-number space basis [18] accelerate the
3D symmetry-unrestricted HFB calculations. In comparison, one- [6, 19, 20] or
two-dimensional [21, 22, 23] HFB calculations in coordinate space are less
expensive due to their restricted spatial degrees of freedom.

In realistic applications, one frequently encounters situations where one- or
two-dimensional degrees of freedom play a major role in the static or
dynamical phenomena which do not possess any spatial symmetry. For example, a
fissioning nucleus assumes rather elongated deformation near or beyond the
scission point. The inner core of a neutron star may become so dense that the
“nuclei” are packed closed to each other with various exotic deformations
favoring only one or two spacial degrees of freedom. These are known as the
nuclear pasta [24, 25] which have been discussed with various Hartree-Fock
calculations in coordinate space [26, 27, 28]. In such situations it is more
sensible to adopt a coordinate basis in the elongated direction(s) and use HO
bases for the remaining directions where the matter is more restricted.

Indeed, a similar idea of using a mixed basis has been realized in the nuclear
finite-range Gogny HFB method and its time-dependent extension [29, 30]. The
current work intends to describe the implementation of the nuclear Skyrme HFB
method using the HO plus finite-difference (HO+FD) mixed basis. Throughout
this work, we refer to the basis which consists of the two HO bases in the
$x$\- and $y$-directions, and the FD treatment in the $z$-direction as the
“mixed basis.”

In Sec. II, we describe the nuclear HFB problem and the numerical details for
the HO+FD mixed-basis HFB calculation. Section III presents a systematic
comparison for the results using the mixed-basis code and those using a 3D HO
based code (HO3D) developed in this work and the hfodd codes (version 2.49t)
[31]. The summary and perspective are included in Sec. IV.

##  II The model

In this section, we provide a short description of the nuclear HFB theory. We
then provide the numerical details for the implementations of the HF(B)mix
code and the parameters chosen for the calculations.

###  II.1 The Skyrme HFB formalism

In the nuclear Skyrme HF mean-field theory, the total energy $E_{\rm Total}$
of a nucleus can be decomposed into the kinetic, Skyrme, pairing, and Coulomb
terms:

$\displaystyle E_{\rm Total}$ $\displaystyle=$ $\displaystyle E_{\rm
kin+c.m.}+E_{\rm Skyrme}+E_{\rm pair}+E_{\rm Coul}$ (1) $\displaystyle=$
$\displaystyle\int
d^{3}r\big{[}\mathcal{K}(\mbox{\boldmath$r$})+\mathcal{E}_{\rm
Skyrme}(\mbox{\boldmath$r$})$ $\displaystyle\quad+\mathcal{E}_{\rm
pair}(\mbox{\boldmath$r$})+\mathcal{E}_{\rm
Coul}(\mbox{\boldmath$r$})\big{]}.$

The energy densities, $\mathcal{K}(\mbox{\boldmath$r$})$, $\mathcal{E}_{\rm
Skyrme}(\mbox{\boldmath$r$})$, $\mathcal{E}_{\rm pair}(\mbox{\boldmath$r$})$,
and $\mathcal{E}_{\rm Coul}(\mbox{\boldmath$r$})$, are functionals of various
local densities (scalar density $\rho$, kinetic density $\tau$, spin-current
density $\mathrm{J}$, and pair density $\tilde{\rho}$). The kinetic-energy
density of the nucleus with the one-body center-of-mass correction is given by

$\displaystyle\mathcal{K}(\bm{r})=\frac{\hbar^{2}}{2m}\tau\quantity(1-\frac{1}{A}),$
(2)

where $m$ is the nucleon mass and $A$ is the mass number. The time-even part
of the Skyrme energy density functional reads

$\displaystyle\mathcal{E}_{\rm Skyrme}(\bm{r})$ $\displaystyle=$
$\displaystyle\frac{b_{0}}{2}\rho^{2}-\frac{b^{\prime}_{0}}{2}\sum_{q}\rho^{2}_{q}+b_{1}\rho\tau-b^{\prime}_{1}\sum_{q}\rho_{q}\tau_{q}$
(3)
$\displaystyle\quad-\frac{b_{2}}{2}\rho\nabla^{2}\rho+\frac{b^{\prime}_{2}}{2}\sum_{q}\rho_{q}\nabla^{2}\rho_{q}+\frac{b_{3}}{3}\rho^{\alpha+2}-\frac{b^{\prime}_{3}}{3}\rho^{\alpha}\sum_{q}\rho^{2}_{q}$
$\displaystyle\quad-
b_{4}\rho\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath$\mathrm{J}$}-b^{\prime}_{4}\sum_{q}\rho_{q}\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath$\mathrm{J}$}_{q}.$

The index $q$ denotes neutrons ($n$) or protons ($p$). The local densities
without the index $q$ indicate the sum of neutron and proton densities. The
current work assumes the time-reversal symmetry and discusses only even-even
systems. Therefore, the time-odd part of the energy density functional does
not contribute in Eq. (3). The Coulomb energy is decomposed into the direct
and exchange terms

$\mathcal{E}_{\rm Coul}(\mbox{\boldmath$r$})=\mathcal{E}_{\rm Coul}^{\rm
Dir.}(\mbox{\boldmath$r$})+\mathcal{E}_{\rm Coul}^{\rm
Exc.}(\mbox{\boldmath$r$}),$ (4)

where the direct term is given by

$\mathcal{E}_{\rm Coul}^{\rm Dir.}(\mbox{\boldmath$r$})=\frac{e^{2}}{2}\int
d^{3}r^{\prime}\frac{\rho_{p}(\mbox{\boldmath$r$})\rho_{p}(\mbox{\boldmath$r$}^{\prime})}{|\mbox{\boldmath$r$}-\mbox{\boldmath$r$}^{\prime}|},$
(5)

and the exchange term is treated in the Slater approximation as

$\mathcal{E}_{\rm Coul}^{\rm
Exc.}(\mbox{\boldmath$r$})=-\frac{3e^{2}}{4}\quantity(\frac{3}{\pi})^{1/3}\rho_{p}^{4/3}(\mbox{\boldmath$r$}).$
(6)

The volume-type pairing energy density is

$\displaystyle\mathcal{E}_{\rm
pair}(\mbox{\boldmath$r$})=\sum_{q}\frac{V^{q}_{0}}{4}\tilde{\rho}^{2}_{q}(\mbox{\boldmath$r$}).$
(7)

The various local densities are constructed using the quasi-particle wave
functions that are solutions of the HFB equation [2, 6]

$\matrixquantity(h_{\uparrow\uparrow}-\lambda&h_{\uparrow\downarrow}&\tilde{h}_{\uparrow\uparrow}&\tilde{h}_{\uparrow\downarrow}\\\
h_{\downarrow\uparrow}&h_{\downarrow\downarrow}-\lambda&\tilde{h}_{\downarrow\uparrow}&\tilde{h}_{\downarrow\downarrow}\\\
\tilde{h}_{\uparrow\uparrow}&\tilde{h}_{\uparrow\downarrow}&-h_{\uparrow\uparrow}+\lambda&-h_{\uparrow\downarrow}\\\
\tilde{h}_{\downarrow\uparrow}&\tilde{h}_{\downarrow\downarrow}&-h_{\downarrow\uparrow}&-h_{\downarrow\downarrow}+\lambda)\matrixquantity(u_{\uparrow,k}\\\
u_{\downarrow,k}\\\ v_{\uparrow,k}\\\
v_{\downarrow,k})=E_{k}\matrixquantity(u_{\uparrow,k}\\\ u_{\downarrow,k}\\\
v_{\uparrow,k}\\\ v_{\downarrow,k}),$ (8)

where we have omitted the superscripts, $q$’s, from all quantities for
simplicity. The arrows, $\uparrow$ and $\downarrow$, denote the nucleon spin
$\sigma=1/2$ and $\sigma=-1/2$, respectively. For the four-component quasi-
particle wave functions, we have abbreviated that $u_{\uparrow,k}\equiv
u_{k}(\mbox{\boldmath$r$}\sigma=\frac{1}{2})$, for instance.

In Eq. (8) the particle-hole part of the mean-field Hamiltonian is

$\displaystyle h^{q}_{\sigma\sigma^{\prime}}(\mbox{\boldmath$r$})$
$\displaystyle=$
$\displaystyle\quantity[-\mbox{\boldmath$\nabla$}\cdot\frac{\hbar^{2}}{2m^{*}_{q}(\mbox{\boldmath$r$})}\mbox{\boldmath$\nabla$}+U_{q}(\mbox{\boldmath$r$})+U_{\rm
Coul}(\mbox{\boldmath$r$})\,\delta_{q,p}]\delta_{\sigma\sigma^{\prime}}$ (9)
$\displaystyle-\quantity[i\mathbf{B}_{q}(\mbox{\boldmath$r$})\cdot\quantity(\mbox{\boldmath$\nabla$}\crossproduct\mbox{\boldmath$\sigma$})]_{\sigma\sigma^{\prime}},$

where $\bm{\sigma}$ is a Pauli matrix. The effective mass is defined through

$\frac{\hbar^{2}}{2m^{*}_{q}(\mbox{\boldmath$r$})}=\frac{\hbar^{2}}{2m}+b_{1}\rho-
b_{1}^{\prime}\rho_{q},$ (10)

and

$\mathbf{B}_{q}(\mbox{\boldmath$r$})=b_{4}\mbox{\boldmath$\nabla$}\rho(\mbox{\boldmath$r$})+b^{\prime}_{4}\mbox{\boldmath$\nabla$}\rho_{q}(\mbox{\boldmath$r$}).$
(11)

The potential due to Skyrme force in Eq. (9) reads

$\displaystyle U_{q}(\mbox{\boldmath$r$})$
$\displaystyle=b_{0}\rho-b^{\prime}_{0}\rho_{q}+b_{1}\tau-b^{\prime}_{1}\tau_{q}-b_{2}\nabla^{2}\rho+b^{\prime}_{2}\nabla^{2}\rho_{q}$
$\displaystyle\quad+\frac{b_{3}}{3}\quantity(\alpha+2)\rho^{\alpha+1}-\frac{b^{\prime}_{3}}{3}\sum_{q}\quantity(\alpha\rho^{\alpha-1}\rho^{2}_{q}+2\rho^{\alpha}\rho_{q})$
$\displaystyle\quad-
b_{4}\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath$\mathrm{J}$}-b^{\prime}_{4}\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath$\mathrm{J}$}_{q}.$
(12)

For protons ($q=p$), the Coulomb potential in Eq. (9) is composed of the
direct and exchange parts

$U_{\rm Coul}(\mbox{\boldmath$r$})=U_{\rm Coul}^{\rm
Dir.}(\mbox{\boldmath$r$})+U_{\rm Coul}^{\rm Exc.}(\mbox{\boldmath$r$}),$ (13)

where

$U_{\rm Coul}^{\rm
Exc.}(\mbox{\boldmath$r$})=-e^{2}\quantity(\frac{3}{\pi})^{1/3}\quantity[\rho_{p}(\mbox{\boldmath$r$})]^{1/3}.$
(14)

The direct part of the Coulomb potential is obtained by solving the 3D Poisson
equation

$\nabla^{2}U_{\rm Coul}^{\rm Dir.}(\mbox{\boldmath$r$})=-4\pi
e^{2}\rho_{p}(\mbox{\boldmath$r$}).$ (15)

For details, see Ref. [32]. Upon obtaining $U_{\rm Coul}^{\rm Dir.}$, the
direct Coulomb energy density (5) can be calculated through

$\mathcal{E}_{\rm Coul}^{\rm
Dir.}(\mbox{\boldmath$r$})=\frac{1}{2}\rho_{p}(\mbox{\boldmath$r$})\,U_{\rm
Coul}^{\rm Dir.}(\mbox{\boldmath$r$}).$ (16)

The pairing mean-field Hamiltonian in Eq. (8) reads

$\tilde{h}^{q}_{\sigma\sigma^{\prime}}(\mbox{\boldmath$r$})=\frac{1}{2}V^{q}_{0}\tilde{\rho}_{q}(\mbox{\boldmath$r$})\delta_{\sigma\sigma^{\prime}}.$
(17)

The local densities appearing in Eqs. (2), (3), (5), (6), (7), (9), and (17)
are expressed in terms of the quasi-particle wave functions of the HFB
equation (8)

$\displaystyle\rho_{q}(\mbox{\boldmath$r$})=\sum_{k}\sum_{\sigma=\pm\frac{1}{2}}\quantity|v_{k,q}(\mbox{\boldmath$r$}\mathbf{\sigma})|^{2},$
$\displaystyle\tau_{q}(\mbox{\boldmath$r$})=\sum_{k}\sum_{\sigma=\pm\frac{1}{2}}\quantity|\bm{\nabla}v_{k,q}(\mbox{\boldmath$r$}\mathbf{\sigma})|^{2},$
$\displaystyle\tilde{\rho}_{q}(\mbox{\boldmath$r$})=-\sum_{k}\sum_{\sigma=\pm\frac{1}{2}}v_{k,q}(\mbox{\boldmath$r$}\sigma)u_{k,q}^{*}(\mbox{\boldmath$r$}\sigma),$
$\displaystyle\mathbf{J}_{q}(\mbox{\boldmath$r$})=-i\sum_{k}v^{\dagger}_{k,q}(\mbox{\boldmath$r$})\quantity(\mbox{\boldmath$\nabla$}\crossproduct\mbox{\boldmath$\sigma$})v_{k,q}(\mbox{\boldmath$r$}).$
(18)

When evaluating the spin-current density $\mathbf{J}_{q}$ on the right hand
side of Eq. (18), we assume the form of
$v_{k,q}(\mbox{\boldmath$r$})\equiv\matrixquantity(v_{k,q}(\mbox{\boldmath$r$},\frac{1}{2})\\\
v_{k,q}(\mbox{\boldmath$r$},-\frac{1}{2}))$. The summations over $k$ are
limited to the states within the energy window of $0<E_{k}<60$ MeV. This is to
make sure the lowest hole state is included.

In practical calculations, one starts with a set of quasi-particle wave
functions [$u_{k}(\mbox{\boldmath$r$}\sigma)$’s and
$v_{k}(\mbox{\boldmath$r$}\sigma)$’s] calculated assuming a Woods-Saxon
potential in Eq. (II.1). The pairing field (17), which is diagonal in the
$\sigma$ spcace, is initialized with
$\tilde{h}^{q}_{\uparrow\uparrow}(\mbox{\boldmath$r$})=\tilde{h}^{q}_{\downarrow\downarrow}(\mbox{\boldmath$r$})=-3.0$
MeV for $|\mbox{\boldmath$r$}|\leq R$, and
$\tilde{h}^{q}_{\uparrow\uparrow}(\mbox{\boldmath$r$})=\tilde{h}^{q}_{\downarrow\downarrow}(\mbox{\boldmath$r$})=0$
for $|\mbox{\boldmath$r$}|>R$, where $R=1.2A^{1/3}$ fm. The obtained quasi-
particle wave functions are inserted into Eq. (II.1) to calculated the
densities. One then assembles various potentials and obtains the HFB matrix
[Eq. (8)] in some suitable representation. Diagonalizing it, one obtains the
new $u_{k}(\mbox{\boldmath$r$}\sigma)$’s and
$v_{k}(\mbox{\boldmath$r$}\sigma)$’s. This process is continued until
convergence is achieved. To ensure stable convergence, we use linearly mixed
densities with 75% of the densities from the last iteration.

The particle-hole [Eqs. (II.1) and (13)] and pairing [Eq. (17)] potentials are
diagonal in the coordinate-space representation. Their matrix forms in the HO
and FD basis can be obtained conveniently. In Appendix A, We show how the
matrix elements of various densities are calculated in the mixed basis. The
kinetic and the SO parts in Eq. (9) contain $\nabla$ and $\nabla^{2}$
operators, which are not diagonal in the coordinate-space representation. In
Appendices B and C, we will give expressions for the kinetic terms and SO
terms in the mixed basis, respectively.

###  II.2 Constraints

For the constrained calculation, we follow the procedure in Ref. [33]. The
essential part of the constrained calculation is to add the following
contribution to the single-particle potential (9)

$C\quantity(\langle\hat{O}\rangle-\mu)\hat{O},$ (19)

where $C$ is the stiffness coefficient; $\hat{O}$ is the multipole operator in
question; $\mu$ is changing in each iteration with

$\mu^{(i+1)}=\mu^{(i)}-0.02\quantity(\langle\hat{O}\rangle^{(i+1)}-\mu_{0}),$
(20)

where $\mu_{0}$ is the constraint value at the first iteration. This procedure
is similar to the linear constraint method, and turns out to give fairly
accurate target constrained values at convergence.

In addition to the principal-axes constraint operators described in the next
subsection, for the multipole operators appears in this work, we use the
quadrupole and octupole operators:

$\displaystyle\hat{Q}_{20}$ $\displaystyle=$ $\displaystyle
2\hat{z}^{2}-\hat{x}^{2}-\hat{y}^{2},$ $\displaystyle\hat{Q}_{22}$
$\displaystyle=$ $\displaystyle\sqrt{3}\quantity(\hat{x}^{2}-\hat{y}^{2}),$
$\displaystyle\hat{Q}_{30}$ $\displaystyle=$
$\displaystyle\hat{r}^{3}Y^{*}_{30}(\hat{\Omega}).$ (21)

The $\beta_{2}$ and $\gamma$ deformations which measure the degree of axial
and triaxial deformations are defined as

$\displaystyle\beta_{2}$ $\displaystyle=$
$\displaystyle\frac{4\pi}{3AR^{2}}\sqrt{Q_{20}^{2}+Q_{22}^{2}},$
$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle\atan{Q_{20}/Q_{22}},$

where $R=1.2A^{1/3}$ fm, $Q_{20}=\langle\hat{Q}_{20}\rangle$, and
$Q_{22}=\langle\hat{Q}_{22}\rangle$.

###  II.3 Numerical details

As noted at the end of Sec. II.1, after solving the HFB problem in the mixed
basis, one obtains a set of eigenvalues and respective eigenfunctions in the
HO+FD mixed-basis representation,
$(E_{k},u^{n_{x},n_{y},i_{z}}_{k},v^{n_{x},n_{y},i_{z}}_{k})$, where we have
ignored the spin and isospin for simplicity. The superscripts of $u$ and $v$
indicate the fact that they are vectors in the $|n_{x}n_{y}i_{z}\rangle$
space. The eigenvectors in the 3D coordinate space can then be obtained using
the following transform

$\displaystyle
u_{k}(i_{x},i_{y},i_{z})=\sum_{n_{x}}\sum_{n_{y}}\psi_{n_{x}}(x_{i_{x}})\psi_{n_{y}}(y_{i_{y}})\,u^{n_{x},n_{y},i_{z}}_{k},$
(22) $\displaystyle
v_{k}(i_{x},i_{y},i_{z})=\sum_{n_{x}}\sum_{n_{y}}\psi_{n_{x}}(x_{i_{x}})\psi_{n_{y}}(y_{i_{y}})\,v^{n_{x},n_{y},i_{z}}_{k},$
(23)

where

$\mu_{i_{\mu}}=(i_{\mu}-0.5)\times d\mu,~{}~{}\mu=x,y,z.$ (24)

We have used the discretized coordinates $i_{x},i_{y},i_{z}$ to mirror the
situation in the $x$ and $y$ directions. In Eq. (22), $\psi_{n}(\mu)$ are the
HO basis functions [34]. We will call this representation “grid basis” later.

In the implemented HFB code with the mixed basis we use two numbers to
characterize the dimension of the basis $|n_{x}n_{y}i_{z}\rangle$, which are
$N_{\rm max}$ and $N_{z}$. 111For HO3D and hfodd calculations, we also use
$N_{\rm max}$ to specify the largest $N$ for the spherical HO basis
$(n_{x}+n_{y}+n_{z}\leq N_{\rm max})$. In some hfodd calculations, we use
$N^{x,y,z}_{\rm max}$ to specify the largest $N$ in the $x$-, $y$-, and
$z$-direction of the non-uniform HO basis. The HO basis number is determined
by requiring $n_{x}+n_{y}\leq N_{\rm max}$. The value of $N_{z}$ specifies the
number of grid points in the $z$ direction. Using this basis, the dimension is
$\frac{(N_{\rm max}+1)\times(N_{\rm max}+2)}{2}\times N_{z}\times 2$ for the
HF calculations. For HFB calculation, the dimension should be doubled. The
spacing between grid points is given by $dz$. It is clear that $N_{\rm max}$,
$N_{z}$, and $dz$ uniquely specify our mixed basis. We always select the value
of $N_{z}$ to make sure the simulating box is large enough so that the neutron
density at the edge is smaller than 10-5 fm-3. We use nine-point formulae for
both the first- and second-derivative operators with respect to $z$ which will
be detailed in Appendix A. For the HO basis, we use an oscillator constant of

$b_{\mu}=\sqrt{\frac{41\times 0.6}{A^{\frac{1}{3}}\times 20.7355}}~{}{\rm
fm}^{-1},~{}~{}~{}\mu=x,y,z.$ (25)

To minimize the difference when comparing the results from HFBmix with those
from hfodd, we will readjust the pairing strengths of hfodd calculations to
match the pairing energy [Eq. (7)] of HFBmix. To check consistency of such a
procedure between the two calculations, another useful quantity is the average
pairing gap, which is defined as

$\Delta_{q}=-\frac{1}{2N_{q}}\int
d^{3}r\sum_{\sigma=\pm\frac{1}{2}}\tilde{h}^{q}_{\sigma\sigma}(\mbox{\boldmath$r$})\rho_{q}(\mbox{\boldmath$r$}),$
(26)

where $N_{q}$ denotes the proton or neutron number.

For the Coulomb potential, we solve the Poisson equation with two points
outside the simulating box as the boundary condition. In the current work, we
adopt the Dirichlet boundary condition. Using very similar choices of boundary
conditions and discretization method (FD method) in the $z$ direction, we note
a good description of spatial vibrations of nucleus [35].

Ideally, one can impose $y$-simplex symmetry for the wave functions, rendering
a block diagonal structure of the HFB matrix in Eq. (8). However, at large
deformation, like the case that will be shown in Sec. III.2.4 for the
fissioning 240Pu, the broken symmetries are encountered often. To allow for
future time-dependent extensions capable of complex simulations at large
deformation, the currently implemented mixed-basis HFB problem does not impose
any spatial symmetry.

If the initial quasi-particle wave functions possess certain symmetry, for
instance reflection symmetry, the iteration procedure will preserve the
symmetry (known as the self-consistent symmetry [2]). Due to the self-
consistent symmetry, most of our calculations which are started from a set of
wave functions conserving parity will continue to conserve parity symmetry. We
always make sure that the nucleus’ principal axes coincide with the Cartesian
system by checking the expectation values of the $\hat{x}$, $\hat{y}$,
$\hat{z}$, $\hat{x}\hat{y}$, $\hat{x}\hat{z}$, and $\hat{y}\hat{z}$ operators
are close to zero. For instance, the center of mass of the triaxially deformed
ground state of 110Mo deviates from the origin by only $\langle
x,y,z\rangle<10^{-13}$ fm at convergence. Thus, except for the calculation of
the fissioning path of 240Pu where the parity symmetry is explicitly broken by
the $\hat{Q}_{30}$-constraining term, we do not add constraints to fix the
orientation or the center of mass of the nucleus. When the axial octupole
moment ($Q_{30}$) of 240Pu deviates from zero, the dipole moment (with
operator $\hat{z}$) is constrained to be zero in order to prevent the center
of mass of the nucleus from moving around.

##  III Results and Discussions

The HF results provide perhaps the most suitable testing cases for
benchmarking codes using different representations. This is because, in the HF
case, the occupied single-particle orbits, which are bound, can be made
identical for the codes using different representations, whereas in the HFB
case, this can hardly be achieved. Indeed, to include all the deep-hole states
in HFB, one has to include the continuum state in the configuration space. The
description of the quasi-particle states in the HFB depends on how the
continuum state is discretized in the employed basis representation. This
point will be illustrated in Sec. III.2. For instance, the HFB calculation for
120Sn with the mixed basis includes 676 proton and 772 neutron quasi-particle
states for $0<E_{k}<60$ MeV. In contrast, there are 574 proton and 616 neutron
quasi-particle states using the 3D HO basis.

###  III.1 HF Results

In this subsection, we perform calculations for two spherical nuclei (16O and
208Pb), one prolately deformed nucleus (24Mg), and one triaxially deformed
nucleus (64Ge) using the HFmix code. We compare the results with those of the
hfodd code to show the accuracy of the developed mixed-basis code on the HF
level. Before evaluating the FD part of the HFBmix code, careful comparisons
are made for 16O between HO3D and hfodd codes to check the accuracy of the HO
part of the HFBmix code.

####  III.1.1 HO3D v.s. hfodd

Before constructing the HO+FD mixed-basis code, we first implement a 3D HO-
basis code, HO3D. This effort is not to have another 3D HO-basis code similar
to hfodd, but the HO3D code should be considered as an intermediate byproduct
leading to the final mixed-basis code. In this subsection, we first compare
the results of HO3D with those of hfodd code. The purpose is to check the
implementation of the HO part of the mixed-basis code. In Table 1, we list the
results of 16O calculated using HO3D and hfodd codes in the same HO model
space with $N_{\rm max}=14$.

In HO3D, we used the simple trapezoidal rule for the integration over the
coordinates in the HO basis. We did not employ the Gauss-Hermite quadrature
which is employed in hfodd code, because the basis in the $z$ direction of the
mixed-basis code is expressed in the equal-distance grid points, and it is
more straightforward to implement the same numerical technique for the
integrals and derivatives in $x$, $y$, and $z$ directions.

Even with the above-mentioned differeces, we see that the agreement between
the two codes on the total energy is less than $1$ keV for $dz=0.6$ fm. In an
earlier study [32], the HO wave functions are represented in the FD method and
the integrals are calculated using the trapezoidal rule. Comparing Table 1
with Table II of Ref. [32], we observe a similar convergence pattern.

Table 1: The total energy of 16O and its decomposition into $E_{\rm kin+c.m.}$, $E_{\rm Skyrme}$, and $E_{\rm Coul}$ (in MeV) calculated in $N_{\rm max}=14$ model space using hfodd and HO3D with different spatial discretizations ($dz$) for integration. The $E_{\rm Skyrme}$ is further decomposed into $E_{\rho^{2}}$, $E_{\rho\tau}$, $E_{\rho^{2+\alpha}}$, $E_{\rho\nabla^{2}\rho}$, and $E_{\rm SO}$ which correspond to the space integrals of terms with coefficients ($b_{0}$, $b_{0}^{\prime}$), ($b_{1}$, $b_{1}^{\prime}$), ($b_{3}$, $b_{3}^{\prime}$), ($b_{2}$, $b_{2}^{\prime}$), and ($b_{4}$, $b_{4}^{\prime}$) in Eq. (3), respectively. Skyrme SLy4 force is employed. | HO3D ($N_{\rm max}=14$) | hfodd
---|---|---
|  $dz=0.9$ fm |  $dz=0.833$ fm |  $dz=0.75$ fm |  $dz=0.6$ fm | ($N_{\rm max}=14$)
$E_{\rm Total}$ |  $-$128.40483 |  $-$128.46119 |  $-$128.46798 |  $-$128.46573 |  $-$128.46614
$E_{\rm kin+c.m.}$ | 222.23732 | 222.16985 | 222.18856 | 222.17432 | 222.18044
$E_{\rho^{2}}$ |  $-$1304.7338 |  $-$1303.3176 |  $-$1303.4058 |  $-$1303.3125 |  $-$1303.3417
$E_{\rho\tau}$ | 50.50790 | 50.45808 | 50.46788 | 50.46340 | 50.46476
$E_{\rho^{2+\alpha}}$ | 829.76110 | 828.75045 | 828.81217 | 828.74865 | 828.76895
$E_{\rho\nabla^{2}\rho}$ | 61.18723 | 60.84091 | 60.83204 | 60.82523 | 60.82684
$E_{\rm SO}$ |  $-$0.95678 |  $-$0.95063 |  $-$0.94971 |  $-$0.94994 |  $-$0.94962
$E_{\rm Coul}^{\rm Dir.}$ | 16.40879 | 16.40333 | 16.40257 | 16.40072 | 16.39990
$E_{\rm Coul}^{\rm Exc.}$ |  $-$2.81656 |  $-$2.81562 |  $-$2.81571 |  $-$2.81562 |  $-$2.81566

####  III.1.2 HFmix v.s. HFODD: Spherical nuclei

In Table 2 we present a set of HF results of 16O using HFmix and hfodd codes.
The model space for hfodd is $N_{\rm max}=11$, while it is $N_{\rm max}=11$
for the HO part of HFmix. Five $dz$ values (1.0, 0.9, 0.833, 0.75, and 0.6 fm)
are used to see the convergence property of the code with decreasing grid
spacing. For each grid spacings, the $N_{z}$ values are listed in the table.
This corresponds to $z_{\rm max}=(N_{z}-1)\times dz/2\approx 8$ fm. Because at
the edge of the boxes the neutron or proton densities are well below 10-6
fm-3, the size of $N_{z}$ is not the main factor affecting the precision of
these calculations.

For the HFmix results, we see a good convergence of the total energy with
decreasing $dz$. For $dz=0.9$ fm, the total energy overbinds by only $<50$ keV
compared to the result of $dz=0.6$ fm. These mixed-basis results show very
small quadrupole moments ($Q_{20}<0.1$ fm2) due to the breaking of the
spherical symmetry in the basis. The smallness of the quadrupole deformation
calculated for 16O can be considered as a good check for the correctness of
the code using a mixed basis. Comparing the HFmix results for 16O with hfodd
ones, we notice a good agreement between them. Indeed, the total energy
calculated by the hfodd code differs from those of HFmix ($dz=0.9$, 0.833,
0.75, and 0.6 fm) by only less than $50$ keV.

Table 2: The total energy of 16O and its decomposition into $E_{\rm kin+c.m.}$, various Skyrme SLy4 energies (defined in Table 1), and $E_{\rm Coul}$ as well as the quadrupole moments calculated using HFmix ($N_{\rm max}=11$) with different spatial discretizations ($dz$) and hfodd ($N_{\rm max}=11$). All values are in MeV, except for $Q_{20}$ and $\beta_{2}$ which are in fm2 and dimensionless. | HFmix($N_{\rm max}=11$) | hfodd
---|---|---
|  $dz=1.0$ fm |  $dz=0.9$ fm |  $dz=0.833$ fm |  $dz=0.75$ fm |  $dz=0.6$ fm | $N_{\rm max}=11$
| $N_{z}=18$ | $N_{z}=18$ | $N_{z}=22$ | $N_{z}=22$ | $N_{z}=30$ |
$E_{\rm Total}$ |  $-$128.570 |  $-$128.504 |  $-$128.490 |  $-$128.482 |  $-$128.476 |  $-$128.450
$E_{\rm kin+c.m.}$ | 222.523 | 222.315 | 222.276 | 222.236 | 222.199 | 222.348
$E_{\rho^{2}}$ |  $-$1305.807 |  $-$1304.498 |  $-$1303.864 |  $-$1303.577 |  $-$1303.272 |  $-$1303.658
$E_{\rho\tau}$ | 50.614 | 50.545 | 50.498 | 50.478 | 50.456 |  $-$364.383
$E_{\rho^{2+\alpha}}$ | 830.487 | 829.602 | 829.127 | 828.920 | 828.701 | 828.944
$E_{\rho\nabla^{2}\rho}$ | 60.942 | 60.880 | 60.825 | 60.819 | 60.805 | 60.809
$E_{\rm SO}$ |  $-$0.931 |  $-$0.942 |  $-$0.942 |  $-$0.946 |  $-$0.949 |  $-$0.940
$E_{\rm Coul}^{\rm Dir.}$ | 16.418 | 16.410 | 16.406 | 16.403 | 16.400 | 16.404
$E_{\rm Coul}^{\rm Exc.}$ |  $-$2.818 |  $-$2.817 |  $-$2.816 |  $-$2.816 |  $-$2.816 |  $-$2.817
$Q_{20}$ | 0.0238 | 0.0417 | 0.0315 | 0.0713 | 0.0889 | 0.0000
$\beta_{2}$ | 0.0002 | 0.0004 | 0.0003 | 0.0006 | 0.0008 | 0.0000

Table 3 lists the results for the heaviest doubly magic nucleus 208Pb using
the HFmix and the hfodd codes. For the HFmix calculations, the integration is
performed over the cubic box ($\approx[-12,+12]^{3}$ fm3). This model space is
large enough that the densities at the edges are smaller than 10-6 fm-3. The
results are shown for $dz=0.9$ and $0.833$ fm with $N_{\rm max}=12$ and 14.
One can see that the HFmix calculation with $dz=0.9$ fm provides reasonably
converged results, with the energy differences between the two grid spacings
being $\approx 50$ keV. The total energy calculated with the HFmix is
approximately 3 MeV (1 MeV) more bound for $N_{\rm max}=12$ $(14)$ than that
obtained with the hfodd. This is due to the limited number of HO bases used.
It can be expected that both codes would give closer results when $N_{\rm
max}$ is increased, which is similar to the case of 16O in Table 2.

Table 3: Similar to Table 2, but for 208Pb.  | HFmix($N_{\rm max}=12$) | hfodd | HFmix($N_{\rm max}=14$) | hfodd
---|---|---|---|---
|  $dz=0.9$ fm |  $dz=0.833$ fm | $N_{\rm max}=12$ |  $dz=0.9$ fm |  $dz=0.833$ fm | $N_{\rm max}=14$
| $N_{z}=26$ | $N_{z}=30$ |  | $N_{z}=26$ | $N_{z}=30$ |
$E_{\rm Total}$ |  $-$1633.868 |  $-$1633.816 |  $-$1630.639 |  $-$1635.052 |  $-$1634.979 |  $-$1634.148
$E_{\rm kin+c.m.}$ | 3864.117 | 3863.738 | 3861.319 | 3864.062 | 3863.647 | 3860.845
$E_{\rho^{2}}$ |  $-$22378.402 |  $-$22375.298 |  $-$22359.370 |  $-$22376.794 |  $-$22373.432 |  $-$22359.968
$E_{\rho\tau}$ | 1330.876 | 1330.586 | 1331.301 | 1329.914 | 1329.599 | 1328.865
$E_{\rho^{2+\alpha}}$ | 14534.801 | 14532.581 | 14523.995 | 14532.424 | 14530.027 | 14521.523
$E_{\rho\nabla^{2}\rho}$ | 314.241 | 314.176 | 311.518 | 314.988 | 314.921 | 314.602
$E_{\rm SO}$ |  $-$96.053 |  $-$96.126 |  $-$95.774 |  $-$96.174 |  $-$96.247 |  $-$96.375
$E_{\rm Coul}^{\rm Dir.}$ | 827.807 | 827.783 | 827.606 | 827.789 | 827.765 | 827.608
$E_{\rm Coul}^{\rm Exc.}$ |  $-$31.256 |  $-$31.255 |  $-$31.233 |  $-$31.260 |  $-$31.259 |  $-$31.248
$Q_{20}$ | 6 | 8 | 0 |  $-$3 |  $-$2 | 0
$\beta_{2}$ | 0.001 | 0.001 | 0.000 | 0.000 | 0.000 | 0.000

Figure 1 plots the neutron density profiles along the three Cartesian
directions for the cross sections that are closest to the center of the 208Pb
nucleus ($dz=0.45$ fm, $N_{\rm max}=12$). For instance, the $x$ profile
denotes the neutron density with $y=z=0.45$ fm ($dz=0.9$ fm), which is the
closest point to the center of the nucleus included in the discretization. For
this heavy nucleus, one notices rather small difference (The largest
difference appears at $z=1.35$ fm, $\Delta\rho\approx 0.0004$ fm-3.) between
the density profile along $z$ direction and those along the $x$ and $y$
directions. The difference is probably because the density is still converging
with increasing $N_{\rm max}$ in the $x$ and $y$ directions.

Figure 1: The neutron density profiles along $x$, $y$, and $z$ directions for
208Pb. The density curves of $\rho_{n}(x,a,a)$, $\rho_{n}(a,y,a)$, and the
discrete points $\rho_{n}(a,a,z)$ are shown, where $a=0.45$ fm.

In Figure 2, we plot the single-neutron energies for 208Pb using HFmix
($dz=0.9$ fm, $N_{\rm max}=12$) and compare them with the results using the
hfodd code ($N_{\rm max}=12$). We notice good agreement between the two codes.
For the HFmix result we see small energy splittings for the single-neutron
energies due to the breaking of the spherical symmetry in the basis, whereas
the spherical results of hfodd show perfect degeneracy. For the $0h_{9/2}$
levels, the energy splitting in HFmix is about 0.03 MeV. Even in the zoomed-in
plot shown in the inset of Fig. 2, the widths are negligible compared to the
large gaps which are of a few MeV.

Figure 2: The single-neutron energies for 208Pb calculated with HFmix and
hfodd codes.

####  III.1.3 Prolate nucleus: 24Mg

In Table 4, we present a set of calculations for a light deformed nucleus 24Mg
performed for fixed $(N_{\rm max},dz,N_{z}$) = (13, 0.833 fm, 22) and (13, 0.9
fm, 18), with $z$ and $x$ being the symmetry axes. It can be seen that the
results for both variants agree on the level of a few tens of keV ($\approx
30$ keV). Decreasing $dz$ to 0.833 fm results in the absolute binding energies
a few tens of keV smaller than those corresponding to $dz=0.9$ fm. Comparing
the HFmix results ($N_{\rm max},dz$)=(13, 0.833 fm) with those of hfodd, we
see the agreement is on the order of tens of keV. This better agreement
compared to 208Pb might be because the energies have been saturated to the
basis-number parameters, $N_{\rm max}$ and $dz$, for this lighter nucleus. It
is satisfying to note that the two codes agree on the calculated $\beta_{2}$
values up to the first two significant digits.

Table 4: The calculated energies (in MeV), quadrupole moments (in fm2), and deformation parameters for 24Mg using the HFmix and hfodd codes. We include results with the symmetry axis of the nucleus being aligned along the $z$ and $x$ directions. |  $N_{\rm max}=13$, $dz=0.9$ fm |  $N_{\rm max}=13$, $dz=0.833$ fm | hfodd
---|---|---|---
| $N_{z}=18$ | $N_{z}=22$ | $N_{\rm max}=13$
| along $z$ | along $x$ | along $z$ | along $x$ |
$E_{\rm Total}$ |  $-$195.732 |  $-$195.698 |  $-$195.670 |  $-$195.701 |  $-$195.660
$E_{\rm kin+c.m.}$ | 380.5068 | 380.539 | 380.757 | 380.582 | 380.462
$E_{\rho^{2}}$ |  $-$2164.306 |  $-$2164.325 |  $-$2163.843 |  $-$2163.843 |  $-$2161.812
$E_{\rho\tau}$ | 97.389 | 97.397 | 97.374 | 97.352 | 97.195
$E_{\rho^{2+\alpha}}$ | 1395.482 | 1395.493 | 1395.401 | 1395.109 | 1393.607
$E_{\rho\nabla^{2}\rho}$ | 89.039 | 89.008 | 88.960 | 88.946 | 88.804
$E_{\rm SO}$ |  $-$22.369 |  $-$22.335 |  $-$22.418 |  $-$22.370 |  $-$22.389
$E_{\rm Coul}^{\rm Dir.}$ | 32.913 | 32.914 | 32.919 | 32.912 | 32.903
$E_{\rm Coul}^{\rm Exc.}$ |  $-$4.388 |  $-$4.388 |  $-$4.389 |  $-$4.388 |  $-$4.387
$Q_{20}$; $Q_{22}$ | 112; 0 |  $-$56; 97 | 111; 0 |  $-$56; 96 | 112; 0
$\beta_{2}$; $\gamma$ | 0.515; 0∘ | 0.515; 120∘ | 0.510; 0∘ | 0.511; 120∘ | 0.515; 0∘

####  III.1.4 Triaxial nucleus: 64Ge

Finally, in Table 5 we list the calculated energies and quadrupole moments for
the triaxially deformed ground states of 64Ge, using the HFmix and hfodd
codes. For both variants, we use $N_{\rm max}=10$ and 12\. For fixed $N_{\rm
max}$, decreasing $dz$ from 0.9 fm to 0.83 fm results in a $\approx 40$ keV
less absolute binding energy. For $N_{\rm max}=10$, the total energy of HFmix
is $\approx 1.2$ MeV more bound compared to that of hfodd ($N_{\rm max}=10$).
Increasing $N_{\rm max}$ to 12 results in the difference of total energies
between HFmix and hfodd being only $\approx 0.5$ MeV. The difference is
expected to decrease with increasing $N_{\rm max}$. Again, the calculated
$\beta_{2}$ and $\gamma$ for this medium-heavy triaxial nucleus, are almost
identical for all the variants of the calculations.

Table 5: Similar to Table 3, but for 64Ge. | HFmix($N_{\rm max}=10$) | hfodd | HFmix($N_{\rm max}=12$) | hfodd
---|---|---|---|---
|  $dz=0.9$ fm |  $dz=0.83$ fm | $N_{\rm max}=10$ |  $dz=0.9$ fm |  $dz=0.83$ fm | $N_{\rm max}=12$
| $N_{z}=22$ | $N_{z}=26$ |  | $N_{z}=22$ | $N_{z}=26$ |
$E_{\rm Total}$ |  $-$542.795 |  $-$542.756 |  $-$541.610 |  $-$543.086 |  $-$543.045 |  $-$542.564
$E_{\rm kin+c.m.}$ | 1107.085 | 1106.908 | 1107.583 | 1107.070 | 1106.877 | 1107.151
$E_{\rho^{2}}$ |  $-$6571.053 |  $-$6569.677 |  $-$6568.314 |  $-$6569.662 |  $-$6568.087 |  $-$6570.121
$E_{\rho\tau}$ | 339.740 | 339.635 | 339.847 | 339.471 | 339.354 | 339.788
$E_{\rho^{2+\alpha}}$ | 4296.263 | 4295.267 | 4294.664 | 4294.943 | 4293.799 | 4295.701
$E_{\rho\nabla^{2}\rho}$ | 159.256 | 159.221 | 158.793 | 159.143 | 159.089 | 159.172
$E_{\rm SO}$ |  $-$39.230 |  $-$39.234 |  $-$39.387 |  $-$39.176 |  $-$39.182 |  $-$39.411
$E_{\rm Coul}^{\rm Dir.}$ | 177.485 | 177.465 | 177.550 | 177.465 | 177.444 | 177.499
$E_{\rm Coul}^{\rm Exc.}$ |  $-$12.342 |  $-$12.340 |  $-$12.347 |  $-$12.341 |  $-$12.339 |  $-$12.342
$Q_{20}$; $Q_{22}$ | 259; 135 | 259; 135 | 258; 135 | 258; 136 | 258; 136 | 255; 135
$\beta_{2}$; $\gamma$ | 0.262; 28∘ | 0.262; 28∘ | 0.261; 28∘ | 0.261; 28∘ | 0.261; 28∘ | 0.259; 28∘

###  III.2 HFB Results

In this section, we present the calculated HFB results for one spherical
nucleus (120Sn), one prolately deformed nucleus (34Mg), one triaxially
deformed nucleus (110Mo) using the HFBmix code. We compare them with those
calculated using the hfodd code. In addition, we include HFBmix results for
the typical configurations on the fission path of 240Pu.

####  III.2.1 Spherical nucleus: 120Sn

Table 6 lists the results of the HFB calculations for 120Sn using the HFBmix
and hfodd codes. The same $N_{\rm max}$ is used for both calculations. The
neutron pairing strength $V_{0}^{n}$ is set to $-200.0$ MeV fm3. In hfodd, the
results with a reduced neutron pairing strength that provides a similar
neutron pairing energy to HFBmix are also shown. The total energy for 120Sn of
HFBmix is $\approx 2$ MeV lower than that of the hfodd calculation. This may
be due to the relatively small $N_{\rm max}$ (and large $dz$) used and the
different ways of discretizing the continuum for HFBmix and hfodd.

Table 6: The calculated energies, Fermi energies, and pairing gaps (in MeV), quadrupole moments (in fm2), and $\beta_{2}$ for 120Sn using the HFBmix and hfodd codes. The Skyrme SLy4 force is used. The neutron pairing strengths are also listed in the first row, which are in MeV fm3. | HFBmix($N_{\rm max}=10$) | hfodd
---|---|---
|  $dz=1.0$ fm |  $dz=0.9$ fm | $N_{\rm max}=10$
| $N_{z}=22$ | $N_{z}=22$ |  |
$V_{0}^{n}$ |  $-$200.0 |  $-$200.0 |  $-$194.8 |  $-$200.0
$E_{\rm Total}$ |  $-$1017.506 |  $-$1017.398 |  $-$1015.458 |  $-$1015.818
$E_{\rm kin+c.m.}$ | 2174.025 | 2173.8967 | 2170.700 | 2174.099
$E_{\rho^{2}}$ |  $-$12683.665 |  $-$12680.942 |  $-$12655.007 |  $-$12661.993
$E_{\rho\tau}$ | 717.670 | 717.416 | 716.451 | 717.525
$E_{\rho^{2+\alpha}}$ | 8264.085 | 8261.9656 | 8244.251 | 8249.274
$E_{\rho\nabla^{2}\rho}$ | 224.724 | 224.619 | 222.092 | 222.655
$E_{\rm SO}$ |  $-$49.251 |  $-$49.284 |  $-$48.554 |  $-$49.914
$E_{\rm Coul}^{\rm Dir.}$ | 366.623 | 366.592 | 366.375 | 366.414
$E_{\rm Coul}^{\rm Exc.}$ |  $-$19.104 |  $-$19.103 |  $-$19.083 |  $-$19.085
$E^{\rm n}_{\rm pair}$ |  $-$12.614 |  $-$12.559 |  $-$12.682 |  $-$14.793
$\lambda_{\rm n}$ |  $-$7.974 |  $-$7.972 |  $-$7.959 |  $-$7.951
$\Delta_{\rm n}$ | 1.436 | 1.433 | 1.434 | 1.572
$Q_{20}$ | 7 |  $-$1 | 0 | 0
$\beta_{2}$ | 0.002 |  $-$0.0003 | 0.0 | 0.0
Figure 3: The calculated occupation probabilities $v^{2}_{n}$’s against the
quasi-particle energies for 120Sn.

Figure 3 plots the occupation probabilities

$v_{k,q}^{2}=\sum_{\sigma=\pm\frac{1}{2}}\int
d^{3}r|v_{k}^{q}(\mbox{\boldmath$r$},\sigma)|^{2},$ (27)

as a function of the quasi-particle energies for 120Sn, calculated with HFBmix
($dz=1.0$ fm) and hfodd codes ($V_{0}^{n}=-194.8$ MeV fm3). We see a nice
agreement between the two codes, especially for the occupied deep-hole states.
A small deviation of $E_{k}$ and $v^{2}_{k}$ from a perfect degeneracy of the
HFBmix calculation can be seen, which is associated with the breaking of the
spherical symmetry with the small quadrupole deformation, as can be assessed
from Table 6. For those states with $v^{2}_{n}$’s close to zero, we see the
HFBmix gives much more states inside the $E_{k}<60$ MeV window than the hfodd
does. This is expected because the two codes discretize the continua for
$E_{k}>-\lambda$ differently.

####  III.2.2 Prolate nucleus: 34Mg

In Table 7, we list the results of the HFB calculations with HFBmix ($N_{\rm
max}=14$, $dz=1.0$ fm) for 34Mg and compare the results with those of the
hfodd code. For the hfodd results we include one column with a deformed HO
basis using $(N_{\rm max}^{x,y},N_{\rm max}^{z})=(13,16)$. The last column
shows the hfodd results with spherical HO basis $N_{\rm max}=14$. For the two
hfodd results, the pairing strengths have been adjusted to match the $E_{\rm
pair}^{\rm{n,p}}$ of HFBmix calculations.

We see that the total energy calculated using HFBmix are lower by $<100$ keV
compared to those of hfodd. The quadrupole moments are rather close to each
other between the three variants, with the differences being $\approx 3\%$. It
has to be noted that the ground-state quadrupole moment of this nucleus
depends on the pairing strengths sensitively. The good agreement between the
three variants is thus rather satisfactory. As the pairing strengths of the
hfodd results are readjusted to give similar $E^{\rm n,p}_{\rm pair}$, the
extracted averaged pairing gaps are also close between the two codes.

Table 7: The calculated energies (in MeV), pairing gaps (in MeV), and quadrupole moments (in fm2) together with the pairing strengths (in MeV fm3) for 34Mg using HFBmix and hfodd codes. | HFBmix ($N_{\rm max}=14$) | hfodd
---|---|---
|  $dz=1.0$ fm, $N_{z}=18$ | $(N_{\rm max}^{x,y},N_{\rm max}^{z})=(13,16)$ | $N_{\rm max}^{x,y,z}=14$
$V_{0}^{n}$; $V_{0}^{p}$ |  $-$218.5; $-$218.5 |  $-$217.5; $-$204.0 |  $-$218.5; $-$205.0
$E_{\rm Total}$ |  $-$257.531 |  $-$257.471 |  $-$257.479
$E_{\rm kin+c.m.}$ | 557.036 | 557.070 | 556.983
$E_{\rm Skyrme}$ |  $-$835.730 |  $-$835.681 |  $-$835.614
$E_{\rm Coul}^{\rm Dir.}$ | 31.752 | 31.733 | 31.733
$E_{\rm Coul}^{\rm Exc.}$ |  $-$4.217 |  $-$4.217 |  $-$4.217
$E^{\rm n}_{\rm pair}$ |  $-$3.870 |  $-$3.905 |  $-$3.883
$E^{\rm p}_{\rm pair}$ |  $-$2.501 |  $-$2.472 |  $-$2.480
$\lambda_{\rm n}$ |  $-$3.221 |  $-$3.258 |  $-$3.260
$\lambda_{\rm p}$ |  $-$20.126 |  $-$20.121 |  $-$20.112
$\Delta_{\rm n}$ | 1.325 | 1.328 | 1.328
$\Delta_{\rm p}$ | 1.262 | 1.213 | 1.215
$Q_{20}$ | 99 | 102 | 102
$\beta_{2}$ | 0.255 | 0.262 | 0.262

####  III.2.3 Triaxial nucleus: 110Mo

Table 8 lists the calculated ground states of 110Mo using the HFBmix and hfodd
codes ($N_{\rm max}=9$) with SkM* force [36]. For HFBmix, we included results
with two orientations where the longest axes of the nucleus coincide with the
$z$ ($z>x>y$) and $x$ ($x>y>z$) axes. For these HFBmix runs, the dimension of
the HFB matrix [Eq. (8)] is 4840. Each iteration takes about 20 min to finish
on Intel(R) Xeon(R) CPU E5-2697 v4 processor.

The total energy obtained with $z>x>y$ is about 0.6 MeV lower than the result
with the longest axis in the $x$ direction. This is because the $z$ direction
is treated with the FD method allowing for better description of the spatially
extended density compared to the HO basis. It is very interesting to note that
the $\beta_{2}$ values of all the calculated triaxial ground states differ
only from the third significant digits.

Table 8: The calculated energies, pairing gaps (in MeV), and total quadrupole moments (in fm2) of 110Mo calculated with HFBmix and hfodd codes. The Skyrme force used is SkM*. The neutron pairing strengths are also shown in unit of MeV fm3, and for the proton pairing strength we use $V_{0}^{p}=-170.0$ MeV fm3. The calculated nucleus is placed in the principal axis, where $z>x>y$ indicates that the longest, medium, and shortest axes of the nucleus are aligned in the $z$, $x$, and $y$ directions, respectively. The $\beta_{2}$ and $\gamma$ deformations are defined in Sec II.2. | HFBmix($N_{\rm max}=9$, $dz=1.0$ fm, $N_{z}=22$) | hfodd
---|---|---
| $z>x>y$ | $x>y>z$ | $N_{\rm max}=9$
$V_{0}^{n}$ |  $-$170.0 |  $-$170.0 |  $-$172.7
$E_{\rm Total}$ |  $-$920.571 |  $-$919.965 |  $-$918.532
$E_{\rm kin+c.m.}$ | 2010.785 | 2011.839 | 2008.9216
$E_{\rho^{2}}$ |  $-$12268.080 |  $-$12265.004 |  $-$12235.894
$E_{\rho\tau}$ | 368.900 | 368.671 | 367.849
$E_{\rho^{2+\alpha}}$ | 8600.508 | 8597.598 | 8575.708
$E_{\rho\nabla^{2}\rho}$ | 193.730 | 192.776 | 190.496
$E_{\rm SO}$ |  $-$73.691 |  $-$73.000 |  $-$72.646
$E_{\rm Coul}^{\rm Dir.}$ | 266.902 | 266.847 | 266.618
$E_{\rm Coul}^{\rm Exc.}$ |  $-$15.749 |  $-$15.740 |  $-$15.718
$E_{\rm pair}^{\rm n}$ |  $-$3.875 |  $-$3.885 |  $-$3.867
$\lambda_{\rm n}$ |  $-$5.414 |  $-$5.407 |  $-$5.379
$\Delta_{\rm n}$ | 0.762 | 0.763 | 0.771
$Q_{20}$; $Q_{22}$ | 949; 361 |  $-$780; 632 | 943; 365
$\beta_{2}$; $\gamma$ | 0.369; 21∘ | 0.364; $-$39∘ | 0.367; 21∘

####  III.2.4 Selected Points on the Fission pathway of 240Pu

Finally, in Table 9, we present the calculated results for the fission pathway
of 240Pu using HFBmix (SkM*). For this calculation, we use $(N_{\rm
max},dz,N_{z})$=(8, 1.2 fm, 26). The neutron and proton pairing strengths are
$V_{0}^{n}=V_{0}^{p}=-170.0$ MeV fm3. These pairing strengths give reasonable
pairing energy for the ground state of 240Pu. We select several configurations
along the fission pathway, that is, the results corresponding to the ground
state (g.s), the fission isomer (Fis. Iso.), the first/inner fission barrier
(1st Fis. Bar.), the second/outer fission barrier (2nd Fis. Bar.), and one
deformation point beyond the second fission barrier at $Q_{20}=150$ b. Here,
we have assumed that the fission pathway is described in the space spanned by
the several multipole deformations. Ideally, one should obtain the solutions
corresponding to the barriers by calculating the full potential energy surface
before the scission area. However, due to the limitation of the computational
resources, we restrict our search within the points near the two fission
barriers. For example, to find the first fission barrier, we perform
calculations with $Q_{20}$ being constrained from 50 to 70 b at a step of 5 b.
We do not impose axial nor reflection symmetries; we allow $Q_{22}$ and
$Q_{30}$ to vary freely when $Q_{20}$ is constrained. The important
deformations which lower the fission barrier and the information on the
barrier position in the potential energy surface in the case of SkM*
functional are known from the literature such as Ref. [37].

Table 9: The calculated energies (in MeV), pairing gaps (in MeV), quadrupole moments [in barn (b)], and octupole moment (in b3/2) for selected deformation points on the fission path of 240Pu using HFBmix (SkM* force). | g.s. | Fis. Iso. | 1st Fis. Bar. | 2nd Fis. Bar. |
---|---|---|---|---|---
$E_{\rm Total}$ |  $-$1794.741 |  $-$1792.300 |  $-$1787.055 |  $-$1788.195 |  $-$1790.132
$E_{\rm kin+c.m.}$ | 4447.739 | 4441.681 | 4442.885 | 4427.123 | 4431.331
$E_{\rm Skyrme}$ |  $-$7230.115 |  $-$7186.925 |  $-$7192.441 |  $-$7135.954 |  $-$7118.399
$E_{\rm Coul}^{\rm Dir.}$ | 1025.625 | 992.332 | 1006.738 | 963.173 | 939.697
$E_{\rm Coul}^{\rm Exc.}$ |  $-$35.687 |  $-$35.603 |  $-$35.620 |  $-$35.537 |  $-$35.540
$E_{\rm pair}^{\rm n}$ |  $-$1.399 |  $-$2.937 |  $-$7.532 |  $-$6.579 |  $-$6.241
$E_{\rm pair}^{\rm p}$ |  $-$0.904 |  $-$0.848 |  $-$1.085 |  $-$0.652 |  $-$0.979
$\lambda_{\rm n}$ |  $-$6.383 |  $-$6.626 |  $-$6.257 |  $-$6.310 |  $-$6.428
$\lambda_{\rm p}$ |  $-$4.484 |  $-$4.715 |  $-$4.574 |  $-$5.155 |  $-$5.180
$\Delta_{\rm n}$ | 0.318 | 0.463 | 0.747 | 0.687 | 0.665
$\Delta_{\rm p}$ | 0.264 | 0.262 | 0.297 | 0.228 | 0.251
$Q_{20}$ | 29.9 | 83.2 | 60.6 | 120.0 | 150.3
$Q_{22}$ | 0 | 0 | 8.6 | 4.9 | 3.0
$Q_{30}$ | 0 | 0 | 0 | 6.1 | 9.4

The excitation energy of the fission isomer extracted from Table 9 is 2.4 MeV;
the heights of the inner and outer fission barriers are 7.7 and 6.5 MeV,
respectively. Comparing these three energies with those of Ref. [37]
calculated with hfodd, we notice that our results are smaller than theirs by
$\leq 0.3$ MeV. These small deviations can be attributed to the rather small
$N_{\rm max}=8$ used in our calculations. Another source of discrepancies is
the differences in the pairing treatment (both pairing strengths and pairing
types chosen).

Figure 4: The density profiles of the fission isomer of 240Pu calculated using
HFBmix (Table 9). The density values correspond to the lines with other
coordinates closest to the center of the nucleus.

Figure 4 shows the density profiles of the fission isomer along the three
axes. In Table 9 we include the calculated results for a deformation beyond
the second fission barrier. In the future, the HFBmix code will be extended to
allow for the full time-dependent HFB simulation which will be started from a
HFB solution corresponding to the deformation point beyond the second fission
barrier.

Although the HFBmix is capable of describing the extremely elongated
configuration, such as those near or beyond the scission point of the
fissioning nucleus, it is not the intention of the current paper to perform
the full fission-path calculation of 240Pu. Here we are content to demonstrate
the feasibility of the HFBmix code for describing some of the typical
configurations of 240Pu fission pathway. We leave a more realistic and
systematic survey of the potential energy surface for the future project with
larger basis size and better pairing treatment.

##  IV Summary and Perspective

We solve the numerically demanding three-dimensional (3D) coordinate-space
nuclear Hartree-Fock-Bogoliubov (HFB) problem using a mixture of bases, which
consist of two harmonic-oscillator (HO) bases in the $x$ and $y$ directions,
and one finite-difference basis in the $z$ direction.

We implement the 3D HO and the mixed-basis codes. Using the mixed-basis code,
we perform systematic calculations for spherical and deformed nuclei using the
Hartree-Fock (HF) and HFB methods. The results obtained are compared with
those of the 3D HO-based code, hfodd.

The differences of the total energies using the HF method for lightest nucleus
16O is in the order of a few tens of keV between the two HO-based codes, HO3D
and hfodd. For the lightest nuclei, 16O and 24Mg, where the total energies
tend to be converged with increasing $N_{\rm max}$, the total energy
differences are only a few tens of keV. For heavier nuclei ($A>100$), the
differences in total energies are about 1-2 MeV. The differences decrease with
the increase of the largest $N_{\rm max}$ of the HO bases. Hence, for the
heavier nuclei, the difference between mixed-basis code and hfodd can be
attributed to the fact that the HF result is still converging with $N_{\rm
max}$.

Comparisons of the HFB results using the HFBmix and the hfodd codes show
similar convergence behavior with respect to $N_{\rm max}$ as the HF case.
Specifically, the differences of the total energies using HFBmix and hfodd
decrease with increasing $N_{\rm max}$.

For spherical nuclei, the calculated quadrupole moments are negligibly small
for the mixed-basis HF and HFB codes. For deformed nuclei, the differences in
the quadrupole moments calculated using the mixed-basis codes are almost
identical to those of the hfodd code. A good agreement of single/quasi-
particle energies is also seen between the mixed-basis codes and that of the
hfodd code.

As a demonstration of the applicability of the HFBmix code to extremely
elongated deformations, selected configurations along the fission pathway in
240Pu are calculated. With future optimizations of the code, the HFBmix will
be advantageous in a systematic survey of the potential energy surfaces with
elongated configurations as well as the deformed drip-line nuclei.

Our future work will be focused on the time-dependent extension of the current
mixed-basis HFB method. Such a mixed-basis strategy would make the full 3D
time-dependent HFB calculations much more affordable and flexible.

###### Acknowledgements.

The current work is supported by National Natural Science Foundation of China
(Grant No. 12075068 and No. 11705038), the Fundamental Research Funds for the
Central Universities (Grant No. HIT.BRET.2021003), JSPS KAKENHI Grant No.
16K17680 and No. 20K03964, and the JSPS-NSFC Bilateral Program for the Joint
Research Project on “Nuclear mass and life for unravelling mysteries of
r-process.” YS thanks the HPC Studio at Physics Department of Harbin Institute
of Technology for computing resources allocated through<EMAIL_ADDRESS>
##  Appendix A The densities in the mixed basis

In this section, we specifically provide the detailed formulae of various
densities [Eq. (II.1)] using the mixed basis.

To compute the densities, $\tau$, $\divergence{\mathbf{J}}$, and
$\gradient{\rho}$, which are needed for the construction of the potentials and
mean-field Hamiltonian, one needs the following first derivatives

$\displaystyle\partialderivative{x}v_{k}(x_{i_{x}},y_{i_{y}},z_{i_{z}})$
$\displaystyle=$
$\displaystyle\sum_{n_{x}}\sum_{n_{y}}\psi^{\prime}_{n_{x}}(x_{i_{x}})\psi_{n_{y}}(y_{i_{y}})\,v^{n_{x},n_{y},i_{z}}_{k},$
(28) $\displaystyle\partialderivative{y}v_{k}(x_{i_{x}},y_{i_{y}},z_{i_{z}})$
$\displaystyle=$
$\displaystyle\sum_{n_{x}}\sum_{n_{y}}\psi_{n_{x}}(x_{i_{x}})\psi^{\prime}_{n_{y}}(y_{i_{y}})\,v^{n_{x},n_{y},i_{z}}_{k},$
(29) $\displaystyle\partialderivative{z}v_{k}(x_{i_{x}},y_{i_{y}},z_{i_{z}})$
$\displaystyle=$
$\displaystyle\sum_{n_{x}}\sum_{n_{y}}\psi_{n_{x}}(x_{i_{x}})\psi_{n_{y}}(y_{i_{y}})\sum_{d=-\lambda}^{+\lambda}m^{(1)}_{d}v_{k}^{n_{x},n_{y},i_{z}+d},$
(30)

where the $m^{(1)}_{d}$ is the FD coefficients for the first derivative. For
the nine-point FD case ($\lambda=4$), the coefficients for the first
($m^{(1)}_{d}$) and second ($m^{(2)}_{d}$) derivatives are listed in Table 10.
These are used throughout this work.

Table 10: The nine-point FD weights of the first and second derivatives. $d$ |  $-$4 |  $-$3 |  $-$2 |  $-$1 | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---|---|---|---|---
$m^{(1)}_{d}$ | $\frac{1}{280}$ | $-\frac{4}{105}$ | $\frac{21}{105}$ | $-\frac{84}{105}$ | 0 | $\frac{84}{105}$ | $-\frac{21}{105}$ | $\frac{4}{105}$ | $-\frac{1}{280}$
$m^{(2)}_{d}$ | $-\frac{1}{560}$ | $\frac{8}{315}$ | $-\frac{1}{5}$ | $\frac{8}{5}$ | $-\frac{205}{72}$ | $\frac{8}{5}$ | $-\frac{1}{5}$ | $\frac{8}{315}$ | $-\frac{1}{560}$

As for the term involving Laplacian ($\laplacian{\rho}$), one needs
$\laplacian{v_{k}}$. This can be done in a similar manner as Eqs. (28) - (30),

$\displaystyle\partialderivative[2]{x}v_{k}(x_{i_{x}},y_{i_{y}},z_{i_{z}})$
$\displaystyle=$
$\displaystyle\sum_{n_{x}}\sum_{n_{y}}\psi^{\prime\prime}_{n_{x}}(x_{i_{x}})\psi_{n_{y}}(y_{i_{y}})\,v^{n_{x},n_{y},i_{z}}_{k},$
(31)
$\displaystyle\partialderivative[2]{y}v_{k}(x_{i_{x}},y_{i_{y}},z_{i_{z}})$
$\displaystyle=$
$\displaystyle\sum_{n_{x}}\sum_{n_{y}}\psi_{n_{x}}(x_{i_{x}})\psi^{\prime\prime}_{n_{y}}(y_{i_{y}})\,v^{n_{x},n_{y},i_{z}}_{k},$
(32)
$\displaystyle\partialderivative[2]{z}v_{k}(x_{i_{x}},y_{i_{y}},z_{i_{z}})$
$\displaystyle=$
$\displaystyle\sum_{n_{x}}\sum_{n_{y}}\psi_{n_{x}}(x_{i_{x}})\psi_{n_{y}}(y_{i_{y}})\sum_{d=-\lambda}^{+\lambda}m^{(2)}_{d}v_{k}^{n_{x},n_{y},i_{z}+d}.$
(33)

In the case of 3D HO basis, the calculation of gradient and Laplacian of
$v_{k}$ is more straightforward where the $z$ direction is calculated in a
similar way as those of the $x$ and $y$ directions.

After obtaining the gradient and the Laplacian of $v_{k}$, one assembles the
densities [38]

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle
D^{\uparrow\uparrow}_{00}+D^{\downarrow\downarrow}_{00},$ $\displaystyle\tau$
$\displaystyle=$
$\displaystyle\sum_{\mu}(D^{\uparrow\uparrow}_{\mu\mu}+D^{\downarrow\downarrow}_{\mu\mu}),$
$\displaystyle\laplacian{\rho}$ $\displaystyle=$ $\displaystyle
2(L^{\uparrow\uparrow}+L^{\downarrow\downarrow})+2\tau,$
$\displaystyle\divergence{\mathbf{J}}$ $\displaystyle=$
$\displaystyle-2\Im(D^{\uparrow\downarrow}_{23}-D^{\uparrow\downarrow}_{32})$
$\displaystyle-2(D^{\uparrow\downarrow}_{31}-D^{\uparrow\downarrow}_{13})$
$\displaystyle+2\Im(D^{\uparrow\uparrow}_{21}-D^{\downarrow\downarrow}_{21}),$
$\displaystyle\nabla_{k}{\rho}$ $\displaystyle=$ $\displaystyle
2(D^{\uparrow\uparrow}_{k0}+D^{\downarrow\downarrow}_{k0}),~{}~{}k=1,2,3,$

where

$\displaystyle D^{\sigma\sigma^{\prime}}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\sum_{k}[\nabla_{\mu}v_{k}(\mbox{\boldmath$r$}\sigma)][\nabla_{\nu}v^{*}_{k}(\mbox{\boldmath$r$}\sigma^{\prime})],$
(34) $\displaystyle L^{\sigma\sigma^{\prime}}$ $\displaystyle=$
$\displaystyle\sum_{k}v_{k}(\mbox{\boldmath$r$}\sigma)[\laplacian
v^{*}_{k}(\mbox{\boldmath$r$}\sigma^{\prime})],$ (35)

and $\nabla_{\mu}\equiv(1,\bm{\nabla})$. Using Eqs. (28) - (33), one can also
obtain the Laplacian of the effective mass in Eq. (36).

##  Appendix B The kinetic terms in HO and mixed bases

In this subsection, we give the matrix elements of the kinetic term in the HO
and mixed-basis representations. Before showing the kinetic term in the
coordinate representation, one notices that, when applied to a function, the
kinetic part of Eq. (9) can be written as [16]

$\displaystyle\hat{h}^{\rm
kin.}\psi(\mbox{\boldmath$r$})=-\frac{1}{2}\Bigg{\\{}\frac{\hbar^{2}}{2m^{*}(\mbox{\boldmath$r$})}\nabla^{2}\psi(\mbox{\boldmath$r$})+\nabla^{2}\quantity[\frac{\hbar^{2}}{2m^{*}(\mbox{\boldmath$r$})}\psi(\mbox{\boldmath$r$})]-\quantity(\nabla^{2}\frac{\hbar^{2}}{2m^{*}(\mbox{\boldmath$r$})})\psi(\mbox{\boldmath$r$})\Bigg{\\}}.$
(36)

The matrix element of the third term of right hand side of Eq. (36) is
diagonal in the FD basis.

###  B.1 The matrix element of kinetic terms in the 1D grid basis

Next we examine the matrix form of the Laplacian operator $\nabla^{2}$ in the
1D grid representation. The second derivative of a function in the 1D
coordinate space is expressed as

$\derivative[2]{f}{z}\Big{|}_{z_{i}}=\sum_{d=-\lambda}^{+\lambda}f(z_{i+d})\,m^{(2)}_{d}\times(dz)^{-2}.$
(37)

Using the matrix form of the Laplacian operator, together with the third local
term in Eq. (36), we have the kinetic matrix elements in the 1D FD
representation

$T_{ij}=-\frac{1}{2}[h^{\rm 1D}(z_{i})+h^{\rm
1D}(z_{j})]\,\nabla^{2}_{ij}-\left[\nabla^{2}h^{\rm
1D}(z)|_{z=z_{i}}\right]\delta_{ij},$ (38)

where $\nabla^{2}_{ij}=m_{d}^{(2)}\times(dz)^{-2}\delta_{j,i+d}$,
$(|d|\leq\lambda)$ denotes the matrix elements for the Laplacian operator in
the FD basis. The function $h^{\rm 1D}(z)$ is the 1D conterpart of
$\displaystyle\frac{\hbar^{2}}{2m^{*}(\mbox{\boldmath$r$})}$ in Eq. (36); $dz$
is the spacing in the $z$ direction.

###  B.2 Matrix elements in the 3D coordinate-space HO basis

The matrix elements in the 3D HO basis for the kinetic part (36) are in the
following form

$\displaystyle\langle n_{x}n_{y}n_{z}|\hat{h}^{\rm
kin.}|n^{\prime}_{x}n^{\prime}_{y}n^{\prime}_{z}\rangle$ $\displaystyle=$
$\displaystyle\int dx\int dy\int
dz\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\psi^{*}_{n_{z}}(z)\hat{h}\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)\psi_{n^{\prime}_{z}}(z)$
(39) $\displaystyle=$ $\displaystyle-\frac{1}{2}\int dx\int dy\int
dz\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\psi^{*}_{n_{z}}(z)\frac{\hbar^{2}}{2m^{*}(\mbox{\boldmath$r$})}\nabla^{2}\quantity[\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)\psi_{n^{\prime}_{z}}(z)]$
$\displaystyle-\frac{1}{2}\int dx\int dy\int
dz\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\psi^{*}_{n_{z}}(z)\nabla^{2}\quantity[\frac{\hbar^{2}}{2m^{*}(\mbox{\boldmath$r$})}\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)\psi_{n^{\prime}_{z}}(z)]$
$\displaystyle+\frac{1}{2}\int dx\int dy\int
dz\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\psi^{*}_{n_{z}}(z)\quantity[\nabla^{2}\frac{\hbar^{2}}{2m^{*}(\mbox{\boldmath$r$})}]\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)\psi_{n^{\prime}_{z}}(z),$

where $\psi_{n_{\mu}}(\mu)$ is the HO wave functions in the $\mu=x$, $y$, and
$z$. The form of the Cartesian-coordinate HO basis $|n_{x}n_{y}n_{z}\rangle$
is identical to those explained in Ref. [34].

###  B.3 Matrix elements in the mixed basis

We write down the matrix elements in the mixed basis before a brief
explanation of each term:

$\displaystyle\langle n_{x}n_{y}i_{z}|\hat{h}^{\rm
kin.}|n^{\prime}_{x}n^{\prime}_{y}i^{\prime}_{z}\rangle$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int dx\int
dy\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\quantity[\frac{\hbar^{2}}{2m^{*}(x,y,z_{i_{z}})}]\quantity(\frac{d^{2}}{dx^{2}}+\frac{d^{2}}{dy^{2}})[\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)]\delta_{i_{z},i^{\prime}_{z}}(dz)^{-2}$
(40) $\displaystyle-\frac{1}{2}\int dx\int
dy\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\quantity(\frac{d^{2}}{dx^{2}}+\frac{d^{2}}{dy^{2}})\quantity[\frac{\hbar^{2}}{2m^{*}(x,y,z_{i_{z}})}\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)]\delta_{i_{z},i^{\prime}_{z}}(dz)^{-2}$
$\displaystyle-\frac{1}{2}\int dx\int
dy\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\frac{1}{2}\quantity[\frac{\hbar^{2}}{2m^{*}(x,y,z_{i_{z}})}+\frac{\hbar^{2}}{2m^{*}(x,y,z_{i^{\prime}_{z}})}][\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)]m^{(2)}_{i_{z}-i^{\prime}_{z}}(dz)^{-2}$
$\displaystyle+\frac{1}{2}\int dx\int
dy\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\quantity[\nabla^{2}\frac{\hbar^{2}}{2m^{*}(\mbox{\boldmath$r$})}]\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)\delta_{i_{z},i^{\prime}_{z}}(dz)^{-2},$

where $|n_{x}n_{y}i_{z}\rangle$ is a 3D basis which uses the Cartesian
coordinate HO bases $\psi_{n_{x}}(x)$ and $\psi_{n_{y}}(y)$ in the $x$ and $y$
directions, respectively. For the $z$ direction, we use the FD basis, for
which the $\nabla^{2}$ and $\nabla$ operators have been shown in Sec. B.1.

The first two terms in the right-hand side of Eq. (40) correspond to the first
two terms in Eq. (36) for the $x$ and $y$ directions. In these two terms, the
matrix is diagonal in the $z$-direction, due to the local property. The third
term corresponds to the first two terms in Eq. (36) for the $z$ direction. It
has been symmetrized in the same way as explained in Sec. B.1. The
coefficients $m^{(2)}_{i_{z}-i^{\prime}_{z}}$ are non-zero only for
$|i_{z}-i^{\prime}_{z}|\leq\lambda$. The last term corresponds to the last
term in Eq. (36). It is diagonal for the basis index in the $z$ direction due
to its locality.

##  Appendix C The SO terms in various basis

For the SO term, the matrix form can be obtained in a very similar way. The SO
part of the Hamiltonian reads

$-i\begin{pmatrix}\displaystyle\partialderivative{f(\mbox{\boldmath$r$})}{x}\partialderivative{y}-\partialderivative{f(\mbox{\boldmath$r$})}{y}\partialderivative{x}&\displaystyle\quantity[\partialderivative{f(\mbox{\boldmath$r$})}{y}\partialderivative{z}-\partialderivative{f(\mbox{\boldmath$r$})}{z}\partialderivative{y}]-i\quantity[\partialderivative{f(\mbox{\boldmath$r$})}{z}\partialderivative{x}-\partialderivative{f(\mbox{\boldmath$r$})}{x}\partialderivative{z}]\\\
\displaystyle\quantity[\partialderivative{f(\mbox{\boldmath$r$})}{y}\partialderivative{z}-\partialderivative{f(\mbox{\boldmath$r$})}{z}\partialderivative{y}]+i\quantity[\partialderivative{f(\mbox{\boldmath$r$})}{z}\partialderivative{x}-\partialderivative{f(\mbox{\boldmath$r$})}{x}\partialderivative{z}]&\displaystyle\partialderivative{f(\mbox{\boldmath$r$})}{y}\partialderivative{x}-\partialderivative{f(\mbox{\boldmath$r$})}{x}\partialderivative{y}\end{pmatrix},$

where $f(\mbox{\boldmath$r$})\equiv
b_{4}\rho(\mbox{\boldmath$r$})-b^{\prime}_{4}\rho_{q}(\mbox{\boldmath$r$})$ is
a scalar function of position.

Next, we list the two typical terms in the SO term:
$\displaystyle\partialderivative{f(\mbox{\boldmath$r$})}{x}\partialderivative{y}$
and
$\displaystyle\partialderivative{f(\mbox{\boldmath$r$})}{y}\partialderivative{z}$,
which contain operators with and without derivative in the $z$ direction:

$\displaystyle\langle
n_{x}n_{y}i_{z}|\partialderivative{f(\mbox{\boldmath$r$})}{x}\partialderivative{y}|n^{\prime}_{x}n^{\prime}_{y}i^{\prime}_{z}\rangle$
$\displaystyle=$ $\displaystyle\int dx\int
dy\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\partialderivative{f(\mbox{\boldmath$r$})}{x}\psi_{n^{\prime}_{x}}(x)\derivative{\psi_{n^{\prime}_{y}}(y)}{y}\delta_{i_{z},i^{\prime}_{z}},$
(41) $\displaystyle\langle
n_{x}n_{y}i_{z}|\partialderivative{f(\mbox{\boldmath$r$})}{y}\partialderivative{z}|n^{\prime}_{x}n^{\prime}_{y}i^{\prime}_{z}\rangle$
$\displaystyle=$ $\displaystyle\int dx\int
dy\,\psi^{*}_{n_{x}}(x)\psi^{*}_{n_{y}}(y)\partialderivative{f(\mbox{\boldmath$r$})}{y}\psi_{n^{\prime}_{x}}(x)\psi_{n^{\prime}_{y}}(y)m^{(1)}_{i^{\prime}_{z}-i_{z}}(dz)^{-1}.$
(42)

## References

  * Nakatsukasa  _et al._ [2016] T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, and K. Yabana, Rev. Mod. Phys. 88, 045004 (2016).
  * Ring and Schuck [1980] P. Ring and P. Schuck,  _The Nuclear Many-Body Problem_ (Springer-Verlag, Berlin, 1980).
  * Brack  _et al._ [1972] M. Brack, J. Damgaard, A. S. Jensen, H. C. Pauli, V. M. Strutinsky, and C. Y. Wong, Rev. Mod. Phys. 44, 320 (1972).
  * Bender  _et al._ [2003] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * Vautherin and Brink [1972] D. Vautherin and D. M. Brink, Phys. Rev. C 5, 626 (1972).
  * Dobaczewski  _et al._ [1984] J. Dobaczewski, H. Flocard, and J. Treiner, Nucl. Phys. A 422, 103 (1984).
  * Gogny [1975] D. Gogny, Nucl. Phys. A 237, 399 (1975).
  * Dechargé and Gogny [1980] J. Dechargé and D. Gogny, Phys. Rev. C 21, 1568 (1980).
  * Younes and Gogny [2009] W. Younes and D. Gogny, Phys. Rev. C 80, 054313 (2009).
  * Serot and Walecka [1992] B. D. Serot and J. D. Walecka, Relativistic nuclear many-body theory, in  _Recent Progress in Many-Body Theories: Volume 3_ , edited by T. L. Ainsworth, C. E. Campbell, B. E. Clements, and E. Krotscheck (Springer US, 1992) pp. 49–92.
  * Ring [1996] P. Ring, Prog. Part. Nucl. Phy. 37, 193 (1996).
  * Meng  _et al._ [2006] J. Meng, H. Toki, S. Zhou, S. Zhang, W. Long, and L. Geng, Prog. Part. Nucl. Phys. 57, 470 (2006).
  * Shen  _et al._ [2019] S. Shen, H. Liang, W. H. Long, J. Meng, and P. Ring, Prog. Part. Nucl. Phys. 109, 103713 (2019).
  * Bonche  _et al._ [1985] P. Bonche, H. Flocard, P. Heenen, S. Krieger, and M. Weiss, Nucl. Phys. A 443, 39 (1985).
  * Pei  _et al._ [2014] J. C. Pei, G. I. Fann, R. J. Harrison, W. Nazarewicz, Y. Shi, and S. Thornton, Phys. Rev. C 90, 024317 (2014).
  * Jin  _et al._ [2017] S. Jin, A. Bulgac, K. Roche, and G. Wlazłowski, Phys. Rev. C 95, 044302 (2017).
  * Kashiwaba and Nakatsukasa [2019] Y. Kashiwaba and T. Nakatsukasa, Phys. Rev. C 100, 035804 (2019).
  * Yamagami [2019] M. Yamagami, Phys. Rev. C 100, 054302 (2019).
  * Bennaceur and Dobaczewski [2005] K. Bennaceur and J. Dobaczewski, Comput. Phys. Commun. 168, 96 (2005).
  * Pöschl  _et al._ [1997] W. Pöschl, D. Vretenar, and P. Ring, Comp. Phys. Commun. 103, 217 (1997).
  * Terán  _et al._ [2003] E. Terán, V. E. Oberacker, and A. S. Umar, Phys. Rev. C 67, 064314 (2003).
  * Pei  _et al._ [2008] J. C. Pei, M. V. Stoitsov, G. I. Fann, W. Nazarewicz, N. Schunck, and F. R. Xu, Phys. Rev. C 78, 064306 (2008).
  * Yoshida and Giai [2008] K. Yoshida and N. V. Giai, Phys. Rev. C 78, 064316 (2008).
  * Ravenhall  _et al._ [1983] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys. Rev. Lett. 50, 2066 (1983).
  * Hashimoto  _et al._ [1984] M. Hashimoto, H. Seki, and M. Yamada, Prog. Theor. Phys. 71, 320 (1984).
  * Magierski and Heenen [2002] P. Magierski and P.-H. Heenen, Phys. Rev. C 65, 045804 (2002).
  * Newton and Stone [2009] W. G. Newton and J. R. Stone, Phys. Rev. C 79, 055801 (2009).
  * Schuetrumpf and Nazarewicz [2015] B. Schuetrumpf and W. Nazarewicz, Phys. Rev. C 92, 045806 (2015).
  * Hashimoto [2012] Y. Hashimoto, Eur. Phys. J. A 48, 55 (2012).
  * Hashimoto [2013] Y. Hashimoto, Phys. Rev. C 88, 034307 (2013).
  * Dobaczewski  _et al._ [2009] J. Dobaczewski, W. Satuła, B. Carlsson, J. Engel, P. Olbratowski, P. Powałowski, M. Sadziak, J. Sarich, N. Schunck, A. Staszczak, M. Stoitsov, M. Zalewski, and H. Zduńczuk, Comput. Phys. Commun. 180, 2361 (2009).
  * Shi [2018] Y. Shi, Phys. Rev. C 98, 014329 (2018).
  * Ryssens  _et al._ [2015] W. Ryssens, V. Hellemans, M. Bender, and P.-H. Heenen, Comput. Phys. Commun. 187, 175 (2015).
  * Dobaczewski and Dudek [1997a] J. Dobaczewski and J. Dudek, Comput. Phys. Commun. 102, 166 (1997a).
  * Shi  _et al._ [2020] Y. Shi, N. Hinohara, and B. Schuetrumpf, Phys. Rev. C 102, 044325 (2020).
  * Bartel  _et al._ [1982] J. Bartel, P. Quentin, M. Brack, C. Guet, and H.-B. Håkansson, Nucl. Phys. A 386, 79 (1982).
  * Schunck  _et al._ [2014] N. Schunck, D. Duke, H. Carr, and A. Knoll, Phys. Rev. C 90, 054305 (2014).
  * Dobaczewski and Dudek [1997b] J. Dobaczewski and J. Dudek, Comput. Phys. Commun. 102, 166 (1997b).