*2.3. Block Lanczos Method*

In a large-scale problem, the number of the Hamiltonian matrix elements is too huge to be stored in the memory. In the KSHELL code, the product of the Hamiltonian matrix and a vector is the most time-consuming part of the Lanczos algorithm, since the Hamiltonian matrix elements are generated on-the-fly for every matrix-vector product. In order to reduce the time of the on-the-fly generation, the block algorithm to the thick-restart Lanczos method, named the thick-restart block Lanczos method [34].

In the block algorithm, the on-the-fly generation is performed on every product of the matrix and block vectors, not a vector. Since the number of the products in the block Lanczos method is usually much smaller than that in the Lanczos method, the number of the generations is reduced in comparison with the Lanczos method. The block Lanczos method is one of the block Krylov subspace methods. The block Krylov subspace is defined with a block of the initial vectors *V*<sup>0</sup> = (*v*0,1, *v*0,2, ··· *v*0,*q*) as

$$\mathcal{K}\_n(H, V\_0) \quad = \quad \text{span}\{V\_0, HV\_0, H^2V\_0, \dots, \cdot, H^{n-1}V\_0\},\tag{6}$$

where *q* is the block size. Although this subspace is spanned by the *nq* vectors the number of on-the-fly generation of the matrix elements is *n*. The eigenvalues of the Hamiltonian in the Krylov block subspace are a good approximation to the exact ones. In practical code, the thick-restart block Lanczos method was introduced in the KSHELL code [34].

As other eigensolvers for the shell-model calculations, the LOBPCG method [38] and the Sakurai–Sugiura method [39] have been tested. In our KSHELL code, the block Sakurai–Sugiura method and its variant, stochastic estimation of eigenvalue density, were also implemented [40,41]. To study resonance states such as the Gamow shell model, the energies of the target many-body resonant states correspond to interior eigenvalues, not the lowest one, and are complex numbers. Since the simple Lanczos method cannot be used in such a case, the Jacobi-Davidson method was adopted in the Gamow shell model [42] and the Sakurai–Sugiura method was tested in the complex scaling method [43].

#### **3. Approximation Methods to Exact Diagonalization**

While shell-model codes for large-scale calculations have been developed vigorously, various approximation schemes to the full model space were proposed. I briefly review these attempts hereafter.

The auxiliary-field Monte Carlo was applied to shell-model calculations and became one of the successful methods in describing *p f*-shell nuclei in the 1990s [44]. However, the notorious sign problem in the quantum Monte Carlo method hampers us from utilizing realistic effective interactions. As a possible solution to the sign problem, the extrapolation method [44], the complex Langevin method [45], and a constrained-path approximation [46] were suggested. Rather schematic interaction without the sign problem has been adopted for the study of nuclear level density [47].

The density matrix renormalization group (DMRG) method is known to be an efficient variational method to solve a one-dimensional quantum many-body problem. It was applied to shell-model calculations coupled to continuum states [48,49]. Its advanced one, the tensor network method, was applied later [50]. The projected shell model has been a useful model with rather schematic interaction to describe high-spin states [51]. The angular-momentum projected CI is another variational method whose basis states are the angular-momentum projected deformed Slater determinants with particle-hole excitation [52]. On top of that, a lot of effort has been paid to developing efficient truncation

schemes, such as the correlated basis method [53], the variational Monte Carlo method with random walkers in *M*-scheme basis states [54,55], nucleon-pair truncation [56,57], generator coordinate method (GCM) [58] and so on.

The variation after mean-field projection in realistic model space (VAMPIR) method is one of the successful methods to describe the *N Z* medium-mass nuclei [59,60]. In the VAMPIR method, the nuclear wave function is expressed as a linear combination of the parity, angular-momentum, number projected quasiparticle vacua. Each basis state is determined one by one by minimizing the energy expectation value of this wave function. In that sense, this method is the variation after the projection and superposition. In the hybrid multi-determinant (HMD) method, the number-projected quasiparticle vacua are replaced by the deformed Slater determinants [61,62]. The HMD has been used mainly for the no-core shell-model approach [63].

These truncation methods provide us with the variational upper limit, which may be controlled by a parameter to define the truncation. The exact energy can be estimated by the extrapolation of the energy function of this parameter. The exponential convergence method was suggested, the energy eigenvalue of the truncated subspace is expected to converge as a function of its dimension [64]. As an alternative, the energy-variance expectation value of the variational wave function can be used as a parameter for extrapolation. The energy function of the energy variance can be fitted as a polynomial function, which enables us to estimate the exact energy precisely. It was first introduced into the shell-model calculation with particle-hole truncation in Ref. [65] and applied also to no-core shell model calculations [66], the MCSM [67], QVSM [24], HMD [68], and the importance-truncated shell model (ITSM) [69].

The ITSM was introduced firstly into the no-core shell model approach by R. Roth and his collaborators successfully [70], and later applied also to the conventional shell model approach [69]. The *M*-scheme basis states are truncated by the importance measure, which is estimated by the many-body perturbation theory. The energy is calculated as a function of the criteria of the importance measure, and it is extrapolated to the full space. One of its results is shown in Section 3.2.

Among these efforts, the MCSM and the QVSM have been introduced. These models are discussed in the following Subsections.

### *3.1. Monte Carlo Shell Model*

Among all the efforts to develop various approximation schemes, the MCSM was suggested in the 1990s [22,71,72] and revised in an advanced form [23,73] in the 2010s. It has been used for extensive studies of neutron-rich nuclei. In this Section, I briefly review the theoretical framework of the MCSM.

