Stochastic estimation of level density in nuclear shell-model calculations

An estimation method of the nuclear level density stochastically based on nuclear shell-model calculations is introduced. In order to count the number of the eigenvalues of the shell-model Hamiltonian matrix, we perform the contour integral of the matrix element of a resolvent. The shifted block Krylov subspace method enables us its efficient computation. Utilizing this method, the contamination of center-of-mass motion is clearly removed.


Introduction
Nuclear level density is a key ingredient of the Hauser-Feshbach calculations [1] to understand neutron-capture processes and compound nucleus microscopically.Theoretically, some phenomenological formulas for level density are well known such as the backshifted Fermi gas model [2].However, microscopic calculation of the nuclear level density is prerequisite for the precise prediction of r-process nuclei and for further microscopic understanding of the statistical properties of nuclei.
Nuclear shell-model calculation is one of the most powerful theoretical frameworks to evaluate the level density.However, the rapid growth in the dimension of the shell-model Hamiltonian matrix remains an inevitable obstacle.An efficient numerical method to solve the eigenvalue problem of the huge sparse matrix and to estimate its eigenvalue distribution has been demanded.Much effort has been paid to develop efficient methods to estimate nuclear level densities based on shell model, e.g. the shell model Monte Carlo (SMMC) [3], the moment-based methods [4,5], and the Lanczos-based method [6].In this proceedings, we review stochastic estimation of the eigenvalue distribution based on shifted Krylov-subspace method [7] and its application to nuclear shell model calculations [8].We also discuss its validity and feasibility.
The Hamiltonian in nuclear shell-model calculations is written in second quantization as where c † i denotes the creation operator of the single particle state i.The many-body wave function is described as a linear combination of a vast number of Slater determinants, which are the antisymmetrized products of the single-particle wave functions.The simplest representation for the Slater determinant is called "M-scheme" basis state: where A and |− are the number of active nucleons and an inert core, respectively.M (α) i denotes the single-particle state which is occupied by the α-th active nucleon.Thus, M i = {M (1)  i , M (2)  i , ..., M (A) i } is called "configuration" and means that 1st, 2nd, ..., A-th particles occupy M (1)  i , M (2)  i , ..., M (A) i singleparticle states, respectively.
Since the model space is fully spanned by the M-scheme basis states, the shell-model wave function is expressed as a linear combination of them: where D is the number of the M-scheme basis states.The Schrödinger equation is transformed into the eigenvalue problem of the Hamiltonian matrix H i j , where the matrix H i j = M i |H|M j is real symmetric.There are two symmetries which the Hamiltonian matrix in the M-scheme basis explicitly has: invariant under rotational symmetry around z-axis and parity inversion.The operators concerning these symmetry are referred by J z and Π and the corresponding quantum numbers are M and π.Owing to these symmetries, the Hamiltonian matrix is a block-diagonal form.We need to solve only a block matrix specified by the M and π due to the degeneracy.The largest dimension of this block matrix is often called the M-scheme dimension.Moreover, there is an additional symmetry to be conserved, namely total angular momentum J 2 .However, it does not contribute to the block form of the matrix since the M-scheme state is not a good eigenstate of J 2 .
Furthermore, this block matrix is quite sparse, namely the fraction of the non-zero matrix elements is quite small, since the shell-model Hamiltonian in Eq.( 3) consists only of one-body and two-body interactions.Figure 1 shows the non-zero matrix elements of the Hamiltonian matrix of 44 Ti with p f shell and M π = 0 + .Its M-scheme dimension is 4000.In this case, the fraction of the non-zero matrix elements is around 5%.The fraction is expected to be small in a larger system.
Owing to this sparsity, the matrix-vector product is performed efficiently with explicitly avoiding the treatment of the zero matrix elements, and the Krylov-subspace method is quite advantageous to  solve the eigenvalue problem.The Krylov subspace is spanned by an initial vector, u, and its products with the matrix H to the power: In usual shell-model calculations the Lanczos method, one of the Krylov-subspace methods, is used to obtain low-lying eigenvalues.The convergence of the eigenvalue as a function of the number of Lanczos iteration, n, is shown in Fig. 2. The y axis shows the eigenvalues of the Hamiltonian matrix in the Krylov subspace.While the lowest eigenvalue converges quite fast, the high-lying eigenvalue converges slower as the eigenvalue goes higher.At a rough guess, the current limitation of the Lanczos method is that the product of the M-scheme dimension and the number of eigenvectors does not exceed O(10 10 ).While the largest dimension for a few low-lying states are O(10 10 ), the limitation for the calculation of the level density is much strict: since it requires O(10 3 ) eigenvectors, the largest dimension becomes O(10 7 ).The largest dimension is extended up to O(10 10 ) by utilizing our proposed method of stochastic estimation, the detail of which is described in the following section.

