# Microscopic calculations of nuclear level densities with the Lanczos method

W. E. Ormand  Lawrence Livermore National Laboratory, P.O. Box 808, L-414,
Livermore, California 94551, USA  Department of Physics and the National
Superconducting Cyclotron Laboratory,
Michigan State University, East Lansing, MI 42284-1321, USA  B. A. Brown
Department of Physics and the National Superconducting Cyclotron Laboratory,
Michigan State University, East Lansing, MI 42284-1321, USA

###### Abstract

A new method for computing the density of states in nuclei making use of an
extrapolated form of the tri-diagonal matrix obtained from the Lanczos method
is presented. It will be shown that the global, average properties of the
entire Lanczos matrix can be predicted from just four Lanczos iterations. The
extrapolated Lanczos matrix (ELM) approach provides for an accurate
computation of the density of states described within the configuration space,
which, in some cases, is sufficient to accurately calculate the density of
states at, or near, the neutron separation energy. Comparisons between theory
and experiment are shown for 57Fe, 74Ge, and 76Ge. In addition, we show
results for the $J$-dependence of moments and the level density for these
three nuclei.

###### pacs:

21.10.Ma,21.60.Cs,27.40.$+$z

##  I Introduction

The density of states is a fundamental property of nuclear structure and plays
a key role in nuclear reactions. An important example is the radiative capture
of neutrons on short-lived nuclei, which, through the r-process r-process  in
supernovae and/or neutron-star mergers merg , are thought to be responsible
for the synthesis of the elements heavier than iron. Ideally, these reactions
can be measured or constrained by experiment. Unfortunately, in most cases,
the target nuclei are so short lived that direct measurement is not possible,
and the only alternative is to rely on theoretical calculations or indirect
measurements such as surrogates surr , which are themselves reliant on
theoretical input.

Nuclear reaction approaches such as Hauser-Feshbach Hauser-Feshbach  can give
an accurate description of the neutron-capture cross section. However, the
Hauser-Feshbach model requires accurate knowledge of the density of states up
to the neutron-decay threshold. A challenge in nuclear theory is to accurately
compute the density of states. This is difficult because of the sheer number
of levels and configurations and the strong nature of the nuclear Hamiltonian.
One microscopic approach is to account for correlations at the Hartree-Fock
level and to “count” non-interacting levels within the corresponding mean-
field single-particle space Goriely . Another is to use the Shell-Model Monte
Carlo (SMMC) AFMC ; AFMC-2 , which utilizes auxiliary fields to compute the
thermal trace for the energy, from which, the density of states can be
extracted from the inverse Laplace transform of the partition function AFMC-
rho . A limitation of the SMMC is the sign problem, which primarily limits the
approach to schematic interactions AFMC-2 . Moments methods, derived from
random matrix theory and statistical spectroscopy, can be used to construct
spin and parity dependent level densities for realistic Hamiltonians Mon75 ;
mom ; Horoi . Moments method, however, have been limited by the ability to
compute higher moments of the Hamiltonian, the overall structural form density
of states, and must be matched to the exact energies for low-lying states. The
stochastic estimation method shi  has a computational cost that is almost the
same order as the Lanczos method used here and requires a special computer
code to apply the shifted Krylov-subspace method 26 ; 27 .

In this article, we report on a new framework to provide an accurate
description of the statistical properties of a model Hamiltonian. Our specific
application is the calculation of the nuclear density of states within the
configuration-interaction approach using fully realistic nuclear Hamiltonians.
From universal properties of the Lanczos algorithm, we will demonstrate that
the first eight moments of the Hamiltonian can be obtained from just four
Lanczos iterations, which, in turn, can provide an accurate description of the
averaged, or global, properties of the nuclear system within the defined
Hilbert space. Several procedures to extract the density of states for model
Hamiltonians are presented here: 1) extrapolating the tri-diagonal Lanczos
matrix well beyond what is computationally viable, leading to an extrapolated
Lanczos method (ELM) to efficiently compute compute the density of states
within the configuration-interaction method; 2) an analytic continuation of
the ELM method; and 3) an approximation of the level density based on the
binomial distribution.

##  II Nuclear Structure Model

The principal goal behind nuclear-structure models is to find energy
eigenvalues and wave functions for the nuclear Hamiltonian within a well-
defined Hilbert space. In the nuclear shell model shell-model , or
configuration interaction, the Hilbert space is defined by a set of orbits,
usually denoted by the principal quantum number $n$, orbital angular momentum
$l$, and angular momentum $j$. The nuclear wave functions are constructed
through a set of basis states obtained by filling these orbits following the
Pauli principle. The basis states can consist of a set of Slater determinants
with well defined $z$-projection of angular momentum, $J_{z}=M$, in the so-
called $M$-scheme, or by projecting angular momentum (and possibly isospin)
onto the $M$-scheme Slater determinants. The $N$ many-body basis states,
$|\psi_{i}\rangle$, spanning the Hilbert space are used to construct the full
solution, i.e., $|\Psi\rangle=\sum_{i}c_{i}|\psi_{i}\rangle$. The coefficients
$c_{i}$ are found by computing the matrix elements of the Hamiltonian,
$H_{ij}=\langle\psi_{i}|\hat{H}|\psi_{j}\rangle$, and diagonalizing the
resulting Hermitian matrix. One of the most effective methods to find the
lowest eigenvalues is the Lanczos algorithm Lanczos , which starts with an
arbitrary vector $|v_{1}\rangle$ in the Hilbert space, and through successive
operations of $\hat{H}$, the matrix H is transformed into tri-diagonal form.
The first three terms are

$\displaystyle\hat{H}|v_{1}\rangle$
$\displaystyle=\alpha_{1}|v_{1}\rangle+\beta_{1}|v_{2}\rangle,$
$\displaystyle\hat{H}|v_{2}\rangle$
$\displaystyle=\beta_{1}|v_{1}\rangle+\alpha_{2}|v_{2}\rangle+\beta_{2}|v_{3}\rangle,$
$\displaystyle\hat{H}|v_{3}\rangle$ $\displaystyle=\hskip
39.12253pt\beta_{2}|v_{2}\rangle+\alpha_{3}|v_{3}\rangle+\beta_{3}|v_{4}\rangle,$
(1)

and the $|v_{i}\rangle$ form an orthonormal set. In practice this amounts to
applying $\hat{H}$ to the Lanczos vectors, and extracting the matrix elements
through subsequent dot-product operations and reorthogonalization, e.g.,
$\alpha_{1}=\langle v_{1}|\hat{H}|v_{1}\rangle$, and $\beta_{1}^{2}=\langle
v_{1}|(\hat{H}^{\dagger}-\alpha_{1})(\hat{H}-\alpha_{1})|v_{1}\rangle$ (note
that the phase of any of the $\beta_{i}$ is arbitrary). The power of the
Lanczos algorithm is that following successive applications of $\hat{H}$
(iterations), the eigenvalues of the tri-diagonal matrix quickly converge to
the extreme eigenvalues of the full matrix. Typically, the lowest energy in
the model space, $E_{0}$, is obtained in approximately 30 iterations
regardless the matrix dimension.

Of particular interest is the behavior of the tri-diagonal matrix elements
with increasing iterations. After several iterations, the diagonal elements,
$\alpha_{i}$, are roughly constant and nearly equal to the first moment
$H_{1}=\frac{1}{N}{\rm{Tr}}[\hat{H}]=\frac{1}{N}\sum_{i}H_{ii}$. At the same
time, the off-diagonal elements, $\beta_{i}$, generally decrease to zero as
$i\rightarrow N$, and exhibit a Gaussian-like behavior zuker .

In this work, we will examine the level density for selected Cr, Fe, and Ge
isotopes within the framework of the nuclear shell model. All shell-model
calculations were performed using angular momentum projected basis states with
the NuShellX shell-model code nushellx  framework. For the Fe isotopes, the
model space is comprised of the $0f_{7/2}$, $0f_{5/2}$, $1p_{3/2}$, and
$1p_{1/2}$ orbitals and the Hamiltonian is defined by the one- and two-body
matrix elements of the GXPF1A interaction of Ref. gxpf1a . The model space for
the Ge isotopes consists of the $0f_{5/2}$, $1p_{3/2}$, $1p_{1/2}$, $0g_{9/2}$
orbitals. For the Ge isotopes, we present results for two different empirical
Hamiltonians: 1) $jj44b$ defined in the appendix of Ref. Muk  and 2) jun45 of
Ref. Homna_2009 . Note that there are no spurious center-of-mass excitations
in either of these model spaces.

##  III Computing the Hamiltonian moments with Lanczos

At its core, the Lanczos algorithm is a really moment method; efficiently
computing $2n$ moments of $\hat{H}$ with respect to the initial pivot vector
$|v_{1}\rangle$ after $n$ iterations. With the choice of
$|v_{1}\rangle=\frac{1}{\sqrt{N}}\sum_{i}\phi_{i}|\psi_{i}\rangle$, where
$\phi_{i}$ is a random phase, we find it is possible to efficiently compute
several moments of the Hamiltonian with just a few Lanczos iterations. This is
illustrated by the first Lanczos matrix element $\alpha_{1}$ given by

$\alpha_{1}=\frac{1}{N}\sum_{i}H_{ii}+\sum_{i\neq
j}\frac{\phi_{i}\phi_{j}}{N}H_{ji}=H_{1}+\sum_{i\neq
j}\frac{\phi_{i}\phi_{j}}{N}H_{ji}.$ (2)

The remainder in Eq. (2) is generally small due to cancellations caused by the
random phases and a diminishing magnitude due to the large factor $N$ in the
denominator. Thus, for systems with large dimensions $\alpha_{1}\approx
H_{1}$. If needed, higher accuracy can be obtained by using different random
initial pivots and averaging. A small remainder in Eq. (2) then suggests a
strategy to compute even higher moments $\hat{H}$ via

$M_{k}=\frac{1}{N}{\rm{Tr}}[(\hat{H}-H_{1})^{k}]\approx\langle
v_{1}|(\hat{H}-\alpha_{1})^{k}|v_{1}\rangle.$ (3)

To compute the moments with Lanczos iterations, we note the recurrence
relation for the $n^{th}$ Lanczos vector

$|v_{n}\rangle=\frac{\hat{h}-\alpha_{n-1}+\alpha_{1}}{\beta_{n-1}}|v_{n-1}\rangle-\frac{\beta_{n-2}}{\beta_{n-1}}|v_{n-2}\rangle,$
(4)

with $\hat{h}=\hat{H}-\alpha_{1}$ and
$|v_{2}\rangle=\frac{\hat{h}}{\beta_{1}}|v_{1}\rangle$. In the case that the
remainder elements are small, we have the approximation $M_{k}\approx\langle
v_{1}|\hat{h}^{k}|v_{1}\rangle$, which can be extracted from the Lanczos
matrix elements through successive application of the recurrence relation,
collecting powers of $\hat{h}$, and back substituting for previous moments.
From the $n^{th}$ Lanczos iteration, which gives the Lanczos vectors up to
$v_{n+1}$, the moment $M_{n+1}$ can be obtained from the normalization
condition $\langle v_{n+1}|v_{n+1}\rangle=1$, while the moment $M_{n}$ can be
extracted from the orthogonality of the Lanczos vectors, i.e., $\langle
v_{n}|v_{n+1}\rangle=0$. For example, $M_{2}$ can be found from normalizing
$|v_{2}\rangle$

$\langle v_{2}|v_{2}\rangle=\frac{\langle
v_{1}|{\hat{h}}^{2}|v_{1}\rangle}{\beta_{1}^{2}}=\frac{M_{2}}{\beta_{1}^{2}}=1,$
(5)

leading to

$M_{2}=\beta_{1}^{2}.$ (6)

For $M_{3}$, we use the orthogonality condition

