Tests

In this section, we discuss results of two numerical tests. First, by switching off the Local Scaling Transformation (LST) of THO, we run HFBTHO in the axial HO basis and test it against HFODD. For a given Skyrme interaction and zero-range, density-dependent pairing force, both codes should give exactly the same results. Since technical details of the inner structure of both codes are completely different, such calculations constitute an extremely stringent test of both codes.

Second, by switching the LST on, we could test the code HFBTHO against the spherical code HFBRAD [5]. Here, results of both codes cannot be exactly identical, because the phase spaces in which the solutions are obtained are significantly different.

Table 1: (Color online) Benchmark results of the HFB calculations performed for the SLy4 interaction and mixed $\delta$ pairing. All energies are in MeV, lengths in fm, and quadrupole moments in barns. Boldface colored digits differ between the HFBTHO and HFODD/HFBRAD calculations. See text for details.

Nucleus:	$^{208}$Pb		$^{168}$Er		$^{168}$Er		$^{120}$Sn
Code:	HFBTHO	HFODD	HFBTHO	HFODD	HFBTHO	HFODD	HFBTHO	HFBRAD
Basis:	2D-HO	3D-HO	2D-HO	3D-HO	2D-HO	3D-HO	2D-THO	Radial
$N_0$	14	14	14	14	$N_\perp$=13, $N_z$=17	$N_\perp$=13, $N_z$=17	25	n.a.
$N_{{st}}$	680	680	680	680	680	680	3276	n.a.
$N^{{qp}}_n$	532	532	489	489	497	497	924	4260
$N^{{qp}}_p$	481	481	448	448	451	451	855	4003
$b_\perp$	2.2348121	2.2348121	2.1566616	2.1566616	2.0581218	2.0581218	2.0390141	n.a.
$b_z$	2.2348121	2.2348121	2.1566616	2.1566616	2.3681210	2.3681210	2.0390141	n.a.
$\lambda_n$	$-$8.114095	$-$8.114020	$-$6.936061	$-$6.936058	$-$6.943872	$-$6.943858	$-$8.016795	$-$8.018081
$\lambda_p$	$-$8.810501	$-$8.810445	$-$7.156485	$-$7.156477	$-$7.152114	$-$7.152007	$-$11.107284	$-$11.107777
$\Delta_n$	0	0	0.394570	0.394578	0.392326	0.392327	1.244750	1.244648
$\Delta_p$	0	0	0.390601	0.390605	0.397728	0.397746	0	0
$R_n$	5.619758	5.619757	5.357578	5.357578	5.360037	5.360044	4.730466	4.730184
$R_p$	5.460080	5.460090	5.225538	5.225539	5.227218	5.227231	4.593884	4.593653
$Q_n$	$-$0.000022	6.6E-11	11.473921	11.473920	11.567875	11.567983	$-$0.001055	0
$Q_p$	$-$0.000017	4.7E-11	7.880228	7.880224	7.930128	7.930227	$-$0.000631	0
$\epsilon ^{{gs}}_n$	$-$58.001139	$-$58.001145	$-$56.014966	$-$56.014973	$-$55.996356	$-$55.996370	$-$55.756516	$-$55.755837
$\epsilon ^{{gs}}_p$	$-$44.042810	$-$44.042814	$-$44.422148	$-$44.422167	$-$44.486154	$-$44.486271	$-$46.629670	$-$46.631739
$\Sigma^{\epsilon }_n$	$-$3009.265452	$-$3009.264720	$-$2401.023343	$-$2401.023305	$-$2401.701865	$-$2401.698888	$-$1667.965633	$-$1668.063705
$\Sigma^{\epsilon }_p$	$-$1678.791400	$-$1678.790238	$-$1439.480739	$-$1439.480826	$-$1439.922261	$-$1439.913577	$-$1123.812244	$-$1123.857483
$E^{{pair}}_n$	0	0	$-$1.716956	$-$1.717024	$-$1.703028	$-$1.703045	$-$12.467146	$-$12.466964
$E^{{pair}}_p$	0	0	$-$1.528611	$-$1.528643	$-$1.584308	$-$1.584480	0	0
$E^{{kin}}_n$	2525.991268	2525.991925	1974.613878	1974.613824	1973.986024	1973.980663	1340.457995	1340.668648
$E^{{kin}}_p$	1334.854760	1334.854465	1118.313614	1118.313442	1118.495643	1118.487818	830.735396	830.848077
$E_{{cen}}$	$-$6194.978513	$-$6194.978930	$-$4944.027994	$-$4944.027545	$-$4943.869108	$-$4943.856093	$-$3475.705844	$-$3476.043789
$E_{{SO }}$	$-$96.374920	$-$96.375003	$-$80.186775	$-$80.186826	$-$80.216433	$-$80.214900	$-$49.167364	$-$49.196956
$E_{{dir}}$	827.607126	827.607885	602.810399	602.810352	602.694020	602.697867	366.472441	366.503834
$E_{{exc}}$	$-$31.248467	$-$31.248462	$-$25.935909	$-$25.935905	$-$25.935633	$-$25.935528	$-$19.102496	$-$19.103705
$E_{{stab}}$	8.1E-09	3.5E-11	1.0E-08	3.4E-06	9.6E-09	3.8E-06	9.9E-09	8.8E-08
$E_{{tot}}$	$-$1634.148747	$-$1634.148120	$-$1357.658354	$-$1357.658322	$-$1358.132823	$-$1358.127702	$-$1018.777019	$-$1018.790854

