# Mean-field nuclear structure calculations on a Basis-Spline Galerkin lattice

## 1 Introduction

In recent years, the area of nuclear structure physics has shown substantial progress and rapid growth . With detectors such as GAMMASPHERE and EUROGAM, the limits of total angular momentum and deformation in atomic nuclei have been explored, and new neutron rich nuclei have been identified in spontaneous fission studies. Gamma-ray detectors under development such as GRETA  will have improved resolving power and should allow for the identification of weakly populated states never seen before in nuclei. Particularly exciting is the proposed construction of a next-generation ISOL FACILITY in the United States which has been been identified in the 1996 DOE/NSAC Long Range Plan  as the highest priority for major new construction.

These experimental developments as well as recent advances in computational physics have sparked renewed interest in nuclear structure theory. In contrast to the well-understood behavior near the valley of stability, there are many open questions as we move towards the proton and neutron driplines and towards the limits in mass number (superheavy region). While the proton dripline has been explored experimentally up to Z=83, the neutron dripline represents mostly “terra incognita”. In these exotic regions of the nuclear chart, one expects to see several new phenomena: near the neutron dripline, the neutron-matter distribution will be very diffuse and of large size giving rise to “neutron halos” and “neutrons skins”. One also expects new collective modes associated with this neutron skin, e.g. the “scissors” vibrational mode or the “pygmy” resonance. In proton-rich nuclei, we have recently seen both spherical and deformed proton emitters; this “proton radioactivity” is caused by the tunneling of weakly bound protons through the Coulomb barrier. The investigation of the properties of exotic nuclei is also essential for our understanding of nucleosynthesis in stars and stellar explosions (rp- and r-process). Our primary goal is to carry out high-precision nuclear structure calculations in connection with Radioactive Ion Beam Facilities. Some of the topics of interest are the effective N-N interaction at large isospin, large pairing correlations and their density dependence, neutron halos/skins, and proton radioactivity. Specifically, we are interested in calculating observables such as the total binding energy, charge radii, densities $`\rho _{p,n}(𝐫)`$, separation energies for neutrons and protons, pairing gaps, and potential energy surfaces.

There are many theoretical approaches to nuclear structure physics. For lack of space, we mention only three of these: in the macroscopic \- microscopic method, one combines the liquid drop / droplet model with a microscopic shell correction from a deformed single-particle shell model (Möller and Nix , Nazarewicz et al. ). For relatively light nuclei, it is possible to diagonalize the nuclear Hamiltonian in a shell model basis. Barrett et al.  have recently carried out large-basis no-core shell model calculations for p-shell nuclei. A different approach to the interacting nuclear shell model is the Shell Model Monte Carlo (SMMC) method developed by Dean et al.  which does not involve matrix diagonalization but a path integral over auxiliary fields. This method has been applied to fp-shell and medium-heavy nuclei. Finally, for heavier nuclei one utilizes either nonrelativistic  or relativistic  mean field theories.

## 2 Outline of the theory: HFB formalism in coordinate space

As we move away from the valley of stability, surprisingly little is known about the pairing force: For example, what is its density dependence? Large pairing correlations are expected near the drip lines which are no longer a small residual interaction. Neutron-rich nuclei are expected to be highly superfluid due to continuum excitation of neutron “Cooper pairs”. The Hartree-Fock-Bogoliubov (HFB) theory unifies the HF mean field theory and the BCS pairing theory into a single selfconsistent variational theory. The main challenge in the theory of exotic nuclei near the proton or neutron drip line is that the outermost nucleons are weakly bound (which implies a very large spatial extent), and that the weakly-bound states are strongly coupled to the particle continuum. This represents a major problem for mean field theories that are based on the traditional shell model basis expansion method in which one utilizes bound harmonic oscillator basis wavefunctions. As illustrated in Figure 1 a weakly bound state can still be reasonably well represented in the oscillator basis, but this is no longer true for the continuum states. In fact, Nazarewicz et al.  have shown that near the driplines the harmonic oscillator basis expansion does not converge even if $`N=50`$ oscillator quanta are used. This implies that one either has to use a continuum-shell model basis or one has to solve the problem directly on a coordinate space lattice. We have chosen the latter method.