$\displaystyle\langle v_{2}|v_{3}\rangle$ $\displaystyle=\frac{\langle
v_{2}|\hat{h}-(\alpha_{2}-\alpha_{1})|v_{2}\rangle}{\beta_{2}}-\frac{\beta_{1}}{\beta_{2}}\langle
v_{2}|v_{1}\rangle,$ (7) $\displaystyle=\frac{\langle
v_{1}|\hat{h}[\hat{h}-(\alpha_{2}-\alpha_{1})]\hat{h}|v_{1}\rangle}{\beta_{2}\beta_{1}^{2}},$
(8)
$\displaystyle=\frac{M_{3}}{\beta_{2}\beta_{1}^{2}}-\frac{\alpha_{2}-\alpha_{1}}{\beta_{2}}=0,$
(9)

giving

$M_{3}=\beta_{1}^{2}(-\alpha_{1}+\alpha_{2}).$ (11)

Overall, while the derivations are tedious, they are straightforward using the
symbolic manipulation program Mathematica. The first eight moments in terms of
the matrix elements from the first four Lanczos iterations are given by

$\displaystyle H_{1}=$ $\displaystyle\alpha_{1}$ (12) $\displaystyle M_{2}=$
$\displaystyle\beta_{1}^{2}$ (13) $\displaystyle M_{3}=$
$\displaystyle\beta_{1}^{2}(-\alpha_{1}+\alpha_{2})$ (14) $\displaystyle
M_{4}=$
$\displaystyle\beta_{1}^{2}(\alpha_{1}^{2}-2\alpha_{1}\alpha_{2}+\alpha_{2}^{2}+\beta_{1}^{2}+\beta_{2}^{2})$
(15) $\displaystyle M_{5}=$
$\displaystyle\beta_{1}^{2}\Bigl{(}-\alpha_{1}\left(3\alpha_{2}^{2}+2\beta_{1}^{2}+3\beta_{2}^{2}\right)+\alpha_{3}\beta_{2}^{2}+2\alpha_{2}\left(\beta_{1}^{2}+\beta_{2}^{2}\right)-\alpha_{1}^{3}+3\alpha_{2}\alpha_{1}^{2}+\alpha_{2}^{3}\Bigr{)}$
(16) $\displaystyle M_{6}=$
$\displaystyle\beta_{1}^{2}\Bigl{(}3\alpha_{1}^{2}\left(2\alpha_{2}^{2}+\beta_{1}^{2}+2\beta_{2}^{2}\right)-2\alpha_{1}\left(\alpha_{2}\left(3\beta_{1}^{2}+4\beta_{2}^{2}\right)+2\alpha_{3}\beta_{2}^{2}+2\alpha_{2}^{3}\right)+$
$\displaystyle\hskip
17.07182pt\alpha_{3}^{2}\beta_{2}^{2}+2\alpha_{2}\alpha_{3}\beta_{2}^{2}+3\alpha_{2}^{2}\left(\beta_{1}^{2}+\beta_{2}^{2}\right)+\alpha_{1}^{4}-4\alpha_{2}\alpha_{1}^{3}+\alpha_{2}^{4}+\beta_{1}^{4}+\beta_{2}^{4}+2\beta_{1}^{2}\beta_{2}^{2}+\beta_{2}^{2}\beta_{3}^{2}\Bigr{)}$
(17) $\displaystyle M_{7}=$
$\displaystyle\beta_{1}^{2}\Bigl{(}-2\alpha_{1}^{3}\left(5\alpha_{2}^{2}+2\beta_{1}^{2}+5\beta_{2}^{2}\right)+2\alpha_{1}^{2}\left(2\alpha_{2}\left(3\beta_{1}^{2}+5\beta_{2}^{2}\right)+5\alpha_{3}\beta_{2}^{2}+5\alpha_{2}^{3}\right)-$
$\displaystyle\hskip
17.07182pt\alpha_{1}\Bigl{(}3\alpha_{2}^{2}\left(4\beta_{1}^{2}+5\beta_{2}^{2}\right)+10\alpha_{3}\alpha_{2}\beta_{2}^{2}+5\beta_{2}^{2}\left(\alpha_{3}^{2}+\beta_{2}^{2}+\beta_{3}^{2}\right)+5\alpha_{2}^{4}+3\beta_{1}^{4}+8\beta_{1}^{2}\beta_{2}^{2}\Bigr{)}+$
$\displaystyle\hskip
17.07182pt3\alpha_{2}^{2}\alpha_{3}\beta_{2}^{2}+4\alpha_{2}^{3}\left(\beta_{1}^{2}+\beta_{2}^{2}\right)+\beta_{2}^{2}\left(2\alpha_{3}\left(\beta_{1}^{2}+\beta_{2}^{2}+\beta_{3}^{2}\right)+\alpha_{4}\beta_{3}^{2}+\alpha_{3}^{3}\right)+$
$\displaystyle\hskip
17.07182pt\alpha_{2}\left(\beta_{2}^{2}\left(2\alpha_{3}^{2}+3\beta_{2}^{2}+2\beta_{3}^{2}\right)+3\beta_{1}^{4}+6\beta_{2}^{2}\beta_{1}^{2}\right)-\alpha_{1}^{5}+5\alpha_{2}\alpha_{1}^{4}+\alpha_{2}^{5}\Bigr{)}$
(18) $\displaystyle M_{8}=$
$\displaystyle\beta_{1}^{2}\Bigl{(}5\alpha_{1}^{4}\left(3\alpha_{2}^{2}+\beta_{1}^{2}+3\beta_{2}^{2}\right)-20\alpha_{1}^{3}\left(\alpha_{2}\left(\beta_{1}^{2}+2\beta_{2}^{2}\right)+\alpha_{3}\beta_{2}^{2}+\alpha_{2}^{3}\right)+$
$\displaystyle\hskip
17.07182pt\alpha_{1}^{2}\left(15\alpha_{2}^{2}\left(2\beta_{1}^{2}+3\beta_{2}^{2}\right)+30\alpha_{3}\alpha_{2}\beta_{2}^{2}+15\beta_{2}^{2}\left(\alpha_{3}^{2}+\beta_{2}^{2}+\beta_{3}^{2}\right)+15\alpha_{2}^{4}+6\beta_{1}^{4}+20\beta_{1}^{2}\beta_{2}^{2}\right)-$
$\displaystyle\hskip
17.07182pt2\alpha_{1}\Bigl{(}2\alpha_{2}^{3}\left(5\beta_{1}^{2}+6\beta_{2}^{2}\right)+9\alpha_{3}\alpha_{2}^{2}\beta_{2}^{2}+3\alpha_{2}\left(\beta_{2}^{2}\left(2\alpha_{3}^{2}+3\beta_{2}^{2}+2\beta_{3}^{2}\right)+2\beta_{1}^{4}+5\beta_{2}^{2}\beta_{1}^{2}\right)+$
$\displaystyle\hskip
39.83368pt\beta_{2}^{2}\left(\alpha_{3}\left(5\beta_{1}^{2}+6\left(\beta_{2}^{2}+\beta_{3}^{2}\right)\right)+3\alpha_{4}\beta_{3}^{2}+3\alpha_{3}^{3}\right)+3\alpha_{2}^{5}\Bigr{)}+$
$\displaystyle\hskip
17.07182pt3\alpha_{3}^{2}\beta_{2}^{4}+\alpha_{3}^{4}\beta_{2}^{2}+2\alpha_{3}^{2}\beta_{1}^{2}\beta_{2}^{2}+4\alpha_{2}^{3}\alpha_{3}\beta_{2}^{2}+3\alpha_{3}^{2}\beta_{2}^{2}\beta_{3}^{2}+\alpha_{4}^{2}\beta_{2}^{2}\beta_{3}^{2}+2\alpha_{3}\alpha_{4}\beta_{2}^{2}\beta_{3}^{2}+5\alpha_{2}^{4}\left(\beta_{1}^{2}+\beta_{2}^{2}\right)+$
$\displaystyle\hskip
17.07182pt3\alpha_{2}^{2}\left(\beta_{2}^{2}\left(\alpha_{3}^{2}+2\beta_{2}^{2}+\beta_{3}^{2}\right)+2\beta_{1}^{4}+4\beta_{2}^{2}\beta_{1}^{2}\right)+2\alpha_{2}\beta_{2}^{2}\left(\alpha_{3}\left(3\beta_{1}^{2}+3\beta_{2}^{2}+2\beta_{3}^{2}\right)+\alpha_{4}\beta_{3}^{2}+\alpha_{3}^{3}\right)+$
$\displaystyle\hskip
17.07182pt\alpha_{1}^{6}-6\alpha_{2}\alpha_{1}^{5}+\alpha_{2}^{6}+\beta_{1}^{6}+\beta_{2}^{6}+3\beta_{1}^{2}\beta_{2}^{4}+\beta_{2}^{2}\beta_{3}^{4}+3\beta_{1}^{4}\beta_{2}^{2}+2\beta_{2}^{4}\beta_{3}^{2}+2\beta_{1}^{2}\beta_{2}^{2}\beta_{3}^{2}+\beta_{2}^{2}\beta_{3}^{2}\beta_{4}^{2}\Bigr{)},$
(19)

In addition, the scaled moments $R_{k}=M_{k}/\sigma^{k}$ (with
$\sigma^{2}=M_{2}\approx\beta_{1}^{2}$) can easily be computed using these
formulae with the substitutions $\alpha_{i}\rightarrow\alpha_{i}/|\beta_{1}|$
and $\beta_{i}\rightarrow\beta_{i}/|\beta_{1}|$.

The validity of Eqs. (12)-(III) is shown in Table 1, where the moments
extracted from the first four Lanczos (L) iterations from a single random
pivot are compared with the exact (Ex) moments for several nuclei within the
$1p0f$-shell model space using the GXPF1A interaction gxpf1a . These systems
were chosen because they have large dimensions, $N\approx 2-4\times 10^{4}$,
but are still small enough to fully diagonalize. For $M_{3-8}$, we show the
scaled moments $R_{k}=M_{k}/\sigma^{k}$. Overall, good agreement is obtained
between the exact and Lanczos-inferred moments. Some differences exist, which
tend to be larger for the higher moments, and are due to an imperfect
cancellation in the remainder term that propagates further into the higher
moments. We find, however, that the remainders in $H_{1}$ and $M_{2}$ decrease
with increasing model space size. We find that these inferred moments are more
than sufficient to describe the averaged properties of the Hamiltonian matrix
and to model the average properties of the remaining Lanczos matrix elements.

In general, most systems within the $1p0f$ shell have been found to have
$R_{4}\approx 2.8$, $R_{6}\approx 12$, and $R_{8}\approx 65-75$. For the
purpose of comparison, note that for a Gaussian distribution, $R_{4}=3$,
$R_{6}=15$, and $R_{8}=105$.

Table 1: Comparison between exact (Ex) moments and those computed with the first four Lanczos (L) iterations for selected nuclei in the $1p0f$-shell model space using the GXPF1A interaction. $H_{1}$ is in units of MeV, $M_{2}$ is units of MeV2, while $R_{3-8}$ are dimensionless. |  |  47Cr |  47Cr |  48Cr |  48Cr |  72Kr |  73Kr
---|---|---|---|---|---|---|---
|  | $1/2^{-}$ | $3/2^{-}$ | 0+ | 12+ | 0+ | $1/2^{-}$
$H_{1}$ | Ex | -46.326 | -46.402 | -55.004 | -59.195 | -363.738 | -380.331
| L | -46.335 | -46.401 | -54.996 | -59.166 | -363.695 | -380.364
$M_{2}$ | Ex | 94.722 | 94.052 | 111.121 | 76.011 | 110.502 | 95.473
| L | 94.766 | 93.284 | 111.828 | 75.645 | 110.853 | 97.063
$R_{3}$ | Ex | -0.067 | -0.070 | -0.072 | -0.092 | 0.021 | 0.039
| L | -0.089 | -0.066 | -0.067 | -0.100 | 0.026 | 0.071
$R_{4}$ | Ex | 2.756 | 2.753 | 2.803 | 2.737 | 2.768 | 2.723
| L | 2.763 | 2.780 | 2.777 | 2.765 | 2.763 | 2.710
$R_{5}$ | Ex | -0.612 | -0.644 | -0.685 | -0.784 | 0.223 | 0.375
| L | -0.711 | -0.620 | -0.703 | -0.817 | 0.234 | 0.535
$R_{6}$ | Ex | 11.742 | 11.724 | 12.421 | 11.515 | 11.875 | 11.190
| L | 11.700 | 11.866 | 12.387 | 11.894 | 11.817 | 11.331
$R_{7}$ | Ex | -5.359 | -5.706 | -6.505 | -6.533 | 2.217 | 3.325
| L | -5.436 | -5.457 | -7.776 | -6.656 | 1.916 | 3.930
$R_{8}$ | Ex | 65.370 | 65.441 | 74.272 | 63.201 | 66.940 | 59.537
| L | 63.491 | 65.525 | 77.997 | 67.283 | 66.255 | 61.830