Framework of the stochastic estimation of the eigenvalue count
We briefly review the stochastic estimation of eigenvalue density of sparse matrix in this section.This method was introduced and demonstrated in condensed matter physics [7].
For the estimation of the level density, we count the number of the eigenvalues in a specified energy region, Γ.It is given as μ = Tr(P Γ ), where P Γ is a projection operator to the subspace spanned by the eigenvectors whose eigenvalues are inside Γ.It is realized by the Cauchy integral as where Γ is a contour to surround a certain energy region, or an oval Γ 1 in Fig. 3.Although we set Γ i for each energy bin for evaluating the level density, we omit the index of energy bin i for simplicity.The contour integral is approximated by a summation of the discretized points z j shown as blue symbols in Fig. 3 with the corresponding weights w j .Note that this projector also appears in the Sakurai-Sugiura method, which is an eigenvalue solver [10,11].However, since the trace of Eq.( 6) is difficult to calculate directly, we estimate it utilizing Hutchinson's estimator [12] as the following equation: where u s is a sample vector whose matrix elements are taken as −1 or 1 randomly.N s is the number of sample vectors.While this estimation contains stochastic error, the magnitude of the error depends on the property of P Γ .For example, if P Γ is a diagonal matrix, the stochastic error becomes zero.
In preceding works, this estimator gives sufficiently small error even with a small number of sample vectors [7].In this work, we take N s = 16.The next step is to evaluate u T s 1 z−H u s by solving a linear equation (z−H)x = u s .Since we have multiple u s , we adopt the block conjugate gradient (BCG) method to solve the multiple linear equations (z − H)X = V where V is the D × N s matrix whose column is u s .We adopt block CG-rQ (BCG-rQ) method [14] for efficient computation.
With A = z − H, the BCG-rQ algorithm is shown in Algorithm 1 where qr(V) denotes QR decomposition of V.The X k , P k , Q k , are also D × N s matrices.The residual matrix Q k is orthonormalized by QR decomposition every iteration to improve the numerical stability [13].The block algorithm accelerates the convergence and also helps to use a CPU efficiently by increasing contiguous memory access.
However, it requires enormous computational resources to perform the BCG-rQ calculation for each z j individually.This problem is overcome by utilizing the shift invariance of the Krylov subspace.This invariance means that the Krylov subspace of the matrix H is the same as that of z − H:  calculations, the model space is taken as the p f shell and the GXPF1A interaction [9] is adopted for the realistic effective interaction.This interaction was constructed based on G-matrix interaction and chi-square fitted for the experimental data of low-lying states.The M-scheme dimension of 56 Fe reaches 501,113,392.The blue symbols with error bars in Fig. 4 show the experimental level density of 56 Fe obtained by the Oslo method [17].The red line shows the exact shell-model level density of 56 Fe.The exact values are obtained by the direct counting of the 100 lowest states.The energy bin is ΔE = 0.2 MeV, and it is the same as the experimental one.The exact Lanczos value shows good agreement the experimental one up to 6 MeV.However, it becomes rapidly difficult to obtain the exact value in higher-energy region since the level density, and hence the number of eigenvalues to be evaluated, increase exponentially.
The black line in Fig. 4 is obtained by the present stochastic estimation with N s = 16 and 1550 BCG iteration.This estimation reproduces well the exact value within a certain stochastic error.The estimation is feasible beyond the Lanczos limitation and it reproduces well the experimental value up to 10 MeV.The deviation around E x ∼ 9 MeV is considered to be the contribution of negative-parity states, which is beyond the present model space.As the excitation energy increases, the ratio of the stochastic error to the level density would decrease since the level density increases exponentially and its error is expected to be proportional to the square root of the level density.Further discussion of the error estimation remains for future work.B. A. Brown and A. C. Larsen also showed the shell-model results of the level density of 56 Fe by the direct counting of the Lanczos results in Ref. [18].They used the same model space and interaction, but the truncation scheme which they applied in order to save computational resources brought about modest underestimation.