Several years ago, Umar et al.  have developed a three-dimensional HF code in Cartesian coordinates using the Basis-Spline discretization technique. The program is based on a density dependent effective N-N interaction (Skyrme force) which also includes the spin-orbit interaction. The code has proven efficient and extremely accurate; it incorporates BCS and Lipkin-Nogami pairing, and various constraints. The configuration space Hartree-Fock approach has had great successes in predicting systematic trends in the global properties of nuclei, in particular the mass, radii, and deformations across large regions of the periodic table.

So far, our attempts to generalize this 3D code to include self-consistent pairing forces (Hartree-Fock-Bogoliubov theory on the lattice) have proven too ambitious. The reason may be the lack of a suitable damping operator in 3D. We have therefore taken a different approach and developed a new Hartree-Fock + BCS pairing code in cylindrical coordinates for axially symmetric nuclei, based on the Galerkin method with B-Spline test functions . The new code has been written in Fortran 90 and makes extensive use of new data concepts, dynamic memory allocation and pointer variables. Extending this code, we believe that it will be easier to implement HFB in 2D because one can use highly efficient LAPACK routines to diagonalize the lattice Hamiltonian and does not necessarily rely on a damping operator.

We outline now our basic theoretical approach for lattice HFB. As is customary, we start by expanding the nucleon field operator into a complete orthonormal set of s.p. basis states $`\varphi _i`$

$$\widehat{\psi }^{}(𝐫,s)=\underset{i}{}\widehat{c}_i^{}\varphi _i^{}(𝐫,s)$$
(1)   
which leads to the Hamiltonian in occupation number representation

$$\widehat{H}=\underset{i,j}{}<i|t|j>\widehat{c}_i^{}\widehat{c}_j+\frac{1}{4}\underset{i,j,m,n}{}<ij|\stackrel{~}{v}^{(2)}|mn>\widehat{c}_i^{}\widehat{c}_j^{}\widehat{c}_n\widehat{c}_m.$$
(2)   
Like in the BCS theory, one performs a canonical transformation to quasiparticle operators $`\widehat{\beta },\widehat{\beta }^{}`$

$$\left(\begin{array}{c}\widehat{\beta }\\ \widehat{\beta }^{}\end{array}\right)=\left(\begin{array}{cc}U^{}& V^{}\\ V^T& U^T\end{array}\right)\left(\begin{array}{c}\widehat{c}\\ \widehat{c}^{}\end{array}\right).$$
(3)   
The HFB ground state is defined as the quasiparticle vacuum

$$\widehat{\beta }_k|\mathrm{\Phi }_0>=0.$$
(4)   
The basic building blocks of the theory are the normal density

$$\rho _{ij}=<\mathrm{\Phi }_0|\widehat{c}_j^{}\widehat{c}_i|\mathrm{\Phi }_0>=(V^{}V^T)_{ij}$$
(5)   
and the pairing tensor

$$\kappa _{ij}=<\mathrm{\Phi }_0|\widehat{c}_j\widehat{c}_i|\mathrm{\Phi }_0>=(V^{}U^T)_{ij}$$
(6)   
from which one can form the generalized density matrix

$$=\left(\begin{array}{cc}\rho & \kappa \\ \kappa ^{}& 1\rho ^{}\end{array}\right)^{}=,^2=.$$
(7)   
Using the definition of the HFB ground state energy

$$E^{}()=<\mathrm{\Phi }_0|\widehat{H}\lambda \widehat{N}|\mathrm{\Phi }_0>$$
(8)   
we derive the equations of motion from the variational principle

$$\delta [E^{}()\mathrm{tr}\mathrm{\Lambda }(^2)]=0$$
(9)   
which results in the standard HFB equations

$$[,]=0$$
(10)   
with the generalized single-particle Hamiltonian