As mentioned above, higher accuracy can be achieved by computing the moments
stochastically; that is by using $N_{\rm samp}$ different initial pivots
$|v_{1}^{j}\rangle$ and averaging the resulting moments, i.e.,

$M_{k}\approx\frac{1}{N_{\rm samp}}\sum_{j}\langle
v_{1}^{j}|\hat{h}^{k}|v_{1}^{j}\rangle.$ (20)

The variance divided by the square root of the number of samples then provides
an estimate the error. This is shown in Figure 1 for the $J^{\pi}=0^{+}$ basis
in 48Cr for $N_{\rm samp}=10$ different initial random pivots (each sample is
indicated by the black dots connected with the black line and labeled on the
$x$-axis by the index $j$) for $H_{1}$, $\sigma$, and $R_{3-8}$ (labeled as
$R_{3-8}$ in the figure). The solid blue line represents the running average
for each moment, the dashed blue line shows the error in the averaging, and
the solid red line is the exact result. The figure shows that for this
relatively small system, any single initial pivot provides result with an
accuracy of a few percent.

In Figure 2, we show moments extracted for 10 different initial random pivots
for the $J^{\pi}=1/2^{-}$ states in 57Fe. Again, the individual results are
represented by the black points, while the solid and dashed blue lines
represent the running average and the estimated error, respectively. We note
that because of the large dimension of this system, $N=13436903$, the
variation in the individual samples is quite small; amounting to less than one
percent. The exact results for $H_{1}$ and $\sigma^{2}$, as computed with the
computer code of Ref. Horoi , are shown with the red lines. Each of the
initial pivots agree with $H_{1}$ to within 10 keV and $\sigma$ to within 5
keV, and the averaged moments are in excellent agreement with the exact
result. This demonstrates that the Lanczos procedure to compute the moments
improves with dimension.

Figure 1: (color online) Moments ($H_{1}$, $\sigma$, and $R_{3-8}$) computed
with 10 initial random pivots for the $J^{\pi}=0^{+}$ basis in 48Cr. The
results for each initial vector $v_{1}^{j}$ are indicated with the black dots
connected with the black line and labeled on the $x$-axis by $j$. The solid
blue line represents the running average for each moment, the dashed blue
shows the error in the averaging, and the solid red line is the exact result.
Figure 2: (color online) Moments ($H_{1}$, $\sigma$, and $R_{3-8}$) computed
with 10 initial random pivots for the $J^{\pi}=1/2^{-}$ basis in 57Fe. The
results for each initial vector $v_{1}^{j}$ are indicated with the black dots
connected with the black line and labeled on the $x$-axis by $j$. The solid
blue line represents the running average for each moment, the dashed blue
shows the error in the averaging, and the solid red line is the exact result
for $H_{1}$ and $\sigma^{2}$.

In Figures 3 and 4, the dependence on angular momentum [in particular, the
square $J(J+1)$] of the first eight moments is shown for calculations of both
57Fe and 74Ge. The 74Ge results were obtained with the $jj44b$ interaction of
Ref. Muk . A strong dependence on the square of the angular momentum is
demonstrated for both the first and second moments for both nuclei. For 57Fe,
the scaled higher moments $R_{k}$ exhibit a weak additional dependence on
angular momentum. On the other hand, in 74Ge, the higher scaled moments show a
marked decrease with increasing angular momentum. Indeed, the eigenspectrum
transitions to a more Gaussian-like distribution since $R_{8}$ decreases from
a large value of 150 to 100. Also, we note that for low angular momenta
$R_{4}>3$. Lastly, the moments for the positive- and negative-parity spaces
are nearly identical.

Figure 3: (color online) The 57Fe moments ($H_{1}$, $\sigma$, and $R_{3-8}$)
as a function of the square of the angular momentum $J(J+1)$.  Figure 4:
(color online) The 74Ge moments ($H_{1}$, $\sigma$, and $R_{3-8}$) as a
function of the square of the angular momentum $J(J+1)$. The black line shows
the dependence for positive-parity states ($J^{+}$), while the red line shows
the negative-parity states ($J^{-}$).

##  IV Modeling the Lanczos Matrix Elements

For large dimensions (e.g., $>10^{8}$), the computation effort for a shell-
model calculation is determined by the Lanczos method; in particular the
application of the Hamiltonian to the pivot vectors to generate the tri-
diagonal matrix. The resulting tri-diagonal matrix with dimensions of
$10^{1-3}$ can easily be diagonalized in a few seconds, while a tri-diagonal
matrix with a dimension of the order $10^{5}$ can be diagonalized within a few
minutes. Thus, our goal is to develop a method to model the entire tri-
digaonal matrix based on the first eight moments. We propose the polynomial
form defining the Lanczos matrix elements at each iteration $i$ as

$\displaystyle\alpha_{i}=$ $\displaystyle
a_{0}+a_{1}z_{i}+{a_{2}}z_{i}^{2}+{a_{3}}z_{i}^{3}$ (21)
$\displaystyle\beta^{2}_{i}=$ $\displaystyle
b_{1}z_{i}[1+{b_{2}}z_{i}+{b_{3}}z_{i}^{2}+{b_{4}}z_{i}^{3}],$ (22)

where $z_{i}={\rm ln}(i/N)$. We note that this representation is different
from the inverse binomial of Ref. zuker  and the shifted Gaussian of Ref.
ormand . This representation provides the flexibility to accurately model the
Lanczos matrix elements for a wide range of systems including those where the
scaled fourth moment is greater than the Gaussian limit, $R_{4}>3$, as is
encountered with Ge isotopes. In addition, the large $N$ limit leads to useful
analytic formulae for the moments that can be useful to fix the parameters.

The $a$\- and $b$-coefficients can determined by requiring that the moments of
the modeled matrix elements reproduce moments of the Hamiltonian. We note that
while the moments are in general high-order polynomials in the $a$\- and
$b$-parameters, they are, themselves, most sensitive to the odd and even
moments, respectively. Further, the dominant parameter is $b_{1}$, which
effectively determines the second moment $M_{2}$. Also, $a_{0}$ is trivially
constrained by $H_{1}$ since it does not affect any of the higher moments.
Lastly, we note that many systems (although not all as, is observed later for
76Ge) have nearly the same value for $b_{2}$. This is due to the fact, as seen
in Table 1, that $R_{4}\approx 2.7-2.8$, which is close the Gaussian limit of
3.

The first eight moments of the tri-diagonal matrix can be computed via

$\displaystyle H_{1}=$ $\displaystyle\langle\alpha\rangle$ (23) $\displaystyle
M_{2}=$
$\displaystyle\langle(\alpha-\langle\alpha\rangle)^{2}\rangle+2\langle\beta^{2}\rangle$
(24) $\displaystyle M_{3}\approx$
$\displaystyle\langle(\alpha-\langle\alpha\rangle)^{3}\rangle+6\langle(\alpha-\langle\alpha\rangle)\beta^{2}\rangle$
(25) $\displaystyle M_{4}\approx$
$\displaystyle\langle(\alpha-\langle\alpha\rangle)^{4}\rangle+12\langle(\alpha-\langle\alpha\rangle)^{2}\beta^{2}\rangle+6\langle\beta^{4}\rangle$
(26) $\displaystyle M_{5}\approx$
$\displaystyle\langle(\alpha-\langle\alpha\rangle)^{5}\rangle+20\langle(\alpha-\langle\alpha\rangle)^{3}\beta^{2}\rangle+$
$\displaystyle 30\langle(\alpha-\langle\alpha\rangle)\beta^{4}\rangle$ (27)
$\displaystyle M_{6}\approx$
$\displaystyle\langle(\alpha-\langle\alpha\rangle)^{6}\rangle+30\langle(\alpha-\langle\alpha\rangle)^{4}\beta^{2}\rangle+$
$\displaystyle
90\langle(\alpha-\langle\alpha\rangle)^{2}\beta^{4}\rangle+20\langle\beta^{6}\rangle$
(28) $\displaystyle M_{7}\approx$
$\displaystyle\langle(\alpha-\langle\alpha\rangle)^{7}\rangle+42\langle(\alpha-\langle\alpha\rangle)^{5}\beta^{2}\rangle+$
$\displaystyle
210\langle(\alpha-\langle\alpha\rangle)^{3}\beta^{4}\rangle+140\langle(\alpha-\langle\alpha\rangle)\beta^{6}\rangle$
(29) $\displaystyle M_{8}\approx$
$\displaystyle\langle(\alpha-\langle\alpha\rangle)^{8}\rangle+56\langle(\alpha-\langle\alpha\rangle)^{6}\beta^{2}\rangle+$
$\displaystyle
420\langle(\alpha-\langle\alpha\rangle)^{4}\beta^{4}\rangle+560\langle(\alpha-\langle\alpha\rangle)^{2}\beta^{6}\rangle+$
$\displaystyle 70\langle\beta^{8}\rangle,$ (30)

where $\langle...\rangle\rightarrow\frac{1}{N}\sum_{i}...$, which for large
$N$ can be extended to the integral $\frac{1}{N}\int_{1}^{N}...dx$. The
approximate equality arises from the assumption that adjacent matrix elements
$\beta_{i}$, $\beta_{i\pm 1}$, $\beta_{i\pm 2}$, $\beta_{i\pm 3}$ are nearly
equal. With Eqs. (23)-(IV) the $a$\- and $b$-parameters can be “fit” to
reproduce the moments of the Hamiltonian; leading to a modeled tri-diagonal
matrix with the same moments as the original Hamiltonian.

In principle, analytic formulae can be obtained for the moments in the large
$N$ limit since

$\lim_{N\rightarrow\infty}\int_{1}^{N}\ln^{m}xdx=m!.$ (31)

In this limit, the first five moments as defined in Eqs. (23)-(IV) are given
in terms of the $a$\- and $b$-parameters of Eqs. (21) and (22) by