In the MCSM the nuclear wave function is written as a linear combination of the parity, angular-momentum projected deformed Slater determinants:

$$|\Psi^{(N\_b)}\rangle = \sum\_{n=1}^{N\_b} \sum\_{K=-J}^{J} f\_{nK} P\_{MK}^{J\pi} |\phi\_n\rangle. \tag{7}$$


$$|\phi\_n\rangle = \prod\_a \left(\sum\_k D\_{k\alpha}^{(n)} c\_k^\dagger \right) |-\rangle,\tag{8}$$

which is parametrized by the complex matrix *D*(*n*) *<sup>k</sup><sup>α</sup>* . The energy *<sup>E</sup>*(*Nb*) is provided by solving the following generalized eigenvalue problem of the Hamiltonian and norm matrices:

$$\sum\_{nK} H\_{mM,nK} f\_{nK} = \quad E^{(N\_{\mathbb{H}})} \sum\_{nK} H\_{mM,nK} f\_{nK} \tag{9}$$

$$\begin{array}{rcl} H\_{mM,nK} & = & \langle \phi\_{\mathbb{m}} | H P\_{MK}^{J\pi} | \phi\_{\mathbb{m}} \rangle \\ N\_{mM,nK} & = & \langle \phi\_{\mathbb{m}} | P\_{MK}^{J\pi} | \phi\_{\mathbb{m}} \rangle \end{array}$$

where *PJ<sup>π</sup> MK* is the angular-momentum (*J*), parity (*π*) projector. *fiK* is the corresponding eigenvector. In the early stage of the history of the MCSM, many candidates of *Dk<sup>α</sup>* are generated by employing the auxiliary-field Monte Carlo technique and highly selected to lower the energy expectation value *E*(*Nb*) [22]. In the advanced version of the MCSM [23], one uses such selected basis states as initial states and optimize *Dk<sup>α</sup>* by minimizing the energy expectation value utilizing the conjugate gradient method.

The angular-momentum projection in Equation (9) is performed by a three-fold integral of the Euler angles and is time-consuming. The computation time of this part is almost proportional to the number of the discretized mesh points of this integral. Ref. [74] proposed that the Lebedev quadrature would reduce the number of the mesh points and thus the amount of the computation to 2/3. The practical code was tuned employing the technique in Ref. [75].

As a benchmark test, the ground-state energy of 56Ni with the FPD6 interaction [12] and the *p f*-shell model space is shown in Figure 1. It has been used as a target of the benchmark tests in various studies since the exact diagonalization is infeasible in the 1990s and its structure is interesting in terms of the soft closed core and shape coexistence. Its exact diagonalization was achieved in 2006 [76].

**Figure 1.** MCSM results of 56Ni with the FPD6 interaction. (**a**): MCSM energy against the number of the MCSM basis states. (**b**): extrapolation plot with the energy variance. (**c**): magnified view of (**b**). The red is the fitted curve with a second-order polynomial. See text for details.

Figure 1a shows the MCSM energy *E*(*Nb*) of the 56Ni as a function of *Nb*. The energy drops rapidly in the region *Nb* < 20 and it gradually converges giving the variational upper limit, −203.192 MeV (*Nb* = 150), which is only 6 keV higher than the exact one, −203.198 MeV. In our preceding works, the MCSM gave −203.161 MeV with *Nb* = 150 in 2010 [67], −203.152 MeV in 2001 [22], −203.100 MeV in 1998 [77]. Thus, the methodological development and progress in computational resources steadily improve the precision of the MCSM.

In this benchmark test, still the 6 keV difference with the exact value remains. In order to fill this small gap, the extrapolation method utilizing the energy variance was introduced. Figure 1b shows the energy against the corresponding energy variance, <sup>Δ</sup>*E*(*Nb*) <sup>=</sup> Ψ(*Nb*)|*H*2|Ψ(*Nb*)−Ψ(*Nb*)|*H*|Ψ(*Nb*)2. The energy variance is a good measure of how well the approximation works since the energy variance of the exact eigenstate is zero. As *Nb* increases, the energy and energy variance gradually decrease and the point smoothly approaches *y* axis. Figure 1c shows the magnified view of Figure 1b. The red line in Figure 1c is a second order polynomial fitted for these points. The *y*-intercept of the fitted line becomes the extrapolated value, which agrees with the exact one within a keV.

Figure 2 shows the overlap probability between the MCSM wave function and the exact wave function. Even at *Nb* = 1 the probability is 0.95. As *Nb* increases it approaches the unit smoothly and surpasses 0.99 at *Nb* = 7. The final value at *Nb* = 150 is 0.9998. Thus, the ground-state wave function of 56Ni is approximated by the MCSM in high precision, and, moreover, the estimation of its eigenenergy is improved by the energy-variance extrapolation. The MCSM method has been quite successful in medium-heavy nuclei mainly in *p f*-shell and neutron-rich nuclei in the medium-heavy mass region. However, the MCSM tends to underestimate the pairing correlation and 2<sup>+</sup> excitation energies. In order to treat the pairing correlation properly, the QVSM was introduced and is discussed in the next Subsection.

**Figure 2.** Overlap probability between the MCSM wave function and the exact wave function by the Lanczos method against the number of the MCSM basis states *Nb*. These wave functions were computed for the ground state of 56Ni with the FPD6 interaction [12].