$$=\left(\begin{array}{cc}h& \mathrm{\Delta }\\ \mathrm{\Delta }^{}& h^{}\end{array}\right);h=E^{}/\rho ,\mathrm{\Delta }=E^{}/\kappa ^{}.$$
(11)   
Our goal is to transform to a coordinate space representation and solve the resulting differential equations on a lattice. For this purpose, we define two types of quasi-particle wavefunctions $`\varphi _1,\varphi _2`$

$$\varphi _1^{}(E_n,𝐫\sigma )=\underset{i}{}U_{in}\varphi _i(𝐫\sigma ),\varphi _2(E_n,𝐫\sigma )=\underset{i}{}V_{in}^{}\varphi _i(𝐫\sigma )$$
(12)   
in terms of which the particle density matrix for the HFB ground state assumes a very simple mathematical structure 

$`\rho _0(𝐫,\sigma ,𝐫^{},\sigma ^{})=<\mathrm{\Phi }_0|\widehat{\psi }^{}(𝐫^{}\sigma ^{})\widehat{\psi }(𝐫\sigma )|\mathrm{\Phi }_0>`$   
$`={\displaystyle \underset{i,j}{}}\rho _{ij}\varphi _i(𝐫\sigma )\varphi _j^{}(𝐫^{}\sigma ^{})={\displaystyle \underset{E_n>0}{\overset{\mathrm{}}{}}}\varphi _2(E_n,𝐫\sigma )\varphi _2^{}(E_n,𝐫^{}\sigma ^{}).`$ (13)   
In a similar fashion we obtain for the pairing tensor

$`\kappa _0(𝐫,\sigma ,𝐫^{},\sigma ^{})=<\mathrm{\Phi }_0|\widehat{\psi }(𝐫^{}\sigma ^{})\widehat{\psi }(𝐫\sigma )|\mathrm{\Phi }_0>`$   
$`={\displaystyle \underset{i,j}{}}\kappa _{ij}\varphi _i(𝐫\sigma )\varphi _j(𝐫^{}\sigma ^{})={\displaystyle \underset{E_n>0}{\overset{\mathrm{}}{}}}\varphi _2(E_n,𝐫\sigma )\varphi _1^{}(E_n,𝐫^{}\sigma ^{}).`$ (14)   
For certain types of effective interactions (e.g. Skyrme forces) the HFB equations in coordinate space are local and have a structure which is reminiscent of the Dirac equation 

$$\left(\begin{array}{cc}(h\lambda )& \stackrel{~}{h}\\ \stackrel{~}{h}& (h\lambda )\end{array}\right)\left(\begin{array}{c}\varphi _1(𝐫)\\ \varphi _2(𝐫)\end{array}\right)=E\left(\begin{array}{c}\varphi _1(𝐫)\\ \varphi _2(𝐫)\end{array}\right),$$
(15)   
where $`h`$ is the “particle” Hamiltonian and $`\stackrel{~}{h}`$ denotes the “pairing” Hamiltonian.

The various terms in $`h`$ depend on the particle densities $`\rho _q(𝐫)`$ for protons and neutrons, on the kinetic energy density $`\tau _q(𝐫)`$, and on the spin-current tensor $`J_{ij}(𝐫)`$. The pairing Hamiltonian $`\stackrel{~}{h}`$ has a similar structure, but depends on the pairing densities $`\stackrel{~}{\rho }_q(𝐫),\stackrel{~}{\tau }_𝐪(𝐫)`$ and $`\stackrel{~}{J}_{ij}(𝐫)`$ instead. Because of the structural similarity between the Dirac equation and the HFB equation in coordinate space, we encounter here similar computational challenges: for example, the spectrum of quasiparticle energies $`E`$ is unbounded from above and below. The spectrum is discrete for $`|E|<\lambda `$ and continuous for $`|E|>\lambda `$. In the case of axially symmetric nuclei, the spinor wavefunctions $`\varphi _1(𝐫)`$ and $`\varphi _2(𝐫)`$ have the structure

$$\psi ^\mathrm{\Omega }(\varphi ,r,z)=\frac{1}{\sqrt{2\pi }}\left(\begin{array}{c}e^{i(\mathrm{\Omega }\frac{1}{2})\varphi }U(r,z)\\ e^{i(\mathrm{\Omega }+\frac{1}{2})\varphi }L(r,z)\end{array}\right).$$
(16)   
## 3 Computational method: Spline-Galerkin lattice representation