$\displaystyle H_{1}=$ $\displaystyle a_{0}-a_{1}+2a_{2}-6a_{3},$ (32)
$\displaystyle M_{2}=$ $\displaystyle
a_{1}^{2}+\left(36a_{3}-8a_{2}\right)a_{1}+4\left(5a_{2}^{2}-54a_{3}a_{2}+171a_{3}^{2}\right)+2b_{1}\left(2b_{2}-6b_{3}+24b_{4}-1\right)$
(33) $\displaystyle M_{3}=$
$\displaystyle-2\Bigl{[}a_{1}^{3}-18\left(a_{2}-6a_{3}\right)a_{1}^{2}+12\left(10a_{2}^{2}-135a_{3}a_{2}+513a_{3}^{2}\right)a_{1}-$
$\displaystyle\hskip
22.76228pt8\left(37a_{2}^{3}-837a_{3}a_{2}^{2}+7047a_{3}^{2}a_{2}-21897a_{3}^{3}\right)\Bigr{]}+$
$\displaystyle
6b_{1}\Bigl{[}18a_{3}\left(-6b_{2}+38b_{3}-272b_{4}+1\right)+a_{1}\left(-4b_{2}+18b_{3}-96b_{4}+1\right)+4a_{2}\left(5b_{2}-27b_{3}+168b_{4}-1\right)\Bigr{]}$
(34) $\displaystyle M_{4}=$ $\displaystyle
3\Bigl{[}3a_{1}^{4}+\left(552a_{3}-80a_{2}\right)a_{1}^{3}+8\left(113a_{2}^{2}-1746a_{3}a_{2}+7515a_{3}^{2}\right)a_{1}^{2}-$
$\displaystyle\hskip
11.38092pt32\left(158a_{2}^{3}-4059a_{3}a_{2}^{2}+38412a_{3}^{2}a_{2}-132921a_{3}^{3}\right)a_{1}+$
$\displaystyle\hskip
11.38092pt16\left(731a_{2}^{4}-27540a_{3}a_{2}^{3}+427014a_{3}^{2}a_{2}^{2}-3208572a_{3}^{3}a_{2}+9800919a_{3}^{4}\right)\Bigr{]}+$
$\displaystyle
12b_{1}\Bigl{[}a_{1}^{2}\left(14b_{2}-78b_{3}+504b_{4}-3\right)-4a_{1}\Bigl{(}a_{2}\left(44b_{2}-282b_{3}+2064b_{4}-8\right)-9a_{3}\left(32b_{2}-234b_{3}+1928b_{4}-5\right)\Bigr{)}+$
$\displaystyle\hskip
25.6073pt4\Bigl{(}a_{2}^{2}\left(158b_{2}-1146b_{3}+9384b_{4}-25\right)-36a_{3}a_{2}\left(65b_{2}-531b_{3}+4844b_{4}-9\right)+$
$\displaystyle\hskip
25.6073pt9a_{3}^{2}\left(1082b_{2}-9846b_{3}+99144b_{4}-133\right)\Bigr{)}\Bigr{]}+$
$\displaystyle
12b_{1}^{2}\Big{[}12b_{2}^{2}-6\left(20b_{3}-120b_{4}+1\right)b_{2}+360b_{3}^{2}+20160b_{4}^{2}+b_{3}\left(24-5040b_{4}\right)-120b_{4}+1\Bigr{]}$
(35) $\displaystyle M_{5}=$
$\displaystyle-4\Big{[}11a_{1}^{5}+10\left(43a_{2}-342a_{3}\right)a_{1}^{4}-20\left(371a_{2}^{2}-6507a_{3}a_{2}+31410a_{3}^{2}\right)a_{1}^{3}+$
$\displaystyle\hskip
22.76228pt40\left(1756a_{2}^{3}-50625a_{3}a_{2}^{2}+532332a_{3}^{2}a_{2}-2029563a_{3}^{3}\right)a_{1}^{2}-$
$\displaystyle\hskip
22.76228pt80\left(4534a_{2}^{4}-189909a_{3}a_{2}^{3}+3245859a_{3}^{2}a_{2}^{2}-26685153a_{3}^{3}a_{2}+88602417a_{3}^{4}\right)a_{1}+$
$\displaystyle\hskip
22.76228pt32\left(25411a_{2}^{5}-1442205a_{3}a_{2}^{4}+35446860a_{3}^{2}a_{2}^{3}-469283490a_{3}^{3}a_{2}^{2}+3331562805a_{3}^{4}a_{2}-10104948693a_{3}^{5}\right)\Bigr{]}+$
$\displaystyle
20b_{1}\Bigl{[}a_{1}^{3}\left(-64b_{2}+426b_{3}-3216b_{4}+11\right)+$
$\displaystyle\hskip
25.6073pt6a_{1}^{2}\Bigl{(}2a_{2}\left(119b_{2}-891b_{3}+7488b_{4}-18\right)-9a_{3}\left(202b_{2}-1694b_{3}+15792b_{4}-27\right)\Bigr{)}-$
$\displaystyle\hskip
25.6073pt12a_{1}\Bigl{(}a_{2}^{2}\left(988b_{2}-8238b_{3}+76416b_{4}-133\right)-36a_{3}a_{2}\left(466b_{2}-4313b_{3}+44036b_{4}-56\right)+$
$\displaystyle\hskip
54.06006pt9a_{3}^{2}\left(8764b_{2}-89298b_{3}+996336b_{4}-949\right)\Bigr{)}+$
$\displaystyle\hskip
25.6073pt8\Bigl{(}a_{2}^{3}\left(4534b_{2}-41754b_{3}+424416b_{4}-548\right)-27a_{3}a_{2}^{2}\left(4714b_{2}-47818b_{3}+531392b_{4}-513\right)+$
$\displaystyle\hskip
36.98866pt54a_{3}^{2}a_{2}\left(24245b_{2}-268947b_{3}+3246768b_{4}-2395\right)+$
$\displaystyle\hskip
36.98866pt27a_{3}^{3}\left(-181498b_{2}+2187714b_{3}-28528896b_{4}+16391\right)\Bigr{)}\Bigr{]}-$
$\displaystyle
120b_{1}^{2}\Bigl{[}-a_{2}\left(168b_{2}^{2}-6\left(400b_{3}-3240b_{4}+9\right)b_{2}+9720b_{3}^{2}+887040b_{4}^{2}-2400b_{4}-336b_{3}\left(525b_{4}-1\right)+5\right)+$
$\displaystyle\hskip
28.45274pta_{1}\Bigl{(}24b_{2}^{2}-3\left(100b_{3}-720b_{4}+3\right)b_{2}+1080b_{3}^{2}+80640b_{4}^{2}-300b_{4}-24b_{3}\left(735b_{4}-2\right)+1\Bigr{)}+$
$\displaystyle\hskip
28.45274pt9a_{3}\Bigl{(}136b_{2}^{2}-2\left(1100b_{3}-9960b_{4}+19\right)b_{2}+9960b_{3}^{2}+1102080b_{4}^{2}-2200b_{4}-$
$\displaystyle\hskip
54.06006pt272b_{3}\left(735b_{4}-1\right)+3\Bigr{)}\Bigr{]}$ (36)

For $k>5$, these formulae are more complicated with extremely large
coefficients. Nonetheless, the analytic formulae for $M_{3}$ and $M_{5}$ are
useful for providing initial estimates for the parameters $a_{1}$ and $a_{2}$.
An alternative, that is somewhat more efficient for the higher moments ($k\geq
5$), and was used here to determine the parameters, is to evaluate the moment
integrals numerically using $z$ as the integration variable, which involves
integrals of the form

$\frac{1}{N}\int_{\ln(1/N)}^{0}e^{z}z^{m}dz.$ (37)

Sufficient accuracy can be achieved using Simpson’s rule with $10^{5}$ points.
For numerical stability, the integrals can be evaluated by scaling relative to
$M_{2}$ by taking $a_{i}\rightarrow a_{i}/\sqrt{-b_{1}}$ followed by setting
$b_{1}\rightarrow-1$.

The procedure used here to find the $a$\- and $b$\- parameters is discussed in
Appendix A.

The utility of the moment method to describe the nuclear Hamiltonian is
illustrated in Figure 5 where the modeled (colored lines) Lanczos matrix
elements are compared with those obtained from a shell-model calculation
(black lines) for the 48Cr, $J^{\pi}=0^{+}$ (top) and 57Fe, $J^{\pi}=25/2^{-}$
(bottom) systems. For 48Cr the entire Lanczos matrix ($N=41355$) is plotted,
while for 57Fe, $J^{\pi}=25/2^{-}$ ($N=13752093$), 3074 Lanczos iterations
were performed and 100000 modeled matrix elements are shown. The 48Cr system
is somewhat typical where the dominant behavior observed in the Lanczos matrix
can be extracted from just the first four moments, i.e., $M_{3}$ to constrain
$a_{1}$ and $M_{2}$ and $M_{4}$ to constrain $b_{2}$ and $b_{4}$. Still, the
figure shows that using moments up to $M_{8}$ can improve the overall
description of modeled Lanczos matrix. The 57Fe system is different in that
the higher moments are essential. The figure shows that limiting to $M_{3}$ to
constrain $a_{1}$ is clearly inadequate and improvement is achieved only by
including the higher odd moments, and the best overall results are obtained
using all eight moments. The 57Fe case is also interesting as it has a
negative skewness ($M_{3}$), which is correctly captured with the Lanczos
method to compute the moments, but also seemingly contradicts the positive
values of ($\alpha_{i}-H_{1}$) shown for the first few thousand iterations.
Indeed, the diagonal matrix elements show a strong curvature and eventually
turn negative for large iteration number. This is captured in the higher odd
moments leading to quadratic and cubic terms in the modeled $\alpha_{i}$
matrix elements. Lastly, the $\beta_{i}$ at low iteration number are also
influenced by the higher even moment $M_{6}$.

Figure 5: (color online) Comparison between shell model (black) and modeled
Lanczos matrix elements $\alpha$ and $\beta$ for 48Cr, $J^{\pi}=0^{+}$ (top)
and 57Fe, $J^{\pi}=25/2^{-}$ (bottom) within the $1p0f$-model space using the
GXPF1A interaction gxpf1a . The colored curves show modeled Lanczos matrix
elements using Eqs. (21) and (22) with the indicated moments to constrain the
$a$\- and $b$-parameters.

##  V Estimating the Level Density

The density of states is a key nuclear property that has a significant impact
on reaction rates for statistical processes, such as radiative neutron
capture. For the most part, reaction models, such as Hauser-Feshbach Hauser-
Feshbach , have relied on a parameterization of the level density based on a
modified back-shifted Fermi gas approach such as was introduced by Gilbert and
Cameron Gilbert-Cameron . This approach requires knowledge about several
parameters such as the single-particle level-density parameter $a$, which may
depend on excitation energy, the pairing gap $\Delta$, and the spin cutoff
parameter. In addition, the back-shifted Fermi gas density is matched to the
low-lying spectrum where the level density is assumed to follow an exponential
form. The matching is accomplished by requiring that the exponential component
reproduces the cumulative density up to an excitation where the discrete
levels are both known and complete and requiring continuity in the logarithmic
derivate of the level density (equivalent to the inverse temperature) at the
matching energy. A drawback of this procedure is that the level-density
parameters are generally constrained by experimental knowledge, such as the
spacings of $l=0$ ($D_{0}$) and $l=1$ ($D_{1}$) resonances at the neutron
separation energy, $S_{n}$. These quantity are generally known only in systems
based on a stable target. For radiative neutron capture, the level density is
needed essentially up to the neutron separation energy.

One approach to generalize our knowledge of the level density is to use
theoretical structure models based on the microscopic physics involved, such
as the nuclear shell model, where high-quality empirical nuclear Hamiltonians
have been developed that are well-known to reproduce the low-lying spectra of
nuclei. It is important to note that these shell-model calculations are based
on a finite model space, and at some excitation energy, $E_{x}$, they will
fail to adequately enumerate the system due to the presence of so-called
“intruder” states. These intruder states, however, are expected to occur at
higher excitation energies, generally of the order of the shell gap for states
with opposite parity and twice the shell gap for states of the same parity.
Thus, in many cases it is not unreasonable to hope that a large-basis shell-
model calculation contains contains sufficient configurations to adequately
describe the states of a given parity up to excitation energies near the
neutron separation energy. This supposition can be tested in a few cases
through comparison with experimentally measured resonance spacings. For
example, within the $1p0f$-shell, the calculated density of states can be
compared with the $l=1$ spacings $D_{1}$, at which point, the computed level
density can be used to define parameters of the back-shifted Fermi gas needed
to describe the full level density.

The most straight forward approach to compute the density of states within the
shell model would be to simply diagonalize the model Hamiltonian and count the
respective states. In many cases, this is computationally prohibitive since
the number of the configurations within the model space can exceed $10^{9}$.
Instead, since the density of states is more of a statistical property of the
Hamiltonian, we propose to model the Hamiltonian via the moments method
outlined above and to compute the density of states from the modeled matrix.
Another approach would be to use the binomial distribution described in Ref.
zuker , which is constrained with just the first four moments of the
Hamiltonian and is appealing due to its analytic nature. In what follows,
several approaches to determine the density of states as a function of
excitation energy are outlined.

###  V.1 Extrapolated Lanczos Method