Table 1 displays the results of test calculations performed for the SLy4 Skyrme interaction [6] and for the mixed zero-range pairing force [8]: $V(\vec{r})$= $V_0\left(1 - {\rho(\vec{r})}/{\rho_0} \right)$ for $\rho_0$=0.32fm$^{-3}$. The cutoff energy of $\epsilon_{{cut}}$=60MeV was used for summing up contributions of the HFB quasiparticle states to density matrices [2]. For a given phase space, the strength of the pairing force $V_0$ was adjusted so as to reproduce the experimental neutron pairing gap in $^{120}$Sn. The resulting values are $V_0$=$-$285.88, $-$284.10, and $-$284.36MeVfm$^3$ for the HO (THO) bases of 680 and 3276 states, and for the radial box of $R_{{box}}$=30fm, respectively. The radial HFBRAD calculations were performed with 300 points (i.e., the $\Delta{r}$=0.1fm grid spacing), and the wave functions were included up to $j_{{max}}$=39/2. We checked that even with $j_{{max}}$=33/2, all energies were stable within 1eV. The nucleon-mass and elementary-charge parameters were fixed at $\hbar^2/2m$=20.73553MeVfm$^2$ and $e^2$=1.439978MeVfm, respectively.

Table 1 displays the following quantities: $N_0$ is the maximum number of the HO oscillator quanta included in the basis (for the deformed basis we give the numbers of quanta in the perpendicular ($N_\perp$) and axial ($N_z$) directions); $N_{{st}}$ is the number of the lowest deformed HO states included in the basis; $N^{{qp}}_n$ and $N^{{qp}}_p$ are the numbers of (doubly degenerate) neutron and proton quasiparticle states with equivalent single-particle energies [2] below the cutoff energy $\epsilon_{{cut}}$; $b_\perp$ and $b_z$ are the oscillator constants in the perpendicular and axial directions; $\lambda_n$ and $\lambda_p$ are the neutron and proton Fermi energies, which, for vanishing pairing correlations, are taken as the s.p. energies of the last occupied states; $\Delta_n$ and $\Delta_p$ are the average pairing gaps [7]; $R_n$ and $R_p$ are the rms radii; $Q_n$ and $Q_p$ are the quadrupole moments $\langle 2z^2-x^2-y^2\rangle$; $\epsilon ^{{gs}}_n$ and $\epsilon ^{{gs}}_p$ are the s.p. energies of the most bound neutron and proton states; $\Sigma^{\epsilon }_n$ and $\Sigma^{\epsilon }_p$ are sums of the canonical energies weighted by the corresponding occupation probabilities; $E^{{pair}}_n$ and $E^{{pair}}_p$ are the pairing energies; $E^{{kin}}_n$ and $E^{{kin}}_p$ are the kinetic energies; $E_{{cen}}$ and $E_{{SO }}$ are the energies corresponding to the central and spin-orbit parts of the Skyrme energy density functional; $E_{{dir}}$ and $E_{{exc}}$ are the direct and exchange parts of the Coulomb energy; and $E_{{stab}}$ is the stability energy characterizing the level of self-consistency. In the code HFODD, $E_{{stab}}$ is estimated from the sum of s.p. energies [9]; in the code HFBTHO $E_{{stab}}$ is estimated from the maximum difference of all matrix elements of s.p. potentials calculated in two consecutive iterations; and in the code HFBRAD it is calculated as a variance of the total binding energy, $E_{{tot}}$, over the last five iterations.

Calculations for $^{208}$Pb yield a spherical solution with vanishing pairing gaps. HFBTHO and HFODD give the total binding energies that differ by 627eV, and this difference can be (primarily) traced back to the direct Coulomb energy. We have checked that without the Coulomb interaction, this difference decreases to 202eV. The axial-basis HFBTHO calculation gives a very small total quadrupole moment of 39$\mu$b. This suggests that the THO basis generates a slight deviation from the spherical symmetry due to a different numerical treatment of $z$- and $\perp$-direction. In this respect, HFODD calculations should be considered more accurate.

Calculations for $^{168}$Er performed within a spherical HO basis, $b_\perp$=$b_z$, yield a well-deformed and weakly paired prolate ground state. Here, the total binding energies and quadrupole moments obtained within HFBTHO and HFODD differ only by 32eV and 5$\mu$b, respectively. When the same calculation is performed in a deformed HO basis, $b_\perp$$\neq$$b_z$, the differences grow to 5.1keV and 207$\mu$b, respectively. Again, without the Coulomb interaction, the difference in the total binding energy is only 96eV. It is seen that by employing the deformed basis, the binding energy decreases, as expected.

Comparison with the coordinate-space code HFBRAD for $^{120}$Sn shows that $E_{{tot}}$ in HFBTHO is correct up to 14keV for $N_0$=25. However, the kinetic energy still differs by as much as 221keV, which is compensated by a similar difference in the interaction energy. Within the HO basis and $N_0$=25, the corresponding differences are larger: 41 and 337keV. The analogous differences obtained for $N_0$=20 are 142 and 1103keV (THO), and 152 and 964keV (HO), respectively. Nevertheless, the above comparison shows that the $N_0$=20 calculations yield $E_{{tot}}$ with a precision of a couple of hundred keV.