For nuclei near the p/n driplines, we overcome the convergence problems of the traditional shell-model expansion method by representing the nuclear Hamiltonian on a lattice utilizing a Basis-Spline expansion . B-Splines $`B_i^M(x)`$ are piecewise-continuous polynomial functions of order $`(M1)`$. They represent generalizations of finite elements which are B-splines with $`M=2`$. A set of fifth-order B-Splines is shown in Figure 2.

Let us now discuss the Galerkin method with B-Spline test functions. We consider an arbitrary (differential) operator equation

$$𝒪\overline{f}(x)\overline{g}(x)=0.$$
(17)   
Special cases include eigenvalue equations of the HF/HFB type where $`𝒪=h`$ and $`\overline{g}(x)=E\overline{f}(x)`$. We assume that both $`\overline{f}(x)`$ and $`\overline{g}(x)`$ are well approximated by Spline functions

$$\overline{f}(x)f(x)\underset{i=1}{\overset{𝒩}{}}B_i^M(x)a^i,\overline{g}(x)g(x)\underset{i=1}{\overset{𝒩}{}}B_i^M(x)b^i.$$
(18)   
Because the functions $`f(x)`$ and $`g(x)`$ are approximations to the exact functions $`\overline{f}(x)`$ and $`\overline{g}(x)`$, the operator equation will in general only be approximately fulfilled

$$𝒪f(x)g(x)=R(x).$$
(19)   
The quantity $`R(x)`$ is called the residual; it is a measure of the accuracy of the lattice representation. We multiply the last equation from the left with the spline function $`B_k(x)`$ and integrate over $`x`$

$$v(x)𝑑xB_k(x)𝒪f(x)v(x)𝑑xB_k(x)g(x)=v(x)𝑑xB_k(x)R(x).$$
(20)   
We have included a volume element weight function $`v(x)`$ in the integrals to emphasize that the formalism applies to arbitrary curvilinear coordinates. Various schemes exist to minimize the residual function $`R(x)`$; in the Galerkin method one requires that there be no overlap between the residual and an arbitrary B-spline function

$$v(x)𝑑xB_k(x)R(x)=0.$$
(21)   
This so called Galerkin condition amounts to a global reduction of the residual. Applying the Galerkin condition and inserting the B-Spline expansions for $`f(x)`$ and $`g(x)`$ results in

$$\underset{i}{}\left[v(x)𝑑xB_k(x)𝒪B_i(x)\right]a^i\underset{i}{}\left[v(x)𝑑xB_k(x)B_i(x)\right]b^i=0.$$
(22)   
Defining the matrix elements

$$𝒪_{ki}=v(x)𝑑xB_k(x)𝒪B_i(x),G_{ki}=v(x)𝑑xB_k(x)B_i(x)$$
(23)   
transforms the (differential) operator equation into a matrix $`\times `$ vector equation

$$\underset{i}{}𝒪_{ki}a^i=\underset{i}{}G_{ki}b^i$$
(24)   
which can be implemented on modern vector or parallel computers with high efficiency. The matrix $`G_{ki}`$ is sometimes referred to as the Gram matrix; it represents the nonvanishing overlap integrals between different B-Spline functions (see Fig. 2). We eliminate the expansion coefficients $`a^i,b^i`$ in the last equation by introducing the function values at the lattice support points $`x_\alpha `$ including both interior and boundary points.

The upper $`(U)`$ and lower $`(L)`$ components of the spinor wavefunctions defined earlier are represented on the 2-D lattice $`(r_\alpha ,z_\beta )`$ by a product of Basis Splines $`B_i(x)`$ evaluated at the lattice support points