Section IV illustrated that for most cases the global, or averaged, properties
of the Lanczos matrix can be predicted from just four Lanczos iterations. This
offers a strategy to predict the statistical properties of the entire energy
spectrum by performing a set of Lanczos iterations sufficient to describe the
low-lying spectrum and then extrapolate the Lanczos matrix elements with Eqs.
(21) and (22) to an iteration number sufficient to properly estimate the
density of states. We refer to this as ELM($k$,$N_{\rm Lanc}$), where $k$
denotes the maximum moment $M_{k}$ used to extrapolate the Lanczos matrix
elements and $N_{\rm Lanc}$ is the number of actual Lanczos iterations used
prior to extrapolation. In general, the Lanczos iterations can be
computationally expensive for large model spaces, and a key question is just
what value of $N_{\rm Lanc}$ is sufficient and/or optimal. A general
requirement is obtaining sufficient accuracy in the ground-state energy,
$E_{gs}$, to establish the excitation energy scale to measure the density of
states. The accuracy required in $E_{gs}$ is model space and Hamiltonian
dependent. For example, for the model spaces and Hamiltonians studied in this
work, we found that an uncertainty of 10 keV in $E_{gs}$ leads to a 1%
uncertainty in the level density, while a 100 keV uncertainty leads to a 10%
change in the level density. As a general rule, 30 - 40 Lanczos iterations are
needed to determine the ground-state energy with an accuracy better than 10
keV, and more often than not, with an accuracy of 1 keV. To some degree, an
optimal number of Lanczos iterations can be thought of as where a smooth
transition (within the fluctuations of the Lanczos matrix elements) occurs
between the computed and modeled Lanczos matrix elements. This may not always
be practical, and while it is true that too few iterations can lead to
difficulties in the direct computation of the level density at lower
excitation energies, an analytic continuation method, discussed below, can
address this issue. Consequently, it is often possible to achieve excellent
results with the ELM method with $N_{\rm Lanc}$ as low as 40.

In Figure 6, results for the $J^{\pi}=0^{+}$ space in 48Cr are shown. The
shell model calculation was performed using the GXPF1A interaction gxpf1a
within the $1p0f$-shell model space with the shell model-code NuShell. Here,
the full shell-model matrix was diagonalized with the Lanczos algorithm. The
black lines show the results from the shell-model calculation with the Lanczos
matrix elements displayed in the top half of the figure and the level density
and cumulative density shown in the left and right sides, respectively, in the
bottom half of the figure. The level density was computed as function of
excitation energy in steps of 100 keV as a running average within an energy
window of $E_{x}\pm 500$ keV, which smooths out fluctuations in the level
density. The red and blue lines show the results for ELM(8,40) and ELM(8,100),
respectively, where the Lanczos matrix was extrapolated to 50,000 iterations.
The ELM(8,100) calculation is nearly indistinguishable from the shell model
calculation. The ELM(8,40) calculation shows a slight deviation from the exact
shell-model calculation at $E_{x}\approx 6$ MeV. This deviation is primarily
due to a small discontinuity in the matching of the Lanczos matrix elements at
$N_{\rm Lanc}$ and hints at how the ELM($k$,$N_{\rm Lanc}$) approach can break
down.

Figure 6:  (color online) Results for the $J^{\pi}=0^{+}$ space in 48Cr within
the $1p0f$-shell model space with the GXPF1A interaction. The black lines show
the shell-model calculations for the Lanczos matrix elements in the upper half
of the figure and the level density and cumulative level density in the
bottom. In the lower half of the figure the level density and cumulative
density are shown for ELM(8,40) (red) and ELM(8,100) (blue).

In addition to the demonstration for 48Cr, we have also applied and tested the
ELM method to 57Fe for $J^{\pi}=1/2^{-}-25/2^{-}$ and 76Ge for
$J=0^{\pm}-14^{\pm}$. In what follows, representative results for these
systems are shown to demonstrate various features of the ELM method. We note
that applications of the ELM(2,100) method to the Fe region were published
earlier in Ref. CERN .

Shown in Figures 7 and 8 are the results obtained for the $1/2^{-}$ and
$25/2^{-}$ states in 57Fe, while the moments are given in Table 2. Again, the
solid black lines are the results fo the shell-model calculation, while the
red and blue lines represent the ELM(8,40) and ELM(8,100) results,
respectively. The level densities were computed by extrapolating the the
Lanczos matrix elements to 150,000 iterations, diagonalizing the resulting
matrix, and as a running average over an excitation energy window of $E_{x}\pm
500$ keV. The primary difference between the $1/2^{-}$ and $25/2^{-}$ angular
momentum spaces lies with the odd moments. Both systems have nearly identical
negative skewness ($R_{3}$) as is shown in Table 2. The high-spin state,
however, has a large non-linear term, and the ($\alpha_{i}-H_{1}$) are
actually positive for smaller iteration number, and then decrease and become
negative at large iteration number. A signature of this behavior is also
exhibited in the higher odd moments. In particular, when $M_{3}$ dominates the
spectral behavior (linear terms in the $\alpha_{i}$), one often finds
$R_{7}\sim 9.0-9.5R_{5}$ and $R_{5}\sim 9.0-9.5R_{3}$. Instead, for the
$25/2^{-}$ space $R_{7}\sim 7.3R_{5}$ and $R_{5}\sim 8R_{3}$.

Figure 7:  (color online) Results for the $J^{\pi}=1/2^{-}$ space in 57Fe within the $1p0f$ shell-model space with the GXPF1A interaction. The black lines show the shell-model calculations for the Lanczos matrix elements in the upper half of the figure and the level density and cumulative level density in the bottom. In the lower half of the figure the level density and cumulative density are shown for ELM(8,40) (red) and ELM(8,100) (blue).  Figure 8:  (color online) Results for the $J^{\pi}=25/2^{-}$ space in 57Fe within the $1p0f$ shell-model space with the GXPF1A interaction. The black lines show the shell-model calculations for the Lanczos matrix elements in the upper half of the figure and the level density and cumulative level density in the bottom. In the lower half of the figure the level density and cumulative density are shown for ELM(8,40) (red) and ELM(8,100) (blue).  Table 2: Comparison of moments computed with the first four Lanczos iterations for $1/2^{-}$ and $1/2^{-}$ angular momentum configuration space in 57Fe. $H_{1}$ is in units of MeV, $M_{2}$ is units of MeV2, while $R_{3-8}$ are dimensionless. |  57Fe |  57Fe
---|---|---
| $1/2^{-}$ | $25/2^{-}$
$H_{1}$ | -143.314 | -145.213
$M_{2}$ | 179.268 | 140.764
$R_{3}$ | -0.026 | -0.022
$R_{4}$ | 2.839 | 2.828
$R_{5}$ | -0.244 | -0.176
$R_{6}$ | 12.726 | 12.595
$R_{7}$ | -2.229 | -1.287
$R_{8}$ | 75.703 | 74.324

This section demonstrating the ELM approach is concluded with an examination
of the $J^{\pi}=0^{+}$ and $4^{+}$ systems in 76Ge using the $jj44b$
interaction of Ref. Muk . The computed moments are shown in Table 3. The key
features of this system are: 1) the large skewness ($R_{3}\sim 0.2$), which is
an order of magnitude larger than that observed in 57Fe, 2) the large fourth
moment ($R_{4}>3$, which is substantially larger than the Gaussian value of
3), and 3) the dramatic difference in the $8^{th}$ moment between the two
angular momenta.

Table 3: Comparison of moments computed with the first four Lanczos iterations for $0^{+}$ and $4^{+}$ states in 76Ge. $H_{1}$ is in units of MeV, $M_{2}$ is units of MeV2, while $R_{3-8}$ are dimensionless. |  76Ge |  76Ge
---|---|---
| $0^{+}$ | $4^{+}$
$H_{1}$ | -190.500 | -190.544
$M_{2}$ | 47.911 | 46.021
$R_{3}$ | 0.228 | 0.201
$R_{4}$ | 3.266 | 3.135
$R_{5}$ | 3.079 | 2.441
$R_{6}$ | 22.298 | 18.436
$R_{7}$ | 53.417 | 32.914
$R_{8}$ | 310.668 | 180.656

Shown in Figures 9 and 10 are the results for $0^{+}$ and $4^{+}$ states,
respectively, for 76Ge obtained with the $jj44b$ interaction. The level
density was computed by extrapolating the Lanczos matrix to a dimension of
150,000 and computing a running average within the excitation energy of
$E_{x}\pm 500$ keV. For illustrative purposes, approximately 1000 lanczos
iterations were performed in each space to to diagnose the calculation in the
level density. The results for the $J^{\pi}=0^{+}$ space are similar to those
shown earlier for 48Cr and 57Fe where the ELM(8,100) closely matches the
shell-model result. This is not the case, however, for the $J^{\pi}=4^{+}$
where there is a clear discrepancy in the spectrum at $E_{x}\approx 3-5$ MeV.
On the other hand, the ELM(8,$N_{\rm Lanc}$) results agree with the shell
model at higher excitation energies, as would be expected since this is the
regime where the statistical nature of the configuration space should dominate
the spectral behavior. The cause of this discrepancy is evident in the upper
part of the figure where the diagonal $\alpha_{i}$ matrix elements exhibit a
clear transition in their behavior. The figure shows that the modeled matrix
elements capture the overall behavior of the Lanczos matrix elements for large
iteration number, but fail to describe the “step” behavior shown to at
approximately 400 iterations. Thus, the modeled matrix elements lead to a
strong dip in the level density for $E_{x}\approx 3-5$ MeV that is caused by a
strong discontinuity between the modeled and actual matrix elements that is
far larger than scatter, or noise, exhibited in the computed Lanczos matrix
elements. In this case, it would be necessary to perform an ELM(8,400)
calculation in order to more accurately describe the system. It has to be
noted that often times such a calculation can be computationally prohibitive.
In addition, while these calculations for 76Ge are quite different than those
in the $1p0f$ shell, it is not always clear if, or where, a sudden transition
in the computed matrix elements may take place; especially for model spaces
involving orbits in different major shells. As is apparent from the upper part
of Figure 10, the clearest signature of a potential problem with the ELM
procedure is the existence of a strong discontinuity at $N_{\rm Lanc}$ between
the compute Lanczos matrix elements and the modeled matrix elements. This
discontinuity may be present in either the $\alpha_{i}$ matrix elements, the
$\beta_{i}$ matrix elements, or both. If such a discontinuity exists, two
alternatives are suggested: 1) an alternative extrapolation between the
computed and modeled matrix elements that smoothly joins the matrix elements
to within the “noise” in the matrix elements, or 2) a procedure to
analytically continue the level density from the high-energy regime to the
lowest state in the model space. The latter approach will be discussed in
Section V.3.

Figure 9:  (color online) Results for the $J^{\pi}=0^{+}$ space in 76Ge within
the $jj44$ shell-model space with the $jj44b$ interaction. The black lines
show the shell-model calculations for the Lanczos matrix elements in the upper
half of the figure and the level density and cumulative level density in the
bottom. In the lower half of the figure the level density and cumulative
density are shown for ELM(8,40) (red) and ELM(8,100) (blue).  Figure 10:
(color online) Results for the $J^{\pi}=4^{+}$ space in 76Ge within the $jj44$
shell-model space with the $jj44b$ interaction. The black lines show the
shell-model calculations for the Lanczos matrix elements in the upper half of
the figure and the level density and cumulative level density in the bottom.
In the lower half of the figure the level density and cumulative density are
shown for ELM(8,40) (red) and ELM(8,100) (blue).

###  V.2 Binomial Approximation for the Level Density

In Ref. zuker_2 , a binomial form was proposed to describe the density of
states for quantum many-body systems, such as those described by the nuclear
shell model. For a system of dimension $N$, three parameters are required to
define the binomial: $\cal N$ the effective dimension of the system, the
asymmetry $p$, and an energy scale $\epsilon$. The span $S$ (the energy
difference between the lowest and highest states), centroid $E_{c}$, variance
$\sigma^{2}$, and dimensionless energy $x$ are given by

$S={\cal N}\epsilon,\hskip 4.26773ptE_{c}={\cal N}p\epsilon,\hskip
4.26773pt\sigma^{2}={\cal N}pq\epsilon^{2},\hskip 4.26773ptx=\frac{E}{S},$
(38)

where $p+q=1$ and obviously $E_{c}=H_{1}$ and $\sigma^{2}=M_{2}$. The binomial
approximation to the level density is then given by

