Stochastic Estimation of Nuclear Level Density in the Nuclear Shell Model: An Application to Parity-Dependent Level Density in $^{58}$Ni

We introduce a novel method to obtain level densities in large-scale shell-model calculations. Our method is a stochastic estimation of eigenvalue count based on a shifted Krylov-subspace method, which enables us to obtain level densities of huge Hamiltonian matrices. This framework leads to a successful description of both low-lying spectroscopy and the experimentally observed equilibration of $J^\pi=2^+$ and $2^-$ states in $^{58}$Ni in a unified manner.


Introduction
Nuclear level density, which constitutes a very large uncertainty in estimating cross sections of compound nuclear reactions, plays an essential role in describing those reactions usually with the Hauser-Feshbach theory [1]. Systematic data for compound nuclear reactions, including (n, γ) reaction data, are highly needed in many applications including nuclear astrophysics and nuclear engineering. The (n, γ) cross sections for short-lived nuclei, however, cannot be directly measured, and therefore an accurate estimate of nuclear level density is crucial for those needs [2]. While nuclear level density is known to well follow phenomenological formulas, such as the backshifted Fermi-gas model and the constant temperature model, their model parameters depend on nuclear structure and are usually derived from experimental data [3]. Moreover, the phenomenological models assume rather simple spin-parity dependence of level density, whose validity needs to be verified. Microscopic theories of nuclear level density should thus be developed.
The large-scale shell-model calculation provides, in principle, quite reliable level density because this method suitably takes two-body correlations into account, thus well describing level structures. Although the Lanczos diagonalization [4] now makes possible the calculation of low-lying levels for Hamiltonian matrices up to O(10 10 ) dimensions in an M-scheme [5], its applicability to the level density is much more limited because high-lying states are very slow to converge in the Lanczos iterations. Much effort has been paid to overcoming this limitation (e.g. the method based on the moment of the matrix [6,7,8]). The shell-model Monte Carlo (SMMC) method has been proposed and developed [9,10], but it is rather difficult for the SMMC to use realistic residual interactions which cause the sign problem. Recently several realistic interactions have succeeded in systematically describing lowlying levels including those of exotic nuclei [5,11,12,13,14,15,16]. It is thus quite interesting to investigate whether a unified description of low-lying levels and level density can be obtained with such modern realistic interactions [17].
In this Letter, we apply a novel method of counting eigenvalues [18] to the shell-model calculation of level density, and demonstrate its feasibility and usefulness. This method is applicable to any effective interaction and it straightforwardly provides spin-parity dependent level densities. More importantly, its computational cost is almost at the same order as that of the usual Lanczos diagonalization for low-lying states. Taking these advantages, we are successful in calculating spin-parity dependent level densities in 58 Ni with a realistic effective interaction, and in resolving a puzzle that previous microscopic calculations failed to resolve in reproducing no parity dependence of the 2 + and 2 − level densities observed recently [19]. The M-scheme dimension for the 2 − levels reaches 1.5 × 10 10 , which is nearly the current limit of the Lanczos diagonalization.