$$U(r_\alpha ,z_\beta )=\underset{i,j}{}B_i(r_\alpha )B_j(z_\beta )U^{ij},L(r_\alpha ,z_\beta )=\underset{i,j}{}B_i(r_\alpha )B_j(z_\beta )L^{ij}.$$
(25)   
We are also extending our previous B-spline work to include nonlinear grids. Use of a nonlinear lattice should be most useful for loosely bound systems near the proton or neutron drip lines. Non-Cartesian coordinates necessitate the use of fixed endpoint boundary conditions; much effort has been directed toward improving the treatment of these boundaries .

## 4 Numerical tests and results

We expect our Spline techniques to be superior to the traditional harmonic oscillator basis expansion method in cases of very strong nuclear deformation. To illustrate this point, we have performed a numerical test using a phenomenological (Woods-Saxon) deformed shell model potential. We calculate the single-particle energy spectrum for neutrons in $`{}_{}{}^{40}Ca`$ for quadrupole deformations ranging from strong oblate ($`\beta _2=1.25`$) to extreme prolate ($`\beta _2=+2.25`$). The results are shown in Fig. 3. Apparently, for $`\beta _2=0`$ we correctly reproduce the spherical shell structure of magic nuclei. As $`\beta _2`$ approaches large positive values our s.p. potential approaches the structure of two separated potential wells; as expected, we observe pairs of levels with opposite parity that are becoming degenerate in energy. The largest quadrupole deformation corresponds physically to a symmetric fission configuration. Clearly, such configurations cannot be described in a single oscillator basis, which confirms the numerical superiority of the B-Spline lattice technique.

In a second test calculation, we have investigated the properties of a nucleus near the neutron drip line. During the last decade the discovery of a ‘neutron halo’ in several neutron-rich isotopes generated a great deal of interest in the area of weakly bound quantum systems. The halo effect was first observed in $`{}_{3}{}^{}{}_{}{}^{11}`$Li, which consists of three protons and six neutrons in a central core and two planetary neutrons which comprise the halo. By adjusting the depth of the Woods-Saxon potential so that the separation energy of the last bound neutron is only $`10`$ keV, i.e. very close to the continuum, we were able to determine this neutron wavefunction on the lattice which shows a very large spatial extent (see Fig. 4). We conclude that the B-Spline lattice techniques are well-suited for representing weakly bound states near the drip lines; a similar calculation in the basis expansion method would require a large number of oscillator shells.

We now discuss our numerical results for the selfconsistent Hartree-Fock calculations with Skyrme-M interaction and BCS pairing. This is a special case of the HFB equation with a constant pairing matrix element. In Fig. 5 we display the proton density for a heavy nucleus, $`{}_{64}{}^{}{}_{}{}^{154}`$Gd, calculated with our new 2-D (HF+BCS) code. It should be noted that ALL 154 nucleons are treated dynamically (no inert core approximation). The theoretical charge density looks quite similar to the experimental result which is shown on the right hand side.

For several spherical nuclei, we have also compared the selfconsistent s.p. energy levels of our 2-D Spline-Galerkin code with a fully converged 1-D radial calculation. The result is shown in Table 1.

### 4.1 Plans and Future Directions

Having validated our new (HF+BCS) code on a 2D lattice with the Spline-Galerkin method, we plan to proceed as follows: We are currently working on the 2D HFB implementation with a pairing delta-force. After that, we will generalize the code utilizing the full SkP force with mean pairing field and pairing spin-orbit term. We will also add appropriate constraints, e.g. $`Q_{20},Q_{30},\omega j_x`$ for calculating potential energy surfaces and rotational bands. As we compare the observables (e.g. total binding energy, charge radii, densities $`\rho _{p,n}(𝐫)`$, separation energies for neutrons and protons, pairing gaps) with experimental data from the RIB facilities, it will almost certainly be necessary to develop new effective N-N interactions as we move farther away from the stability line towards the p/n drip lines.

## Acknowledgments

This research project was sponsored by the U.S. Department of Energy under contract No. DE-FG02-96ER40975 with Vanderbilt University. Some of the numerical calculations were carried out on CRAY supercomputers at the National Energy Research Scientific Computing Center (NERSC) at Lawrence Berkeley National Laboratory. We also acknowledge travel support from the NATO Collaborative Research Grants Program.