$\rho_{b}(x)=p^{x{\cal N}}q^{\bar{x}{\cal N}}\frac{\Gamma({\cal
N}+1)}{\Gamma(x{\cal N}+1)\Gamma(\bar{x}{\cal N}+1)}\frac{N{\cal N}}{S},$ (39)

with $\bar{x}=1-x$. The binomial parameters $p$ and ${\cal N}$ can be
determined by the 3rd and 4th moments of the Hamiltonian since for the
binomial

$R_{3}=\frac{q-p}{\sqrt{{\cal N}pq}}$ (40)

and

$R_{4}=3+\frac{1-6pq}{{\cal N}pq}.$ (41)

Defining $R=R_{3}^{2}/(R_{4}-3)$, the parameter $p$ becomes

$p=\frac{1}{2}\left[1-{\rm
sgn}(M_{3})\sqrt{1-2\left(\frac{1-R}{2-3R}\right)}~{}\right],$ (42)

from which, $\cal N$ follows directly from Eq. (41). With $p$ and $\cal N$
known, the span is then given by

$S=\sqrt{\frac{{\cal N}\sigma^{2}}{pq}}.$ (43)

In addition, for the binomial, the ground-state energy is $E_{gs}^{b}=-Sp$,
which may not correspond to the actual ground state energy $E_{gs}$. In this
case, the level density in Eq. (39) is shifted by $x-(E_{c}-Sp)/S$ so that the
binomial centroid corresponds to the centroid of the Hamiltonian relative to
the exact ground state. For the most part, the most significant hurdle in
implementing this approach has been the ability to compute $R_{3}$ and
$R_{4}$, which can now be computed using the Lanczos method.

Note from Eq. (42), a real solution with $0\leq p\leq 1$ requires $R\leq 0$,
which implies $R_{4}<3$ and is representative of systems approaching an
asymmetric Gaussian. Note that mathematically a solution for $p$ also exists
when $R>1$, which would imply $R_{4}>3$ with a very large asymmetry. This
solution, however, does not yield a physical solution where the $R_{3}$ and
$R_{4}$ moments of the binomial correspond to the actual moments. Thus, the
binomial is not applicable to the 76Ge results shown in Section V.1.

In Figure 11, results for the level density and cumulative density for the
$J^{\pi}=1/2^{-}$ and $25/2^{-}$ states in 57Fe are shown for the the binomial
approximation (green lines) and are compared to the ELM(8,100) (blue lines)
and the shell model (black lines) obtained with a finite number of Lanczos
iterations as specified in the figures. The figures show that both ELM and the
binomial approximation are in agreement at higher excitation energies where
the density of states is quite high. At lower excitation energies, the
binomial approximation can be poor since it lacks information about the ground
state and the low-lying spectrum, and in the case for the $J^{\pi}=1/2^{-}$
state in 57Fe, the “effective” lowest energy lies above the shell-model state.
This is not surprising since the binomial is limited to only four moments, and
as was already pointed out, one would need of the order 40 moments (20 Lanczos
iterations) for a reasonable calculation of the ground-state energy. In
addition, the low-energy behavior of the binomial is Gaussian-like, and thus,
the level density tends to decrease dramatically at low energy, giving an
effective lowest state so that $E_{0}^{\rm eff}>E_{0}$. For the most part, the
ELM procedure can provide a better description of the low-lying spectrum if
sufficient Lanczos iterations are performed in order to determine the energy
of the lowest state in the specified model space.

Figure 11:  (color online) Results for the level density and cumulative
density for $J^{\pi}=1/2^{-}$ and $25/2^{-}$states in 57Fe within the $1p0f$
shell-model space with the GXPF1A interaction. The black lines show the shell-
model calculation, while the blue and green lines represent ELM(8,100) and the
binomial approximation, respectively.

###  V.3 Analytic Continuation of the Level Density

As is shown in the previous sections, the level density modeled from the
moments and shifted relative to the exact ground-state energy is a good
representation of the exact shell-model level density at higher excitation
energies. The principal question, however, is how to properly describe the
level density in the situations illustrated in Figure 10 where the ELM has a
discontinuity in the Lanczos matrix elements and Figure 11 where the binomial
approximation substantially undershoots the shell-model result. In both cases,
the moments by themselves dramatically miss the lowest energy, $E_{0}$ in the
configuration space, leading to an “effective” $E_{0}^{\rm eff}$ that is too
high in energy. In principle, the ELM(8,$N_{\rm Lanc}$) procedure will work by
ensuring that the modeled and exact Lanczos matrix elements are reasonably
matched so that there isn’t a discontinuity larger than natural noise in the
calculated matrix elements. In some cases, however, the number of Lanczos
iterations, $N_{\rm Lanc}$ required would be prohibitively large, which in
effect negates any advantages in the approach.

A strategy for the case when $E_{0}^{\rm eff}>E_{0}$ is similar to that
outlined by Gilbert and Cameron Gilbert-Cameron  where the goal was to
describe the level density via two components: an exponentially increasing
function at low energy that is then matched to the back-shifted Fermi gas at
higher energies. Here, we take a similar approach by matching an exponentially
increasing level density to the ELM level density at a matching energy
$E_{m}$. Thus, at low energy, the density of states is taken to be

$\rho(E_{x})=\exp\left[(E_{x}-E_{\rm shift})/T\right].$ (44)

Note that $E_{\rm shift}$ specifies that the for cumulative density we have
$N(E_{\rm shift})=1$. The exponential level density of Eq. (44) can then be
matched at energy $E_{m}$ to the ELM or binomial approximation by requiring
continuity in the level density and by defining the temperature $T$ as the
inverse of the logarithmic derivative of $\rho$, i.e.,

$T(E_{x})=\frac{\rho(E_{x})}{\rho^{\prime}(E_{x})}.$ (45)

At a given $E_{m}$, the continuity requirement for the level density specifies
$E_{\rm shift}$ as

$E_{\rm shift}(E_{m})=E_{m}-T(E_{m})\ln\left[\rho(E_{m})\right].$ (46)

Thus, the matching energy can be chosen so that $E_{\rm shift}(E_{m})=E_{0}$.
Practical considerations for finding the matching energy for the ELM procedure
are given in Appendix B

Figure 12:  (color online) Results for the level density and cumulative
density for $J^{\pi}=1/2^{-}$ and $25/2^{-}$ states in 57Fe within the $1p0f$
shell-model space with the GXPF1A interaction. The black lines show the shell-
model calculation, while the red and green lines represent the ELMAC(8,40) and
binomial calculations, as described in the text, respectively. Figure 13:
(color online) Results for the level density and cumulative density for
$J^{\pi}=0^{+}$ and $4^{+}$ states in 76Ge within the $jj44$ shell-model space
with the $jj44b$ interaction. The black lines show the shell-model
calculation, while the red and blue lines represent ELM(8,40) and the
ELMAC(8,40) as described in the text, respectively.

In Figure 12, results for the ELM analytic continuation, ELMAC(8,40), and the
binomial level densities are shown for the $J^{\pi}=1/2^{-}$ and $25/2^{-}$
states in 57Fe, while the ELMAC(8,40) level density for the $J^{\pi}=0^{+}$
and $4^{+}$ states in 76Ge are shown in Figure 13 (note that the binomial
approach is not applicable due to $R_{4}>3$). Overall, the extrapolation works
well; especially when the effective lowest state for the modeled level density
is higher than the actual, i.e., $E_{0}^{\rm eff}>E_{0}$. Under this
condition, it is possible to smoothly match the modeled level density down to
the lowest state. As can be seen in Figures 12 and 13, however, in some cases,
such as the lower spin, the extrapolated level density tends to miss a “gap”
in the excitation spectrum at low excitation energies. This most likely
reflects the effect of pairing.

The case where $E_{0}^{\rm eff}<E_{0}$ is less common and is generally not
possible with the binomial level density due to the high curvature of the
Gaussian, which tends to decrease the level density dramatically at low
excitation energy. However, this can occur for the ELM when a small number of
actual Lanczos iterations, $N_{\rm Lanc}$, is used. On the other hand, for
ELM, better agreement with the low-lying spectrum is achieved with increasing
$N_{\rm Lanc}$, which is also needed in order to obtain a reasonable estimate
of $E_{0}$. Shown in Figure 14 is the case of the $J^{\pi}=15/2^{-}$ space in
57Fe (within the $1p0f$-shell model space with the GXPF1A interaction), where
results for ELM(8,4) (red), ELM(8,40) (green), and ELM(8,100) (blue) are shown
in comparison to the shell model with 219 Lanczos iterations (black) and the
analytically continued binomial approximation (orange). All the modeled level
densities are in agreement at high energy, where the spectrum is dominated by
the statistical properties of the Hamiltonian. The agreement between
ELM(8,$N_{\rm Lanc}$) and the shell model at low excitation energy improves
with increasing $N_{\rm Lanc}$ as is to be expected. Indeed, reasonable
agreement is achieved with ELM(8,40) which is close to the minimum number of
iterations needed to give an accurate energy for $E_{0}$ and the next level.