Removal of the contamination of the center-of-mass motion
Since an atomic nucleus is an isolated system and translational invariance should be taken into account the contamination of center-of-mass motion should be removed in nuclear shell model calculations beyond 0 ω model space [19].In order to remove this contamination, we replace the Hamiltonian by with β ω/A = 10 MeV.H CM is the Hamiltonian of the center-of-mass motion in a harmonic oscillator potential whose energy quanta is ω.We take β large enough so that the expectation value of H CM − 3 2 ω is close to zero and the excitation of the center-of-mass motion is suppressed.The contaminated states are lifted up to high-energy region.This prescription is called the Lawson method [20].
Figure 5 shows the level density of J π = 1 − states of 48 Ca.The stochastic estimation in shell model is performed with sd-p f -sdg shells and 1 ω truncation.The level densities are obtained in 1 MeV energy bin.The realistic interaction presently used was introduced in Ref. [21] and provides us with a good description of collective states such as the 3 − 1 state and the giant dipole resonances of Ca isotopes [22].The exact value obtained by the Lanczos method is well reproduced by the present estimation with N s = 16 and 500 BCG iterations.The estimation with N s = 4 and 100 BCG iterations is shown by the dotted line.Since the number of iterations is small and the BCG does not converge well, it shows pseudo oscillation and disagreement with the exact one.
Figure 6 shows the high-energy region of Fig. 5.The spurious states are lifted up to E x ∼ 520 MeV, corresponding a small peak in Fig. 6, thereby separated clearly from the main peak.Such separation is also clearly seen even in the dotted line, which is obtained before the convergence.

Summary
We proposed a novel stochastic estimation method for calculation level densities in nuclear shell model.This method enables us to estimate level density in a rather small energy bin, such as 200 keV as demonstrated in the level density of 56 Fe.It is also demonstrated that the contamination of the CNR * 15   2003  center-of-mass motion is clearly removed by the Lawson prescription.The present method is feasible up to O(10 10 ) M-scheme dimension, which is a similar order to the limitation to obtain some lowlying states by the Lanczos method.By adopting this framework, we successfully reproduced the parity equilibration of 2 + and 2 − level densities of 58 Ni in low-energy region in Ref. [8].Further applications of this method are in progress.

Figure 1 .
Figure 1.M-scheme Hamiltonian matrix of 44 Ti with p f shell and M = 0.The points represent non-zero matrix elements.Only upper-right part is shown since it is a real symmetric matrix.

Figure 2 .
Figure 2. Convergence of the Lanczos method against the number of iteration, n.The symbols denote the shell-model energies of J π = 0 + 1 (black circles) to J π = 0 + 10 (yellow circles) states of 56 Ni with p f shell and GXPF1A interaction [9].

Figure 3 .
Figure 3. Schematic view of the contour line in complex plane z and its discretized points for the numerical calculations.The red crosses show eigenvalues of H.

Figure 4 .
Figure 4. Level densities of 56 Fe against the excitation energy.These are obtained by stochastic estimation (black line) and by direct counting of the Lanczos results (red line).The experimental values (blue symbols with error bars) are taken from [17].The energy bin is 0.2 MeV.

Figure 5 .¢Figure 6 .
Figure 5. Level densities of J π = 1 − states of 48 Ca against the excitation energy.These are obtained by the stochastic estimation with N s = 4 and 100 CG iterations (black dashed line), that with N s = 16 and 500 CG iterations (blue solid line), and the exact values obtained by direct counting of the Lanczos results (red dashed line).The energy bin is 1 MeV. )