Formalism
We first outline a method for stochastically estimating eigenvalue distribution proposed recently [18]. In nuclear shell model calculations, a manybody wave function is described as a linear combination of a huge amount of many-body configurations, which are called the M-scheme basis states [5]. The shell-model energy is obtained as an eigenvalue of the M-scheme shellmodel Hamiltonian matrix, H. The number of eigenvalues µ k in a specified energy range E (k−1) < E < E (k) is obtained by the contour integral Γ k on the complex plane in Fig.1 as The contour integral on Γ k is numerically obtained by discretizing the contour line with mesh points z ...  The trace of the inverse of matrix in Eq.(1) is stochastically estimated by sampling dozens of random vectors from the whole Hilbert space. An unbiased estimation is given by where N s is the number of sample vectors and v s are vectors whose elements take 1 or −1 with equal probability randomly [20]. Note that the norm of v s equals the whole number of states. Since the Hamiltonian matrix is written in the M-scheme representation, whose basis states are the eigenstates of the z-component of spin and parity, we obtain the total level density with a fixed parity at once. The remaining task is to obtain v T s (z In preceding works, the recursion method with the Lanczos algorithm with was used to compute this value [21], but the Lanczos algorithm requires reorthogonalization procedure [22]. In the present work, we adopt the complex orthogonal conjugate gradient (COCG) method [23], which does not need reorthogonalization.
We briefly explain how to calculate the v T s (z In actual calculations, x (s) (s = 1, 2, ..., N s ) are obtained simultaneously with a block algorithm (BCOCG) [24] with economical implementation for computing block bilinear form [25] for efficient computation.
Second, we solve v s = (z j . If we perform the COCG method for each z (k) j individually, an impractically large amount of computation is required. However, the shifted Krylov-subspace method [26,27] enables us to solve these equations at once. The basic idea behind this method is the invariance property of the shifted Krylov subspace: The Krylov subspace of the matrix z − H and vector v, which is defined as is independent of z. Since the COCG method is one of the Krylov-subspace methods, its solution of z = z ref is described in this subspace. The solution of any z is also obtained in the same subspace simultaneously, if n is large enough. Thus, by utilizing the COCG solution at z = z ref , the shifted method allows us to avoid the time-consuming calculations of matrix-vector product for any z (k) j , and to obtain the level density in the whole energy region at once. In practice, we adopt the shifted block CG-rQ (SBCG-rQ) method [24] for efficient computation. This computation demands memory usage for only three vector blocks (x n , r n , and p n ) and a working area for the matrix-vector products, while all eigenvectors are stored in the Lanczos method.
For spin-dependent (J-dependent) level density, we replace v s by v J s = P J v s , where P J is the J-projection operator. This projection is actually realized by filtering out a specified range of the eigenvalue of the J 2 operator by the Sakurai-Sugiura method [22,28]. Note that the J-projection is not necessary during the COCG iteration in this work, whereas J-projection is performed at every iteration to suppress numerical errors in the J-projected Lanczos method [5,22]. This is because an undesirable-J component caused by numerical error does not affect the inner product (v J s ) T x σ n . Even using such a sophisticated method, the eigenvalue count of a huge matrix requires high-performance computing. We combine the nuclear shell-model code "KSHELL" [29] and the eigenvalue-solver library "z-Pares" [30], which enable us to utilize the latest supercomputer efficiently.
It is worth noting a pioneering work for the level density using the Lanczos method by estimating the traces with random vectors in Ref. [31].

Benchmark test
We demonstrate how this stochastic estimation works well in a small system. Figure 2 shows the level density of 28 Si with the USD interaction [32]. The level density is obtained with N s = 32 and a 200 keV energy bin. We compare the present estimation with the exact values, which are provided by counting the 1000 lowest states calculated with the Lanczos method. The present estimation successfully reproduces the exact values in 8% error at around E x = 15 MeV with such a small bin. Non-empirical error evaluation remains as future work [40]. Figure 3 (a) shows the residual |r n |/|v s | of the BCOCG against the iteration number in the 28 Si case. The residual rapidly converges to zero at low energies of practical interest. However, the convergence becomes worse in high excited energies as the level density increases. While a similar tendency is also seen for the resultant level densities in Fig. 3 (b), their convergence is much faster than the residual and is quite stable. For example, only 28 iterations are required for the convergence of the level density at E x =10 MeV.

Parity-equilibration of 58 Ni level density
Now that the feasibility of the present method has been confirmed, we can move on to its application. Here we will demonstrate that the 2 + and 2 − level  densities in 58 Ni [19] are excellently reproduced with the present method. To describe both parity states of nuclei around 58 Ni with the shell model in a practical way, we take the 0hω and 1hω states in the full sd+pf +sdg valence shell for natural-and unnatural-parity states, respectively. Since 58 Ni is located in the middle of the pf shell, higherhω states are expected to be dominant in relatively high excitation energies, which will be confirmed later. The full 1hω calculation is still impractical for 58 Ni, and we truncate the basis states by excluding configurations involving more than six-particle excitations from the 0f 7/2 filling configurations. We have confirmed that this truncation has a minor effect on the level density in the 0hω calculation. The resulting M-scheme dimension for the 2 − states in 58 Ni, 1.5 × 10 10 , is much beyond the capability of the direct counting of eigenstates by the Lanczos diagonalization, but within the scope of the present method. We remove the spurious center-of-mass contamination associated with 1hω excitation by using the Lawson method [33], whereas this removal is usually not included in SMMC. The spurious states are lifted up to E x ∼ 600 MeV, thereby separated clearly with the Hamiltonian H ′ = H + βH CM with βhω/A = 10 MeV. The effective interaction taken in this study is a natural extension of the SDPF-MU interaction for the sd-pf shell [16] to the one for the sd-pf -sdg shell. Namely, exactly in the same way as the SDPF-MU, the present interaction consists of the USD interaction [32] for the sd shell, the GXPF1B interaction [12] for the pf shell, and a refined monopole-based universal interaction (V MU ) [34,35] for the rest of the two-body matrix elements. It is noted that the SDPF-MU interaction well describes low-lying 1hω levels [36,37] as well as 0hω states [16]. Going back to the present sd-pf -sdg-shell Hamiltonian, the single-particle energies (SPEs) of the sdg orbits are left to be determined, while those of the sd-pf orbits are taken from the SDPF-MU interaction. Here, the SPEs for the 0g 9/2 , 1d 5/2 , and 2s 1/2 orbits are fixed to reproduce experimental spectroscopic strengths around 58 Ni, as shown in the next paragraph. On the other hand, experimental information on the remaining upper sdg orbits is missing near 58 Ni. These SPEs are chosen to fit the empirical neutron effective SPEs on top of 90 Zr [34].
The 1hω shell gaps (the sd-pf and pf -sdg shell gaps) dominate the overall positions of the 1hω states and thus, their level densities. The reliability of those shell gaps in the present Hamiltonian is examined by comparing theoretical and experimental spectroscopic strengths (C 2 S) for the one-proton (±1p) and one-neutron (±1n) transfer reactions in 58 Ni. The results are summarized in Table 1, where the sd-pf shell gaps and the pf -sdg shell gaps are probed by the −1p and −1n reactions ( 57 Co and 57 Ni) and the +1p and +1n reactions ( 59 Cu and 59 Ni), respectively. Aside from some difference in C 2 S fragmentation for the 3/2 + levels in 57 Ni and the 1/2 + levels in 59 Ni, the present shell-model results agree well with the experimental data, thus confirming suitable effective SPEs for upper sd (1s 1/2 and 0d 3/2 ) and lower sdg (0g 9/2 , 1d 5/2 and 2s 1/2 ) orbits. In addition, the 3 − 1 level in 58 Ni is also well reproduced: 4.55 MeV (Cal.) vs. 4.47 MeV (Exp.).
By utilizing this realistic interaction and the new stochastic method for estimating level density, we can calculate 2 + and 2 − level densities in 58 Ni. The 2 − level densities are calculated with N s = 16 and 1000 BCOCG-iterations which are large enough to reach reasonable convergence. Figure 4 compares 2 ± level densities in 58 Ni between the experimental data [19] and the present calculation. What is surprising in this experimental data is that the equilibration of 2 + and 2 − levels occurs already at E x ∼ 8 MeV at variance with SMMC and HFB-based [39] estimates in which the equilibration is realized  Table 1: Excitation energies and one-nucleon spectroscopic factors (C 2 S) of the ground states and low-lying unnatural parity states in 57 Ni, 59 Ni, 57 Co, and 59 Cu compared between the shell-model results and experiments [38]. The C 2 S values are for one-nucleon transfer reactions of 58 Ni, and the transferred nucleon and its orbital are specified in the fifth column.
only at around 20 MeV [19]. Both of the calculations overestimate the 2 + level densities and underestimate the 2 − level densities in low excitation energies. In contrast, the present calculation correctly reproduces the early onset of the equilibration on the basis of good agreement of the 2 + and 2 − level densities there. It is worth pointing out that nucleon excitation from the sd shell accounts for about half of the 2 − level densities. With increasing excitation energy, the present calculation begins to underestimate the

Summary
In summary, we introduced a stochastic estimation of the level density in nuclear shell-model calculations and demonstrated its feasibility. This method is based on the shifted Krylov-subspace method and enables us to obtain the level density of a huge-dimension matrix beyond the limitation of the conventional Lanczos method. The contamination of spurious center-of-mass excitation is clearly removed by the Lawson method. By combining proven SDPF-MU and V MU interactions, we constructed a realistic effective interaction which successfully describes low-lying levels and their spectroscopic factors around 58 Ni. With this realistic interaction and the stochastic method, we obtained the level densities of J π = 2 + and J π = 2 − in 58 Ni. These densities are in excellent agreement with the experimental results and show the equilibration at E x ≥ 8 MeV, whereas the preceding microscopic calculations showed strong parity dependence. The present framework bridges the lowlying spectroscopy and microscopic understanding in statistical region; thus, further studies in this direction should be quite promising. Moreover, the application to ab initio shell-model calculations in light nuclei is expected.
15K05094) from JSPS, Japan. The numerical calculations were performed on K computer at RIKEN AICS (hp140210, hp150224) and COMA supercomputer at the University of Tsukuba.