Figure 14:  (color online) Comparison of results with ELM(8,$N_{\rm Lanc}$ for
$N_{\rm Lanc}=4$ (red), $40$ (green) , and $100$ (blue) for the
$J^{\pi}=15/2^{-}$ space in 57Fe within the $1p0f$-shell model space with the
GXPF1A interaction. In addition the shell model with 219 Lanczos iterations
(black) and the analytically continued binomial (orange) are also shown.

Shown in Figure 15 are results for the various approaches for the summed over
angular momenta with fixed parity. The black lines show the results from the
Lanczos iterations while the red lines are the ELM(8,100) results. The blue
lines show the ELMAC(8,40) results. The green line is the result for the
analytically continued binomial, while the dashed green line is the binomial
(57Fe only). For the most part, both the ELM and binomial agree at high
excitation energy. In general, the most successful approach is ELM(8,$N_{\rm
Lanc}$) where $N_{\rm Lanc}$ is large enough to capture key features of the
Lanczos matrix elements. As discussed previously, this is when the difference
between the modeled and actual Lanczos matrix elements is less than the
natural “noise” in the matrix elements; that is no strong discontinuities. For
57Fe, this is generally achieved with $N_{\rm Lanc}\approx 50-100$. In this
sense, the ELM(8,100) results are likely representative of the full shell
model with 200,000 iterations. For 57Fe, analytically continuing the binomial
to $E_{0}$ is a significant improvement over the binomial itself. It does,
however, tend to underestimate the actual level density in the region
$E_{x}\approx 3-8$ MeV. For 76Ge, one would need $N_{\rm Lanc}\geq 1000$ in
order to avoid the most significant discontinuities. A situation that is less
than ideal. On the other hand, analytically continuing the ELM(8,40) gives a
good overall description of the level density.

Figure 15:  (color online) Comparison of results of the level density summing
all angular momenta of a given parity for 57Fe ($\sum_{J^{-}}$), 76Ge
($\sum_{J^{+}}$), and 76Ge ($\sum_{J^{-}}$). The black lines are from the
Lanczos iterations, the red line is the ELM(8,100) reconstruction, the blue
and green lines are ELMAC(8,40) results. The dashed green line is from the
binomial, while the green line is the analytic continuation of the binomial.
(57Fe only).

##  VI Applications of the ELM: 57Fe, 74Ge, and 76Ge

Figure 16:  (color online) Level densities for 57Fe. The black lines show the
calculated angular-momentum summed level density for negative-parity states up
to the value of $2J^{\pi}_{\rm max}$ as indicated to the right of each line.
The experimental $\ell=0$ and $1$ level densities mug  at $E_{x}=S_{n}$ are
shown by the red cross and circle, respectively, with error bars about the
size of the symbols. Note that the calculated level densities are only for the
negative parity states contained within the $1p0f$ shell model space. Other
data is for the sum of both negative and positive parity states. The red line
shows the experimental level density obtained from the states listed in NNDC
nndc . Level densities inferred from reaction data are shown by the shaded
areas: (green) 55Mn(3He,$\alpha$) reaction Voinov , (blue) 57Fe(3He,3He′)
reaction Algin , and (orange) 57Ni(p,p′) reaction Lar17 .

We now apply the extrapolated Lanczos method to compute the level density for
57Fe within the $0p1f$ model space using the GXPF1A interaction. The level
densities for each negative parity, angular momentum configuration space were
computed with ELM(8,100) and are shown in Figure 16. The black lines show the
angular-momentum summed level density for negative parity states up to the
$2J^{\pi}_{\rm max}$ value indicated to the right of each line. The
experimental $\ell$=1 level density mug  ($\rho_{1/2^{-}}$ \+
$\rho_{3/2^{-}}$) at $E_{x}=S_{n}$ (the neutron decay threshold) is shown with
the red circle (the error bar is approximately equal to the size of the
circle). The experimental value for $\ell$=0 level density mug
($\rho_{1/2^{-}}$) at $E_{x}=S_{n}$ is shown with the red cross (the error bar
is approximately equal to the size of the cross). Other data shown in the
figure is for the sum of both positive and negative parity states. The red
line shows the level density obtained from the experimentally observed states
listed in the NNDC nndc . The shaded areas are the bounds inferred from the
various reaction data (see Figure caption) Lar17 ; Algin ; Voinov . We note
that the $1p0f$ shell model space does not contain any positive parity states
for 57Fe.

The agreement between our calculation (the sum of densities for states with
1/2- and 3/2-) and the $\ell$=1 level density is excellent. In addition, the
level density for $1/2^{+}$ states is nearly the same as that computed for
$1/2^{-}$ states, which indicates that the parity ratio is close to unity at
$E_{x}=S_{n}$. Thus, our estimate of the total level density would be a factor
of two larger than shown in Fig. 16. The level density obtained from NNDC nndc
levels (see Figure 16) becomes about a factor of two larger than that
calculated for negative parity states starting around $E_{x}$ = 3 MeV,
indicating a parity ratio close to unity around 3 MeV. Taking this into
account, the total level density above 3 MeV should be a factor of two larger
than that for negative parity states alone. The overall agreement between the
calculated level density and that inferred from reaction data is reasonable.
However, the differences exhibited between the different reactions and the
fact that the inferred level densities are of the order as those computed here
suggest that each reaction might be more selective than expected and the
analysis is potentially missing states.

A proper treatment of the 1/2+ level density for Fe nuclei must take into
account particle-hole excitations beyond the $1p0f$ model space. For example,
for 57Fe we should consider the coupling of the $\nu(0g_{9/2})$ particle
orbital to the calculated level density of (4,5)+ states of 56Fe, and the
coupling of $\pi(0d_{3/2},1s_{1/2})$ hole orbitals to the calculated level
density of (0,1,2,3)+ states of 58Co. This extension will be explored in the
future.

Figure 17:  (color online) Level densities for 74Ge compared with experimental
values. The black points, labeled Ohio, are inferred from proton evaporation
spectra Voinov-2 , while the brown squares, labeled Oslo, are from the Oslo
method Ge74-Oslo . Level densities are shown for two shell model interactions,
jun45 (upper) and $jj44b$ (lower). The green and blue lines represent the
total level density for positive- and negative-parity states, respectively,
while the red line is the total level density.

Figure 18:  (color online) Level densities for 76Ge compared with experimental
values. The black points, labeled Ohio, are inferred from proton evaporation
spectra Voinov-2 , while the brown squares, labeled Oslo, are from the
$\beta$-Oslo method Ge76-Oslo . Level densities are shown for two shell model
interactions, jun45 (upper) and $jj44b$ (lower). The green and blue lines
represent the total level density for positive- and negative-parity states,
respectively, while the red line is the total level density.

In Figures 17 and 18, the ELMAC(8,100) results are shown for the nuclei 74Ge
and 76Ge within the $jj44$ shell-model space and the $jj44b$ and jun45
interactions in comparison with experimental values inferred from proton
evaporation spectra resulting from the compound nuclear reactions
68,70Zn(7Li,Xp) (black circles) Voinov-2 . In addition, for 74Ge, results
Ge74-Oslo  from the Oslo method are shown (brown squares), while for 76Ge,
results Ge76-Oslo  from the $\beta$-Oslo method are shown. Note that the Oslo
method requires a normalization, which was extracted from the experimental
$D_{0}$ value. Overall, the agreement between the ELMAC(8,100) results and
those inferred from proton-evaporation spectra are excellent up to
$E_{x}\approx 8-9$ MeV. This is well within the expectation that the shell
model provides an accurate representation of the excitation spectrum up to the
point where intruder states appear.

Table 4: Comparison between calculated and experimental mug  level spacings for $\l=1$ neutron resonances ($D_{1}$) for various Fe isotopes. The neutron separation energy, $S_{n}$, for the isotope of listed and the angular momentum $J_{t}^{\pi}$ for the target A-1Fe nucleus are shown. | $J_{t}^{\pi}$ |  $S_{n}$ (MeV) |  $D_{1}^{\rm calc}$ (keV) |  $D_{1}^{\rm exp}$ (keV)
---|---|---|---|---
55Fe | $0^{+}$ | 9.298 | 5.6 | 4.75$\pm$0.15
57Fe | $0^{+}$ | 7.646 | 7.6 | 8.21$\pm$0.48
58Fe | $\frac{1}{2}^{-}$ | 10.044 | 3.3 | 2.58$\pm$0.26
59Fe | $0^{+}$ | 9.298 | 11.6 | 5.03$\pm$0.30
Table 5: Comparison between experimental mug  and calculated (with the $jj44b$ and jun45 interactions) level spacings for $\l=0$ neutron resonances ($D_{0}$) for various Ge isotopes. The neutron separation energy, $S_{n}$, for the isotope of listed and the angular momentum $J_{t}^{\pi}$ for the target A-1Ge nucleus are shown. | $J_{t}^{\pi}$ |  $S_{n}$ (MeV) |  $D_{0}^{\rm calc}$ (keV) |  $D_{0}^{\rm calc}$ (keV) |  $D_{0}^{\rm exp}$ (keV)
---|---|---|---|---|---
|  |  | $jj44b$ | jun45 |
73Ge | $0^{+}$ | 6.782 | 6.6 | 4.3 | 2.07$\pm$0.29
74Ge | $\frac{9}{2}^{+}$ | 10.196 | 0.33 | 0.23 | 0.099$\pm$0.001
75Ge | $0^{+}$ | 6.505 | 8.9 | 5.5 | 3.0$\pm$1.5
77Ge | $0^{+}$ | 6.076 | 18.14 | 10.6 | 4.82$\pm$0.76

To conclude this section, calculated values for the level spacings for Fe and
Ge isotopes are shown in Tables 4 and 5, respectively. For Fe isotopes, level
spacing for $l=1$ neutron resonances, $D_{1}$, are shown, while for Ge
isotopes, the level spacings for $l=0$ neutron resonances are displayed. The
experimental neutron separation energy, $S_{n}$, which is equivalent to the
excitation energy of the system of interest, is tabulated as well as the
angular momentum and parity, $J_{t}^{\pi}$, of the target ${A-1}$ nucleus. The
experimental data are from Ref. mug . For the Ge isotopes, results are shown
for the two shell-model Hamiltonians $jj44b$ and jun45. Overall, good
agreement is achieved for Fe isotopes except for 59Fe, which is likely
signaling an increasing importance of the $0g_{9/2}$ orbit as more neutrons
are added. For the Ge isotopes, the calculated $D_{0}$ values are larger than
experiment. This implies that the computed level densities are too small,
which is in contradiction with the agreement with the level densities inferred
from proton-evaporation spectra as shown in Figs. 17 and 18. The jun45
interaction has a larger level density and generally yields a $D_{0}$ value
within a factor of two from experiment. The exception is 74Ge, but here
$S_{n}=10.196$ MeV, which from Fig. 17, is an excitation energy about 1-2 MeV
above where the model space is valid. On the other hand, we note the overall
good agreement between our Ge calculations and the data from Ref. Voinov-2
shown in Figs. 17 and 18.

##  VII Angular Momentum Dependence of the Level Density

The angular momentum dependence of the level density is key to understanding
many reactions. A commonly used form comes from the original work of Bethe
bethe , where the level density for a given $J$ is

$\rho(E_{x},J)=P(J)\rho(E_{x})$ (47)

with

$P(J)=\frac{(2J+1)}{2\sigma^{2}}\,\,{\rm
exp}\left[-(J+1/2)^{2}/2\sigma^{2}\right],$ (48)

and $\sigma^{2}$ being the so-called spin cutoff parameter, which is energy
dependent. The spin cutoff parameter can be determined at a fixed excitation
energy via

$\sigma^{2}=\langle(J+1/2)^{2}\rangle/2.$ (49)

The calculated spin cutoff parameters for 57Fe, 74Ge, and 76Ge as a function
of excitation energy are shown in Figures 19 and 20. In Figure 20, both the
positive- and negative-parity spin cutoff parameters are shown.

Figure 19: Calculated spin cutoff parameter for 57Fe as a function of
excitation energy.

Figure 20: Calculated spin cutoff parameter for 74Ge and 76Ge as a function of
excitation energy. The black and red lines represent the positive- and
negative-parity spaces, respectively.

Figure 21: Angular momenta probabilities for 57Fe are shown for five
excitation-energy slices across the density of states. The red lines show the
results obtained from Eq. (48) with the spin cutoff parameter shown in Figure
19.

The probability distribution of angular momenta for the three nuclei studied
are shown in Figures 21 \- 23 at five distinct excitation energies. The black
points are the probability distributions from the extrapolated Lanczos method,
while the red lines represent the results from Eq. (48) using the spin cutoff
partameters computed at each excitation energy as shown in Figures 19 and 20.
Overall, the computed angular momenta distributions are in excellent agreement
with the Bethe ansatz of Eq. (48).

Figure 22: Angular momenta probabilities for 74Ge are shown for five
excitation energies across the positive- and negative-parity states. The red
lines show the results obtained from Eq. (48) with the spin cutoff parameter
shown in Figure 20.

Figure 23: Angular momenta probabilities for 76Ge are shown for five
excitation energies for positive- and negative-parity states. The red lines
show the results obtained from Eq. (48) with the spin cutoff parameter shown
in Figure 20.

##  VIII Conclusions

We have discussed the application of the Lanczos method to the calculation of
level densities. We showed that for a given $J$ value, the $\alpha_{1-4}$ and
$\beta_{1-4}$ components of the Lanczos matrix obtained from the first four
Lanczos iterations provide the information sufficient to obtain the lowest
eight moments of the Hamiltonian with an accuracy of approximately 1%. We
derive exact but complex equations that relate these $\alpha$ and $\beta$
matrix elements to the moments. We compare the results to calculations for
matrix dimensions up to $10^{6}$ where exact results from full diagonalization
can be obtained. We also show that the uncertainty of the moments decreases
with increasing matrix dimension.

A method to extrapolate the Lanczos matrix (ELM) to the full space was
presented that made use of the first eight moments of the Hamiltonian. Level
densities were obtained with the ELM method and compared to exact shell-model
results where possible. The ELM procedure was shown to provide an excellent
representation of the asymptotic (high-energy) behavior of the level density,
and with a sufficient number of actual Lanczos iterations, the ELM method was
shown to provide excellent agreement with the exact shell-model level density.
In some cases, a discontinuity exists between the exact Lanczos iterations and
the modeled matrix elements that causes the ELM procedure to miscalculate the
level density at low excitation energies. A procedure to analytically continue
the the level density from the high-energy region to the lowest energy in the
configuration space ($E_{0}$) was presented. A calculated uncertainly of about
100(10) keV in the ground-state energy is enough to obtain the level density
above with an accuracy of approximately 10(1)% for a given model space and
Hamiltonian. The calculation of the ground-state energy to within 100 keV
requires on the order of 20 Lanczos iterations.

We compare the results of the ELM method with those obtained with the binomial
approximation that makes use of the first four moments. In some cases with the
moments close to the Gaussian limit, the two methods give similar results. But
there are other cases with the binomial method cannot be used. Finally, we
compared calculations for the level density with ELM for 57Fe and 74,76Ge
nuclei with those extracted from experiment. In addition, we computed $\ell=0$
and 1 resonance spacing, $D_{0}$ and $D_{1}$, for Fe and Ge isotopes.

###### Acknowledgements.

We gratefully acknowledges several useful discussions with S. Grimes, A. C.
Larsen, Z. Meissel, A. Voinov, and K. Wendt. This work was performed under the
auspices of the U.S. Department of Energy by Lawrence Livermore National
Laboratory under Contract DE-AC52-07NA27344 and NSF grant PHY-1811855.

##  Appendix A Solution for $a$\- and $b$-parameters

Given the set of moments $H_{1}$, $M_{2}$, and $R_{3-8}$, the strategy is then
to find an optimal set of coefficients $a_{i}$ and $b_{i}$ that reproduce
these moments. From Eqs.(32)-(IV), it is clear that the moments are highly
non-linear functions of the parameters $a_{i}$ and $b_{i}$. However, in
general, the dominant parameters will be $a_{0}$ and $b_{1}$. For example, in
the limit of a Gaussian, the odd moments are zero, and $M_{2}=-2b_{1}$. Thus,
one strategy to find the parameters is to assume that $a_{1-3}$ and $b_{2-4}$
are small and that the moments can be linearized relative to small changes in
the parameters. We start with all $a_{i>0}=0$ and solve for $b_{1}$ and
$b_{2}$ using $M_{2}$ and $M_{4}$. Note that $b_{2}$ can be isolated with the
ratio $M_{4}/M_{2}^{2}$, yielding a quadratic equation with two solutions,
with the smallest being the most realistic. With $b_{2}$ found, we then use
$M_{2}$ to fix $b_{1}$. Initial estimates for $a_{1}$ and $a_{2}$ can then be
found from the odd moments $M_{3}$ and $M_{5}$ by truncating the analytic
expressions to the leading linear terms in $a_{1}$ and $a_{2}$, yielding two
coupled linear equations:

$\displaystyle M_{3}\approx$ $\displaystyle
6b_{1}\Big{[}a_{1}(1-4b_{2})-4a_{2}(1-5b_{2})\Bigr{]}$ (50) $\displaystyle
M_{5}\approx$ $\displaystyle
120b_{1}^{2}\Bigl{[}a_{1}(3+24b_{2}^{2})-a_{2}(9+168b_{2}^{2})\Bigr{]}$ (51)

With these initial estimates, we then perform a Taylor expansion for the
moments and truncate to first order. Representing the parameters $a_{i}$ and
$b_{i}$ with the combined parameters, $p_{i}$, and using vector notation
$\vec{p}=\\{\vec{a},\vec{b}\\}$, a set of coupled linearized expressions for
the moments can be obtained, i.e.,

$M_{k}-M_{k}(\vec{p})=\sum_{i}D_{ki}\Delta p_{i},$ (52)

where $M_{k}$ is the moment for shell-model Hamiltonian and $M_{k}(\vec{p})$
is the modeled moment evaluated from Eqs. (23)-(IV) using the modeled Lanczos
matrix elements $\alpha_{i}$ and $\beta_{i}$ from Eqs. (21) and (22).
$D_{ki}=\frac{\partial M_{k}}{\partial p_{i}}$ is the derivative of the
$k^{th}$ moment with respect to parameter $p_{i}$. Under the conditions that
the non-linear terms are small, one can iteratively obtain the optimal
parameters $\vec{p}$ by solving for the shift $\Delta\vec{p}$ and updating the
derivative matrix after each iteration. In order to minimize potential effects
of non-linear terms, at each iteration a fraction of the shift is taken to
update the new values. In practice, half the new value was chosen, and the
procedure typically finds optimal solutions in approximately 20 iterations.

##  Appendix B Finding the Matching Energy for the ELM

Finding the matching energy $E_{m}$ to analytically continue the ELM
calculation of the level density is complicated by local fluctuations in the
level density due to the discrete nature of the spectrum. Thus, it is
necessary to introduce a smoothing procedure in order to make use of Eqs. (45)
and (46). In this work, we made use of a low-pass filter, or Savitsky-Golay
filter savitsky-golay , to both smooth and compute the derivative of the level
density. To first order, the Savitsky-Golay filter is essentially a least-
squares fit of polynomial of order $M$ to the data of interest over a region
of data extending $n_{L}$ and $n_{R}$ points to the left and right of the data
point of interest respectively. Here, satisfactory results were obtained by
smoothing the level density directly (not the logarithm) with $M=4$ over the
interval defined by $n_{L}=n_{R}=10$.

## References

  * (1) M. Burbidge, G. Burbidge, W. Fowler, and F. Hoyle, Rev. Mod. Phys. 29, 547 (1957).
  * (2) D. Kasen, B. Metzger, J. Barnes, E. Quataert, and E. Ramirez-Ruiz, Nature Vol. 551, 80 (2017).
  * (3) J. E. Escher, J. T. Burke, F. S. Dietrich, N. D. Scielzo, I. J. Thompson, W. Younes, Rev. Mod. Phys. 84, 353 (2012).
  * (4) W. Hauser, H. Feshbach, Phys. Rev. 87, 366 (1952).
  * (5) S. Goriely, M. Samyn, J.M. Pearson, Phys. Rev. C 75 (2007) 064312; S. Goriely, S. Hilaire, A. J. Koning, Phys. Rev. C 78 (2008) 064307.
  * (6) C. W. Johnson, S. E. Koonin, G. H. Lang, W. E. Ormand, Phys. Rev. Lett 69, 3157 (1992).
  * (7) G. H. Lang, C. W. Johnson, S. E. Koonin, and W. E. Ormand, Phys. Rev. C 48, 1518 (1993)
  * (8) S. E. Koonin, D. J. Dean, and K. Langanke, Phys. Rep. 278, 1 (1997); H. Nakada and Y. Alhassid, Phys. Rev. Lett. 79, 2939 (1997); W. E. Ormand, Phys. Rev. C 56, R1678 (1997); Y. Alhassid, S. Liu, and H. Nakada, Phys. Rev. Lett. 83, 4265 (1999); H. Nakada and Y. Alhassid, Phys. Rev. C78, 051304(R) (2008); C. Özen, Y. Alhassid, and H. Nakada, Phys. Rev. C 91, 034329 (2015).
  * (9) K. K. Mon and J. B. French, Ann. of Phys. 95,1 (1975); S. S. M. Wong and J. .B. French, Nucl. Phys. A198, 188 (1972); J. B. French and V. K. B. Kota, Ann. Rev. Nucl. Part. Sci. 32, 35 (1982); J. B. French and V. K. B. Kota, Phys. Rev. Lett. 51, 2183 (1983); Z. Pluhar and H. A. Weidenmuller, Phys. Rev. C 38, 1046 (1988); A. .P. Zuker, Phys. Rev. C 64, 021303(R) (2001).
  * (10) R. A. Sen’kov and M. Horoi, Phys. Rev. C 82, 024304 (2010)
  * (11) R. A. Sen’kov, M. Horoi, and V. G. Zelivinsky, Comp. Phys. Comm. 184, 215 (2013).
  * (12) N. Shimizu, Y. Utsuno, Y. Futamura, T. Sakurai, T. Mizusaki, and T. Otsuka, Phys. Lett. B 753, 13 (2016)
  * (13) B. Jegerlehner, arXiv:hep-lat/9612014, (1996).
  * (14) S. Yamamoto, T. Sogabe, T. Hoshi, S.-L. Zhang, T. Fujiwara, J. Phys. Soc. Jpn. 77, 114713 (2008).
  * (15) R. D. Lawson, Theory of the Nuclear Shell Model, (Clarendon, Oxford, 1980); P. J. Brussaard and P. W. M. Glaudemans, Shell Model Applications in Nuclear Spectroscopy, (North-Holland, Amsterdam, 1977).
  * (16) C. Lanczos, J. Res. Nat. Bur. Stand. 45, 252 (1950); J. H. Wilkinson, The Algebraic Eigenvalue Problem, (Clarendon, Oxford, 1965); R. R. Whitehead, A. Watt, B. J. Cole, and I. Morrison, Adv, in Nucl. Phys. 9, 123 (1977).
  * (17) A.P. Zuker, L. Waha Ndeuna, F. Nowacki, and E. Caurier, Phys. Rev. C 64, 021394(R) (2001); E. Caurier, G. Martinez-Pinedo, F. Nowacki, A. Poves, A.P. Zuker, rev. Mod. Phys. 77, 427 (2005).
  * (18) B. A. Brown and W. D. M. Rae, Nuclear Data Sheets 120, 115 (2014).
  * (19) M. Honma, T. Otsuka, B.A. Brown and T. Mizusaki, Euro. Phys. Jour. A 25 Suppl. 1, 499 (2005).
  * (20) S. Mukhopadhyay, B. P. Crider, B. A. Brown, S. F. Ashley, A. Chakraborty, A. Kumar, M. T. McEllistrem, E. E. Peters, F. M. Prados-Estévez, and S. W. Yates, Phys. Rev. C 95, 014327 (2017).
  * (21) M. Honma, T. Otsuka, T. Mizusaki, and M. Hjorth-Jensen, Phys. Rev. C 80, 064323 (2009).
  * (22) W. E. Ormand, Int. J. Mod. Phys. 14, 67 (2005).
  * (23) A. Gilbert and A. G. W. Cameron, Can. J. Phys. 43, 1446 (1965).
  * (24) B. A. Brown and W. E. Ormand, CERN Proc. 1, 21 (2019). DOI: 10.23727/CERN-Proceedings-2019-001
  * (25) A. P. Zuker, Phys. Rev. C 64, 021303 (2001).
  * (26) Atlas of Neutron Resonances Volume 1: Resonance Properties and Thermal Cross Sections Z= 1-60 by Said F. Mughabghab, Elsevier (2018) (doi.org/10.1016/B978-0-44-463769-7.00001-4).
  * (27) Data from the NNDC On-Line Data Service database as of August 2017 http://nndc.bnl.gov/nudat2/
  * (28) E. Algin, et al., Phys. Atomic Nuclei 70, 1634 (2007).
  * (29) A. Voinov, et al., Phys. Rev. Lett. 93, 142504 (2004).
  * (30) A. C. Larsen et al., J. Phys. G 44, 065005 (2017); and private communication (2018).
  * (31) A. V. Voinov, et al., Phys. Rev. C 99, 054609 (2019).
  * (32) T. Renstrøm, H.-T. Nyhus, H. Utsunomiya, R. Schwengner, S. Goriely, A. C. Larsen, D. M. Filipescu, I. Gheorghe, L. A. Bernstein, D. L. Bleuel, T. Glodariu, A. Görgen, M.Guttormsen, T. W. Hagen, B. V. Kheswa, Y.-W. Lui, D. Negi,I. E. Ruud, T. Shima, S. Siem, K. Takahisa, O. Tesileanu, T. G.Tornyi, G. M. Tveten, and M. Wiedeking, Phys. Rev. C93, 064302 (2016).
  * (33) A. Spyrou, S. N. Liddick, A. C. Larsen, M. Guttormsen, K. Cooper, A. C. Dombos, D. J. Morrissey, F. Naqvi, G. Perdikakis, S. J. Quinn, T. Renstrom, J. A. Rodriguez, A. Simon, C. S. Sumithrarachchi, and R. G. T. Zegers, Phys. Rev. Lett. 113, 232502 (2014).
  * (34) H. A. Bethe, Phys. Rev. 50, 332 (1936); Rev. Mod. Phys. 9, 69 (1937).
  * (35) A. Savitsky and M. J. E. Golay, Anal. Chem. 36, 1627 (1964); R. W. Hamming, Digital Filters, 2nd ed. (Englewood Cliffs, NJ: Prentice Hall, 1983); M. U. A. Bromba, Anal. Chem. 53, 1583 (1981); W. H. Press, S. A. Teukolsky, W. T. vettering, and B. P. Flannery, Numerical Recipes in Fortran, 2nd ed. (Cambridge University Press, New York,1992), p. 644-649.