Adaptive Multi-resolution 3D Hartree-Fock-Bogoliubov Solver for Nuclear Structure
Junchen Pei, George Fann, Robert Harrison, Witold Nazarewicz, Yue Shi, Scott Thornton
AAdaptive Multi-resolution 3D Hartree-Fock-Bogoliubov Solver for Nuclear Structure
J.C. Pei ( 裴 俊 琛 ),
1, 2, 3
G.I. Fann, R. J. Harrison,
5, 6
W. Nazarewicz,
2, 7, 8
Yue Shi ( 石 跃 ),
2, 3 and S. Thornton State Key Laboratory of Nuclear Physics and Technology,School of Physics, Peking University, Beijing 100871, China Department of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USA Joint Institute for Nuclear Physics and Applications,Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA Institute for Advanced Computational Science, Stony Brook University, Stony Brook, New York 11794, USA Computational Science Center, Brookhaven National Laboratory, Upton, New York 11973, USA Physics Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA Institute of Theoretical Physics, Faculty of Physics,University of Warsaw, ul. Ho˙za 69, PL-00681 Warsaw, Poland
Background:
Complex many-body systems, such as triaxial and reflection-asymmetric nuclei, weakly-boundhalo states, cluster configurations, nuclear fragments produced in heavy-ion fusion reactions, cold Fermi gases,and pasta phases in neutron star crust, they are all characterized by large sizes and complex topologies, in whichmany geometrical symmetries characteristic of ground-state configurations are broken. A tool of choice to studysuch complex forms of matter is an adaptive multi-resolution wavelet analysis. This method has generated muchexcitement since it provides a common framework linking many diversified methodologies across different fields,including signal processing, data compression, harmonic analysis and operator theory, fractals, and quantum fieldtheory.
Purpose:
To describe complex superfluid many-fermion systems, we introduce an adaptive pseudo-spectralmethod for solving self-consistent equations of nuclear density functional theory in three dimensions, withoutsymmetry restrictions.
Methods:
The numerical method is based on the multi-resolution and computational harmonic analysis tech-niques with multiwavelet basis. The application of state-of-the-art in parallel programming techniques includesophisticated object oriented templates which parses the high-level code into distributed parallel tasks with amultithread task queue scheduler for each multicore node. The inter-node communications are asynchronous.The algorithm is variational and is capable of solving coupled complex-geometric systems of equations adaptively,with functional and boundary constraints, in a finite spatial domain of very large sizes, limited by existing parallelcomputer memory. For smooth functions, user defined finite precision is guaranteed.
Results:
The new adaptive multi-resolution Hartree-Fock-Bogoliubov (HFB) solver madness-hfb is bench-marked against a two-dimensional coordinate-space solver hfb-ax based on B-spline technique and three-dimensional solver hfodd based on the harmonic oscillator basis expansion. Several examples are considered,including self-consistent HFB problem for spin-polarized trapped cold fermions and Skyrme-Hartree-Fock (+BCS)problem for triaxial deformed nuclei.
Conclusions:
The new madness-hfb framework has many attractive features when applied to nuclear andatomic problems involving many-particle superfluid systems. Of particular interest are weakly-bound nuclearconfigurations close to particle drip lines, strongly elongated and dinuclear configurations such as those presentin fission and heavy ion fusion, and exotic pasta phases that appear in the neutron star crust.
PACS numbers: 21.60.Jz,31.15.E-,03.65.Ge,07.05.Tp,67.85.-d,03.75.Hh
I. INTRODUCTION
The roadmap for nuclear structure theory includesQCD-derived (or inspired) nuclear interactions, ab-initiocalculations for light and medium nuclei, configurationinteraction approaches for near-magic systems, and den-sity functional theory and its extensions for heavy, com-plex nuclei [1]. On the road to the quantitative un-derstanding of nuclear structure and reactions, high-performance computing plays an increasingly importantrole. As stated in the recent decadal survey of nuclearphysics [2] “High performance computing provides an-swers to questions that neither experiment nor analytictheory can address; hence, it becomes a third leg support- ing the field of nuclear physics.” Largest collaborationsin computational nuclear structure and reactions involvenuclear theorists, computer scientists, and applied math-ematicians to break analytic, algorithmic, and computa-tional barriers [1, 3]. This paper offers an example of sucha joint collaborative effort in the area of nuclear DensityFunctional Theory (DFT).A key element of any DFT framework is a HFBsolver that computes self-consistent solutions of HFB(or Bogoliubov-de Gennes) equations. Traditionally, theHFB solvers in nuclear physics are based on the ba-sis expansion method, usually employing harmonic os-cillator wave functions [4–7]. These methods are veryefficient but they require huge bases for cases involv-ing weakly-bound systems and large deformations [8, a r X i v : . [ nu c l - t h ] J u l madness-hfb solver for HFB equations and Hartree-Fock (HF) equa-tions, which is a multi-resolution, adaptive spectral ap-proximations based solver, using a multiwavelet basis,with a scalable parallel implementation [18]. The newframework is applied to polarized ultracold Fermi gasesin elongated optical traps as well as triaxial nuclei. Inboth cases, we will demonstrate the capability of verylarge box calculations which is essential for descriptionsof complex geometries and topologies.This paper is organized as follows. Section II briefly in-troduce the multiresolution mathematics, low-separationrank approximation, and parallel runtime environment.The iterative algorithm applied in madness-hfb is pre-sented in Sec. III. In Sec. IV, we benchmark madness-hfb solutions for cold fermions and nuclei. Finally, con-clusions are given in Sec. V. II. MADNESS-HFB FRAMEWORK
Our implementation of madness-hfb uses the Mul-tiresolution Adaptive Numerical Environment for Scien-tific Simulations (MADNESS) framework [18]. MAD-NESS is based on computational harmonic analysis andnonlinear approximations using Alpert’s multiwaveletbasis [17, 19, 20] to represent functions. Fast parallelcode development and scalable performance have beenpossible due to the ease of programming based on object-oriented abstractions for interprocessor communications,multithreading and mathematical operations.
A. Mathematics of MADNESS
The mathematics implemented in the MADNESS soft-ware are based on multiresolution analysis (MRA) [19,20], nonlinear approximations, and pseudo-spectral tech-niques. There are two types of techniques used in MAD-NESS to approximate functions and operators. The firstis the use of multiresolution analysis based on Alpert’smultiwavelets [17]. The second technique is the use of thelow-separation rank approximations of Green’s functionsbased on Gaussian functions [21, 22]. In the following,we follow the notation and derivations of Ref. [20].
1. Multiresolution analysis with wavelets
The application of MRA separates the behavior offunctions and operators at different length scales in asystematic expansion. A consequence of the separationof scales is that each operator and wave function hasa naturally independent adaptive refinement structure,reflected in terms of significant expansion coefficients ofdesired precision. The thresholding and truncation of ex-pansion coefficients below a user-defined error providesadaptive blocks of non-trivial coefficients for a pseudo-spectral expansion. The union of the domains of themultiwavelets with non-zero coefficients provide an adap-tive dyadic spatial localization of the relevent contribu-tions for the corresponding refinement levels. In 1D, thenon-zero sets define an adaptive dyadic refinement andcorrespondingly in 3D a pruned octtree type refinement.The MRA representation used in MADNESS is anal-ogous to that used in an adaptive hp-SEM (spectral el-ement method), which employs elements of variable size h and piecewise-polynomial approximations of degree p .By suitably refining the mesh through h -refinements (di-viding the volume elements into smaller pieces) and p -refinements (increasing the polynomial degree in the ex-pansion within the elements) one can reach exponentialconvergence [23]. In MADNESS, for each function oroperator, the union of the domains of the multiwaveletbasis functions with non-zero coefficients, after threshold-ing, defines an adaptive and heirarchical h -structure andthe associated multiwavelets form the set of the piecewisepolynomials up to order p . Thus, there are multiple h − p refinement structures that are used simultaneously.The basis of scaling functions in 1D is constructed interms of the normalized Legendre polynomials rescaledto the unit interval (0 ,
1) and zero elsewhere. For eachlevel n (defining the volume refinement), the rescaled andtranslated basis function is given by: φ nil ( x ) = 2 n/ φ i (2 n x − l ) , (1)where φ i ( x ) = √ i + 1 P i (2 x − P i ( x ) being theLegendre polynomial on ( − , l = 0 , ..., n −
1. The basis functions (1) at level n havedomain of width 2 − n .Let V n = { φ nil ( x ) , i = 0 , ..., k − } be the span ofthe subspace at level n . Let W = { ψ i ( x ) } denote anorthonormal basis which spans the difference subspace V − V . These functions are called multiwavelets. Aswith the scaling functions let ψ nil ( x ) and W n denote therescaled and shifted multiwavelets and the correspond-ing subspace spanned by these functions at level n . Thedefinition of scaling functions and multiwavelets definesan ascending sequence of subspaces V ⊂ V ⊂ V ... ⊂ V n (2)and V n = V ⊕ W ⊕ ... ⊕ W n − , (3)where the ⊕ denotes orthogonal sum. The dimensionof V i is greater than dimension of the subspace V i − ;thus, the basis functions of V i − and W i − can be writ-ten exactly in terms of the basis functions of V i . Theseheirarchical linear algebraic relations between the basesdefines the 2-scale refinement structure between the co-efficients at level i − i , and fundamentally definesthe adaptive structure with a given threshold truncation.A smooth function f ( x ) in the subspace V n can beapproximated in terms of scaling functions as: f ( x ) = n − (cid:88) l =0 k − (cid:88) j =0 s njl φ njl ( x ) . (4)Represented in the multiwavelet basis, f ( x ) is f ( x ) = k − (cid:88) j =0 s j φ j ( x ) + k − (cid:88) j =0 n − (cid:88) m =1 2 m − (cid:88) l =0 d mjl ψ mjl ( x ) , (5)with s njl = (cid:82) − n ( l +1)2 − n l f ( x ) φ njl dx and d mjl = (cid:82) − m ( l +1)2 − m l f ( x ) ψ mjl dx .In the discussion above, we described the representa-tions based on multiwavelets in 1D. In 3D applications,we use tensor products of 1D multiwavelets as well asscaling functions in non-standard form. Figure 1 illus-trates the multiresolution structure of sample wave func-tions. For smooth functions the computational methodolo-gies are guaranteed to approximate the solutions to thedesired user precision (cid:15) , with respect to the relative norm,with the correct number of digits specified by the error.The estimate is based on truncating the difference coef-ficients in the multiwavelet expansion, || d nl || = (cid:115)(cid:88) j | d njl | ≤ (cid:15) min(1 , − n L ) , (6)where L is the minimum of the width of the computa-tional domain. FIG. 1. Pedagogical illustration of adaptive representations inMADNESS-HFB. Top: (a) the modulus squared of the single-neutron wave function corresponding to the single-particle en-ergy of − .
214 MeV obtained in madness-hf calculations for
Mo (see Sec. IV B for details), and (b) the correspond-ing spectral refinement structure. Bottom: (c) the modulussquared of the single-proton wave function corresponding tothe energy eigenvalue − .
272 MeV in
Mo and its adaptivespectral structure (d). Notice that the refinement structurefor the proton wave function is similar to a truncated octtreetype of refinement but the structure for neutron wave functionis more complicated, especially at the finer level, see insets inpanels (b) and (d).
2. Multiresolution
For the one-body Schr¨odinger equation,( − ∆ + V ) ψ = Eψ, (7)the usual diagonalization approach is also derived andused. In this case, given a basis ψ i , a Hamiltonian matrixis formed H i,j = (cid:104) ψ i | − ∆ + V | ψ j (cid:105) ; S i,j = (cid:104) ψ i | ψ j (cid:105) , (8)to form a generalized eigenproblem Hψ = Sψ.
A generalized eigensolver computes the eigenvaluesand the eigenvectors. The eigenvectors are coefficientswith respect to the multiwavelets basis, and they areconverted back to the spectral representation for furthercomputation. The Laplacian ∆, the potential V , andthe wave-functions ψ i are all in MRA form. The deriva-tives of multiwavelets are expanded in terms of multi-wavelets, and the coefficients are tabulated. By linearity,the derivatives of a function can be computed by matrix-vector products, or tensor-tensor products in higher di-mensions, using only the multiwavelet coefficients.This procedure permits computation of “self-consistent” solutions of DFT equations.
3. Low-separation rank approximation of Green’s functions
Recall that the (one-body) Schr¨odinger equation (7)can be rewritten as a Lippmann-Schwinger equation as(∆ + E ) ψ = V ψ. (9)There are several advantages of using the integral form(9) over the differential form (7). Namely, the integralform provides higher accuracy as high-frequency noiseis attenuated not amplified, builds correct asymptotics,good condition number, and is potentially more compu-tationally efficient. In most bases Green’s function rep-resentation is often dense, and the use of multiresolutionanalysis and multiwavelets provides fast algorithms withsparse structure in finite floating arithmetic with guar-anteed precision. If no controlled truncations of the mul-tiwavelet coefficients are performed the representation ofthe Green’s function and its application will be dense.The formal solution of (9) can be written as: ψ ( r ) = (cid:90) ∞−∞ G ( r − r (cid:48) ) V ψ ( r (cid:48) ) dr (cid:48) = ( G (cid:63) V ψ ) , (10)where the Green’s function G ( r ) is the Helmholtz kerneland the symbol (cid:63) represents convolution. If the eigen-value is bound ( E < − kr ) /r where k = √− E . In gen-eral, one works with G = (∆ + E + iε ) − with ε → +0and specifies how to integrate around the poles.For bound-states, a low-separation rank (LSR) expan-sion [21, 22] of the Yukawa potential is used, e − kr r = (cid:88) l σ l e − τ l r + O ( (cid:15) ) . (11)The LSR approximation represents Green’s function interms of a Gaussian expansion. Such a form reducesthe application of 3D convolutions to an set of uncou-pled 1D convolutions with the number of terms scalinglogarithmically with respect to the relative precision (cid:15) .Since the convolution operator is linear, tables of precom-puted transformation matrices with respect to the mul-tiwavelets enable fast applications of convolutions [24]. The technique described above to solve the Schr¨odingerequation can be directly applied to a HF problem, and –after a minor generalization – to HFB equations. B. MADNESS parallel runtime environment
A novel parallel execution runtime environment hasbeen implemented in the MADNESS software library.MADNESS uses one Message Passing Interface processto communicate between nodes, and POSIX Threadswithin a node to exploit shared memory parallelism witha global addressable view of memory space in software.The MADNESS runtime is based on a parallel task-basedcomputing model with a graph-based scheduler and atask queue on each node, to enable distributed multi-threaded computation. A microparser is used to decoupletasks as much as possible but also to detect data depen-dencies so as many independent and out-of-order taskscan execute simultaneously, ensuring correct and mini-mal number of synchronization and thread termination.Although the dedicated use of a core for inter-nodecommunication and a core for handling thread schedul-ing may be a big sacrifice of computational resources, forsupercomputers with large numbers of cores per node,we are able to obtain more than 50% of peak core per-formance for the remaining cores. Most scientific andengineering codes obtain only about 10% of the peak pro-cessor performance.The flexibility of madness-hfb in its design and pro-gramming style permits the solution of multiphysicsproblems with complex geometric structures and bound-ary conditions in large volumes in the coordinate-spaceformulation – limited only by the size of aggregate com-puter memory. Nuclear fission, exotic topologies insuper- and hyperheavy nuclei, neutron star crusts, andcold atoms in elongated traps are some examples whichcan take advantage of these features.
III. MADNESS-HFB STRATEGY
The general HFB equation for a two-component (e.g.,spin-up ↑ and spin-down ↓ ) system of fermions can bewritten as [25–28]: (cid:20) h ↑ − λ ↑ ∆∆ ∗ − h ↓ + λ ↓ (cid:21) (cid:20) u i v i (cid:21) = E i (cid:20) u i v i (cid:21) , (12)where h ↑ and h ↓ are the Hartree-Fock Hamiltoniansfor the spin-up and spin-down components, respectively.The corresponding chemical potentials are denoted as λ ↑ and λ ↓ , and the pairing potential is ∆.There are two standard approaches to solve the HFBequation (12). In the basis expansion method, eigenvec-tors ( u i , v i ) are expressed in terms of a single-particlebasis and the self-consistent procedure applies the HFBHamiltonian matrix diagonalization. The HFB solvers hfbtho [4] (using the cylindrical transformed deformedharmonic-oscillator basis) and hfodd [7] (using theCartesian deformed harmonic-oscillator basis), employedin this work to benchmark madness-hfb belong to thisclass. A second way is to solve the HFB equations inthe coordinate-space by finite difference or finite ele-ment methods [12, 13, 29] or in the momentum spaceusing fast Fourier transform [30]. The strategy appliedin madness-hfb , described in Sec. II, combines featuresfrom these two approaches. The original method hasbeen developed in the context of HF and DFT problemsin computational chemistry [22, 31].To illustrate the self-consistent procedure, let us con-sider a case of an unpolarized system ( h ↑ = h ↓ ) with con-stant effective mass 1 /α . The mean-field Hamiltonianis h ( r ) = − α ∇ U ( r ) , (13)where U ( r ) is the HF potential. As discussed inSec. II A 3, it is convenient to rewrite the HFB equa-tion (12) in a Lippmann-Schwinger form. To this end,in each step of HFB iterations, we introduce the Green’sfunctions G + and G − : G n ± = 1 ± α ∇ + ( E s ± E ni ) , (14)where E ni is the i -th HFB eigenvalue in the n -th iterationstep, and E s is the energy displacement that shifts thepositive-energy HFB eigenvalues so that the associatedGreen’s function is properly defined.To solve the self-consistent HFB eigenproblem, theHFB wave functions can be updated by as follows, u n +1 i = ( G + (cid:63) [( U − λ ) u ni + ∆ v ni + E s u ni ]) (15a) v n +1 i = ( G − (cid:63) [( U − λ ) v ni − ∆ u ni + E s v ni ]) (15b)Following this strategy, in the following section, we use madness-hfb to solve HFB problems with advanced lo-cal energy density functionals for cold fermions and nu-clei. IV. BENCHMARK PROBLEMS
In this section, the madness-hfb framework is bench-marked by solving HFB equations for the trapped uni-tary Fermi gas and HF-BCS equations for a triaxial nu-cleus. The madness-hfb solutions for atoms and nucleiare compared with results of 2D hfb-ax and 3D hfodd calculations, respectively.
A. HFB solver for unitary Fermi gas
The unitary limit of Fermi gas, is characterized by aninfinite s -wave scattering length. Of particular interest are superfluid phases in spin-imbalanced systems, such asthe Fulde-Ferrell-Larkin-Ovchinnikov [32, 33] phase thatexhibits oscillated pairing gaps. The ultracold Fermionsat the unitary limit can be described by the superfluiddensity functional SLDA [34] and its asymmetric exten-sion ASLDA for spin-polarized systems [25].The single-particle Hamiltonian of the ASLDA forasymmetric systems can be written as [25]: h σ = − (cid:126) m ∇ · ( ∇ α σ ( r )) + U σ ( r ) + V ext ( r ) , (16)where σ = ( ↑ , ↓ ) denotes the spin up and spin down com-ponents. The local polarization is denoted as x ( r ) = ρ ↓ ( r ) /ρ ↑ ( r ) with x ( r ) (cid:54)
1, where ρ ↑ ( r ) and ρ ↓ ( r ) aredensities of spin-up and spin-down atoms, respectively.The total polarization of the system is P = ( N ↑ − N ↓ ) /N .The quantity α σ ( x ( r )) is the local effective mass. TheSLDA formalism can be obtained from ASLDA by as-suming x ( r ) = 1, resulting in identical effective massesand Hartree potentials for spin-up and spin-down species.The cold atoms are trapped in an external potential V ext ( x, y, z ) = V (cid:20) − e − ω x y z /η V (cid:21) , (17)where the trap aspect ratio η denotes the elongation ofthe optical trap potential. The equations are normalizedso that (cid:126) = m = ω = 1 (trap units). All other detailspertaining to our SLDA and ASLDA calculations closelyfollow Ref. [28].We first consider an SLDA case of ten particles in aspherical trap with V = 10 and the quasiparticle en-ergy cutoff E cut = 6. The calculations were performedin a 3D box ( − L, L ) with L = 60 With this box andcutoff, the self-consistent HFB solution involves 296 one-quasiparticle eigenfunctions. In the present SLDA andASLDA benchmark calculations, we adopt wavelet orderof p = 8 with a requested truncation precision of (cid:15) = 10 − (see Eq.6).The madness-hfb results were benchmarked using the2D HFB solver hfb-ax . In the hfb-ax calculation,the maximum mesh size is 0 .
3, the order of B-splinesis k = 12, and the box size is R max = Z max = 14. Theeigenvalues and occupation numbers of some of the low-est and highest states from the two codes are comparedin Table I. The agreement is excellent, also for the totalenergy and chemical potential.Next we consider the functional ASLDA, which wasdeveloped to describe polarized Fermi systems. Becauseof non-zero spin polarization, the corresponding HFB so-lutions break time-reversal symmetry. In the first test,we performed madness-hfb and hfb-ax simulations for10 particles with a total polarization of P = 0 . madness-hfb and hfb-ax . Some ofthe eigenvalues are compared in Table II. Note that thecalculation conditions adopted in Table I and Table II TABLE I. Benchmark comparison of madness-hfb and hfb-ax results for 10 particles in the spherical trap without polar-ization. Displayed are: one-quasiparticle energies E i , occu-pations v i , chemical potential λ , and total energy E t . Eachone-quasiparticle state is labelled by means of orbital quan-tum number (cid:96) and parity π = ( − (cid:96) . Note that every solutionis 2 (cid:96) +1-folded degenerate with respect to the magnetic quan-tum number. The numbers in parentheses denote powers of10. The energy is expressed in trap units ( (cid:126) = m = ω =1). madness-hfb hfb-ax i (cid:96) E i v i E i v i (-2)8 0 2.82496 0.60883 2.8250 (-3)10 1 3.44774 2.3162(-2) 3.4477 (-2)21 7 5.54071 3.1957(-5) 5.5407 (-5)22 2 5.58728 3.6548(-3) 5.58728 3.655 (-3)23 4 5.75254 1.8024(-3) 5.7525 E t = 18 . E t = 18 . λ = 2 . λ = 2 . z-axis (trap units) den s i t y ( t r ap un i t s ) ASLDA N =10, P =0.1 FIG. 2. (Color online) Comparison between density distribu-tions ρ ↓ ( r ) and ρ ↑ ( r ) obtained in ASLDA with madness-hfb and hfb-ax for 10 particles in a spherical trap with polariza-tion P = 0 . are the same. It can been seen that the agreement isgood up to the 4th digit since the calculations of localpolarization x ( r ) = ρ ↓ ( r ) /ρ ↑ ( r ) may lose accuracy inboth approaches when both the spin-up and spin-downdensities are very small. In this case, required preci-sion ALSDA should be significantly greater than thatrequested in SLDA calculations.To demonstrate the capability of madness-hfb for ac-curate simulation of large systems, we carried out SLDA TABLE II. Similar as in Table I but for a polarized systemin ASLDA. madness-hfb hfb-ax i E i v i E i v i − . − . (-1)5 1.0157 0.2749 1.016 (-1)8 1.8346 0.6160 1.834
23 4.6417 0.0155(-1) 4.641 (-1)24 4.8158 0.1689(-5) 4.815 (-5) E t = 19 . E t = 19 . ( λ ↑ + λ ↓ ) / . λ ↑ + λ ↓ ) / . N ↑ − N ↓ = 1 . N ↑ − N ↓ = 1 . z-axis (trap units) den s i t y ( t r ap un i t s ) SLDA N =100, P =0.2, η =16 FIG. 3. (Color online) Comparison between density distribu-tions ρ ↑ + ρ ↓ and ρ ↑ − ρ ↓ obtained in SLDA with madness-hfb and hfb-ax for 100 particles with P = 0 . η = 16. simulations for 100 particles with polarization P = 0 . η = 16. The choice of SLDAwas motivated by the above-mentioned loss accuracy ofASLDA caused by a numerical error on x ( r ) at low den-sities (large distances). The simulation box is ( − L, L ) with L = 120. This computation involves about 2,000eigenstates and 5,000 cores on Titan supercomputer, andtakes about 4 hours to reach convergence. The total andpolarization densities for the madness-hfb and hfb-ax simulations are shown in Fig. 3. The 3D pairing poten-tial is displayed in Fig. 4. The oscillations of the pair-ing field in a spin-polarized system, characteristic of theLarkin-Ovchinnikov phase, are clearly seen (see Ref. [28]for more discussion). FIG. 4. (Color online) The pairing potential of 100 particleswith P = 0 . η = 16 computedwith SLDA. B. Skyrme HF+BCS solver for nuclei
Most of the currently-envisioned applications of madness-hfb pertain to the nuclear many-body prob-lem. To this end, the adaptive multiresolution Skyrme-HFB solver has also been developed. The madness-hfb approach for nuclei is similar to the SLDA for cold atomsbut much more involved due to the continuum discretiza-tion, as the atomic nucleus is an open system and associ-ated boxes are large [35]. Therefore, as an initial step, wecarry out Skyrme HF and Skyrme HF+BCS calculationsand benchmark them with HFODD.For both HF and HF+BCS calculations, we con-sider the neutron-rich nucleus
Mo which is triaxi-ally deformed in its ground state in some models [36].We use SkM* [37] Skyrme parametrization, and take (cid:126) /2m=20.73 MeV fm for benchmarking purpose.In pairing calculations for Mo, due to the smallneutron separation energy, the positive-energy HF lev-els are important as they participate in pairing. Thiscreates a problem when trying to compare BCS or HFBresults based on solvers using coordinate-space frame-work and harmonic-oscillator expansion as the contin-uum representation is different in both approaches. In-deed, coordinate-space solvers madness-hf or madness-hf+bcs , when applied to large boxes, produce a verydense unbound single-neutron spectrum [11, 35, 38]. Onthe other hand, the single-neutron spectrum of oscillator-based hfodd is fairly sparse. Therefore, to minimizethe difference between these two codes for a meaning-ful benchmarking, we switch off neutron pairing and re-tain only bound 70 single-proton orbits in the BCS phasespace. We adopt mixed density dependent delta interac-tion [39]. The proton pairing strength is chosen to be −
500 MeV to obtain paired solution. Our madness-hf and madness-hf+bcs calculations are performed in alarge 3D box ( − L, L ) with L = 50 fm. The wavelet orderis p = 9 with requested truncation precision (cid:15) = 10 − . hfodd calculations are performed with 1140 and 1540spherical harmonic oscillator states, corresponding to 17and 19 shells, respectively. The oscillator constant is0.4975890 fm − .Since MADNESS calculations are numerically exten- sive, it is desirable to warm-start the runs with wavefunctions (or densities) from the converged hfodd solu-tion. We have implemented such an interface between hfodd and madness-hfb .Table III compares madness-hf and hfodd results forthe triaxial ground-state configuration in Mo.
TABLE III. Comparison between results of madness-hf and hfodd for the triaxial ground state of
Mo: total bindingenergy E t , kinetic energy, E kin , Coulomb energy E c , and spin-orbit energy E SO (all in MeV), mass r.m.s. radius R rms (infm), and mass quadrupole moments Q and Q (in fm ).The “0-th iter” column shows madness-hf warm-start num-bers at the beginning of the iteration process with wave func-tions and densities imported from converged hfodd resultsusing 1140 basis states. hfodd hfodd madness-hf madness-hf (1140) (1540) (0-th iter) (converged) E t − − − − E kin E c E SO − − − − R rms Q | Q | The madness-hf results labeled “0-th iter” are warm-start initialization numbers, with densities imported from hfodd (1140). As expected, “0-th iter” and hfodd (1140)values are extremely close. A very small difference ≈ ∼ ρ γ +2 ) produces the largest difference. The ex-cellent agreement between these two calculations indi-cates that the interface between the two solvers has beenimplemented correctly, and that the individual SkyrmeEDF terms have been coded properly in madness-hf .By increasing the basis size in hfodd to 1540 states, thetotal binding energy decreases by ∼
130 keV. However,it is still ∼
190 keV above the madness-hf result. Thisdifference can be traced back to asymptotic behavior ofnucleonic densities obtained in the two solvers. Figure. 5displays the neutron density profiles along x -, y -, and z -axis (moving from the origin) computed in hfodd (1140), hfodd (1540), and madness-hf . When displayed in lin-ear scale, one can hardly see a difference between hfodd and madness-hf predictions. However, when inspectingthe density in a logarithmic scale, one can see a charac-teristic damping at large distances (10-12 fm) in hfodd due to the finite size of oscillator basis. We recall thatthe madness-hf calculations were carried out in a boxextending to 50 fm.Finally, Table IV displays HF+BCS results. Again,the agreement between madness-hfb and hfodd is ex-cellent, with the total binding energy in madness-hfb being ∼
150 keV below that of hfodd (1540). x, y, z (fm) yx z l og ( ρ n ) x, y, z (fm) ρ n ( f m - ) HFODD(1140)HFODD(1540)MADNESSyxz Mo FIG. 5. (Color online) Neutron density distribution for
Mo in madness-hf (solid line), hfodd (1140) (dotted line),and hfodd (1540) (dashed line) along x -, y -, and z -axis, mov-ing from the origin. The inset (in a logarithmic scale) illus-trates the tail behavior of density.TABLE IV. Similar to table III, except that we include BCSpairing for protons, see text for details. hfodd hfodd madness-hf madness-hf (1140) (1540) (0-th iter) (converged) E t − − − − E pair − − − − λ p − − − − E kin E c E SO − − − − R rms Q | Q | V. SUMMARY
In this paper, we introduce nuclear DFT frameworkbased on the adaptive multi-resolution 3D HFB solver madness-hfb . The numerical method employs harmonicanalysis techniques with multiwavelet basis; user-definedfinite precision is guaranteed. The solver applies state-of-the-art in parallel programming techniques that cantake advantage of high performance supercomputers.Applications have been presented for polarized ul-tracold atoms in very elongated traps and for triaxialneutron-rich nuclei. The solver has been benchmarkedagainst other advanced HFB solvers: a 2D coordinate-space solver hfb-ax based on B-spline technique and a3D solver hfodd employing the harmonic oscillator basisexpansion. The advantage of madness-hfb is its abilityto treat large and complex systems without restrictionon symmetries. Examples of future nuclear structure ap-plications include: weakly bound nuclei with large spa-tial extensions, heavy-ion fusion, nuclear fission, complextopologies in super- and hyperheavy nuclei [40–42], andpasta phases in the inner crust of neutron stars [43–46].Future atomic applications of madness-hfb include de-scription of large number of fermions ( ∼ ) in highlyelongated optical traps ( η ∼
50) [47].
ACKNOWLEDGMENTS
Useful discussions with N. Hinohara, J. Sheikh, andN. Schunck are gratefully acknowledged. This work wassupported by the U.S. Department of Energy (DOE)under Contracts No. DE-AC05-00OR22725 (ORNL),No. DE-FG02-96ER40963 (University of Tennessee),No. de-sc0008499 (NUCLEI SciDAC Collaboration), No.DE-FG02-13ER42025 (China-U.S. Theory Institute forPhysics with Exotic Nuclei); by the National NaturalScience Foundation of China under Grant No.11375016,11235001. An award of computer time was providedby the National Institute for Computational Sciences(NICS) and the Innovative and Novel Computational Im-pact on Theory and Experiment (INCITE) program us-ing resources of the OLCF and ALCF facilities. [1] S. Bogner et al. , Comput. Phys. Commun. , 2235(2013).[2]
Nuclear Physics: Exploring the Heart of Matter. Reportof the Committee on the Assessment of and Outlook forNuclear Physics (The National Academies Press, 2012).[3] G. F. Bertsch, D. J. Dean, and W. Nazarewicz, SciDACReview , 42 (2007).[4] M. Stoitsov, J. Dobaczewski, W. Nazarewicz, andP. Ring, Comput. Phys. Commun. , 43 (2005).[5] M. V. Stoitsov, N. Schunck, M. Kortelainen, N. Michel,H. Nam, E. Olsen, J. Sarich, and S. Wild, Comput. Phys.Commun. , 1592 (2013). [6] J. Dobaczewski and P. Olbratowski, Comput. Phys. Com-mun. , 158 (2004).[7] N. Schunck, J. Dobaczewski, J. McDonnell, W. Satu(cid:32)la,J. Sheikh, A. Staszczak, M. Stoitsov, and P. Toivanen,Comput. Phys. Commun. , 166 (2012).[8] N. Michel, K. Matsuyanagi, and M. Stoitsov, Phys. Rev.C , 044319 (2008).[9] M. Kortelainen, J. McDonnell, W. Nazarewicz, P.-G.Reinhard, J. Sarich, N. Schunck, M. V. Stoitsov, andS. M. Wild, Phys. Rev. C , 024304 (2012).[10] J. Dobaczewski, H. Flocard, and J. Treiner, Nucl. Phys. A422 (1984). [11] J. Dobaczewski, W. Nazarewicz, T. R. Werner, J. F.Berger, C. R. Chinn, and J. Decharg´e, Phys. Rev. C , 2809 (1996).[12] J. C. Pei, M. V. Stoitsov, G. I. Fann, W. Nazarewicz,N. Schunck, and F. R. Xu, Phys. Rev. C , 064306(2008).[13] E. Ter´an, V. E. Oberacker, and A. S. Umar, Phys. Rev.C , 064314 (2003).[14] I. Stetcu, A. Bulgac, P. Magierski, and K. J. Roche,Phys. Rev. C , 051309 (2011).[15] Y. Hashimoto, Phys. Rev. C , 034307 (2013).[16] A. Bulgac and M. M. Forbes, Phys. Rev. C , 051301(2013).[17] B. Alpert, SIAM J. Math. Anal. , 264 (1993).[18] MADNESS website:, http://code.google.com/p/m-a-d-n-e-s-s/.[19] S. G. Mallat, Trans. Amer. Math. Soc. , 69 (1989).[20] B. Alpert, G. Beylkin, D. Gines, and L. Vozovoi, J.Comp. Phys. , 149 (2002).[21] G. Beylkin and M. J. Mohlenkamp, Proc. Natl. Acad.Sci. USA , 10246 (2002).[22] R. J. Harrison, G. I. Fann, T. Yanai, Z. Gan, andG. Beylkin, J. Chem. Phys. , 11587 (2004).[23] I. Babuˇska and B. Q. Guo, Adv. Eng. Softw. , 159(1992).[24] G. I. Fann, R. J. Harrison, and G. Beylkin, J. Phys.:Conf. Ser. , 461 (2005).[25] A. Bulgac and M. M. Forbes, arXiv:0808.1436 (2008).[26] A. Bulgac and M. M. Forbes, Phys. Rev. Lett. ,215301 (2008).[27] G. Bertsch, J. Dobaczewski, W. Nazarewicz, and J. Pei,Phys. Rev. A , 043602 (2009).[28] J. C. Pei, J. Dukelsky, and W. Nazarewicz, Phys. Rev.A , 021603 (2010).[29] K. Bennaceur and J. Dobaczewski, Comput. Phys. Com-mun. , 96 (2005). [30] A. Bulgac and K. J. Roche, J. Phys.: Conf. Ser. ,012064 (2008).[31] G. I. Fann, J. Pei, R. J. Harrison, J. Jia, J. Hill, M. Ou,W. Nazarewicz, W. A. Shelton, and N. Schunck, J.Phys.: Conf. Ser. , 012080 (2009).[32] P. Fulde and R. A. Ferrell, Phys. Rev. , A550 (1964).[33] A. I. Larkin and Y. N. Ovchinnikov, Sov. Phys. JETP , 762 (1965).[34] A. Bulgac, Phys. Rev. A , 040502 (2007).[35] J. C. Pei, A. T. Kruppa, and W. Nazarewicz, Phys. Rev.C , 024311 (2011).[36] Y. Shi, C. L. Zhang, J. Dobaczewski, andW. Nazarewicz, Phys. Rev. C , 034311 (2013).[37] J. Bartel, P. Quentin, M. Brack, C. Guet, and H.-B.H˚akansson, Nucl. Phys. A , 79 (1982).[38] P. J. Borycki, J. Dobaczewski, W. Nazarewicz, and M. V.Stoitsov, Phys. Rev. C , 044319 (2006).[39] J. Dobaczewski, W. Nazarewicz, and M. V. Stoitsov,Eur. Phys. J. A , 21 (2002).[40] W. Nazarewicz, M. Bender, S. ´Cwiok, P. Heenen,A. Kruppa, P.-G. Reinhard, and T. Vertse, Nucl. Phys.A , 165 (2002).[41] J. Decharg´e, J.-F. Berger, M. Girod, and K. Dietrich,Nucl. Phys. A , 55 (2003).[42] P. Jachimowicz, M. Kowal, and J. Skalski, Phys. Rev. C , 054302 (2011).[43] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys.Rev. Lett. , 2066 (1983).[44] W. G. Newton and J. R. Stone, Phys. Rev. C , 055801(2009).[45] C. O. Dorso, P. A. Gim´enez Molinelli, and J. A. L´opez,Phys. Rev. C , 055805 (2012).[46] A. S. Schneider, C. J. Horowitz, J. Hughto, and D. K.Berry, Phys. Rev. C , 065807 (2013).[47] G. B. Partridge, W. Li, R. I. Kamar, Y.-A. Liao, andR. G. Hulet, Science311