## Abstract

We present a newly enhanced version of the Monte Carlo shell-model (MCSM) method by incorporating the conjugate gradient method and energy-variance extrapolation. This new method enables us to perform large-scale shell-model calculations that the direct diagonalization method cannot reach. This new-generation framework of the MCSM provides us with a powerful tool to perform very advanced large-scale shell-model calculations on current massively parallel computers such as the K computer. We discuss the validity of this method in ab initio calculations of light nuclei, and propose a new method to describe the intrinsic wave function in terms of the shell-model picture. We also apply this new MCSM to the study of neutron-rich Cr and Ni isotopes using conventional shell-model calculations with an inert ^{40}Ca core and discuss how the magicity of *N*=28,40,50 remains or is broken.

## Introduction

Understanding of the structure of all nuclei from first principles, usually called the ab initio nuclear calculation, is one of the ultimate goals in nuclear theory. For this purpose, one usually starts with the nucleon degree of freedom, i.e., protons and neutrons as the building blocks of a nucleus, assuming the free (or bare) nucleon–nucleon force to be the interaction between nucleons. Recent ab initio calculations have often included not only the two-nucleon force but also the three-nucleon force, which is specific to the interaction between composite particles like nucleons. Several ab initio methods have been proposed, and their accuracy has been investigated in great detail, for instance, in terms of the binding energy of the four-nucleon system [1]. It remains, however, rather difficult or infeasible to go beyond systems with $$A \gtrsim 12$$, where *A* is the number of nucleons. This is largely due to strong nucleon–nucleon correlations that require a large number of single-particle states to be included for the description of many-body states in ground states or low excited states. Although this problem can be resolved to a certain extent by using effective interactions based on various renormalization techniques, the number of *many-body* states to be included remains prohibitively large in most cases, and increases exponentially with the nucleon number. In this paper, we present a new version of the Monte Carlo shell model (MCSM), demonstrating how it can contribute to the ab initio nuclear calculation as well as to conventional but quite huge shell-model calculations.

The nuclear shell model is known to have been successful in describing the structure of atomic nuclei based on nuclear forces. It dates back to 1949 when Mayer and Jensen discovered the magic numbers as a consequence of a mean potential including spin–orbit coupling [2, 3]. While the original concept of the shell model of Mayer and Jensen is basically an independent particle picture, the concept has been modified and extended significantly over decades since then. The modern shell model uses a sufficiently large number of superposed many-body basis states, utilizing properly constructed effective interactions so as to provide us with accurate many-body eigenstates. The many-body basis state is a Slater determinant, usually with single-particle states taken from a harmonic oscillator potential.

In conventional shell-model calculations, an inert core is assumed: all single-particle states below a given magic number are completely occupied, forming a closed shell. Single-particle states between this magic number and the next magic number constitute a *valence shell*, and nucleons in this shell are called valence nucleons. In the conventional shell model, only valence nucleons are activated. The single-particle states of activated nucleons in the shell-model calculation are called the *model space*. The model space is a concept for calculation, and is the same as some valence space in many cases. But, in other cases, the model space can be taken wider or smaller than the relevant valence shell depending on some interest or limitation. We thus distinguish model space from valence shell hereafter. Note that the model space is a more computational concept.

An effective interaction is obtained for valence nucleons, and is defined for each model space. Effects from the inert core, those from virtual excitations from the inert core, and those from virtual excitations to states above the model space are assumed to all be incorporated into this effective interaction by renormalizing it appropriately. This can be a very important issue, but the conventional shell model assumes that there is such an interaction, while its determination can be phenomenological.

It has been shown that many nuclear properties with 8≤*N*(*Z*)≤20 are excellently described with the conventional shell-model calculation in the 1*s*-0*d* (model) space (often called the *sd* shell) [4, 5], and those with $$20\le N(Z) \lesssim 32$$ are also quite well described in the 1*p*–0*f* (model) space (often called the *pf* shell) [6, 7, 8]. Note that the effective interactions used in those calculations are hybrid products of microscopic derivation and empirical adjustments.

While the conventional shell model assumes a relatively limited model space as exemplified above, the same computational procedure is applicable to the ab initio calculation when no inert core is assumed and a large number of single-particle states are taken, so that the calculation becomes close to calculations in the whole Hilbert space. Here, we refer to this ab initio method within the shell model as the ab initio shell model, an example of which is the no-core shell model (NCSM) [9, 10, 11, 12], a famous model that has been used with great success.

Whether the ab initio shell model or the conventional shell model (with a core) is used, all one has to do in computation is to diagonalize a Hamiltonian matrix spanned by all the possible many-body states in a given model space (i.e., single-particle space of activated nucleons). The number of relevant many-body basis states determines the dimension of this matrix. This dimension is often called the shell-model dimension, and causes a serious computational issue, as we shall see later. The number of many-body states is _{Np}*C*_{np}×_{Nn}*C*_{nn} without symmetry consideration, where *N*_{p} (*N*_{n}) is the number of proton (neutron) single-particle states taken and *n*_{p} (*n*_{n}) is the number of protons (neutrons) activated. It roughly increases exponentially with *N*_{p} (or *N*_{n}) and also with *n*_{p} (or *n*_{n}). Hence, without even considering the ab initio shell model with a huge number of *N*_{p} and *N*_{n}, the conventional shell model for heavier nuclei already suffers from the huge dimensionality of the Hamiltonian matrix necessary to tackle the whole nuclear chart, because the number of valence orbits (*N*_{p} or *N*_{n}) increases for heavier nuclei. For instance, the dimension needed for *A*∼80 *N*=*Z* nuclei is estimated to be ∼10^{27} without any symmetry consideration [13], when the 1*p*–0*f*–2*s*–1*d*–0*g* orbits are assumed to be the valence shell. Although this number can be reduced by 1–2 orders of magnitude by choosing only the states of the same *J*_{z} as are usually taken (i.e., *M*-scheme calculation), it is still far from the current computational limit of 10^{10−11} in the *M*-scheme.

In order to go beyond the computational limit of the shell model, a new method for performing the shell-model calculation, called the Monte Carlo shell model (MCSM), has been developed since 1996 [14], guided by the quantum Monte Carlo diagonalization method [15]. The MCSM utilizes the idea of the auxiliary field Monte Carlo method that is taken in the shell model Monte Carlo method [16], but this method and the MCSM are completely different. The MCSM aims to represent many-body states with a small number of highly selected many-body basis states. In this sense, the MCSM is regarded as an “importance truncation” of the entire many-body space [17]. The basis state should be represented in a compact form (i.e., a wave function with a small number of parameters), it should be able to approximate the nuclear many-body state efficiently, and its matrix elements should be calculated easily. To meet this demand, we usually use deformed Slater determinants, and parity and total angular momentum projections are operated on them if needed. Otherwise, a pair-condensed state [18] or a quasi-particle vacuum state, which is used in the Hartree–Fock–Bogoliubov calculation and also taken as the basis state of the VAMPIR (Variational After Mean field Projection In Realistic model spaces) method [19], is a good candidate when required. Note that the VAMPIR method has been developed with more emphasis on calculating the energy spectra rather than improving the energy of a specified state. The MCSM basis states are added one by one, and each basis state to be added is selected among many candidate Slater determinants generated stochastically so that the energy of the state under consideration can be as low as possible. This step is repeated until the energy converges sufficiently. Thus, this original MCSM is characterized as “stochastic” in terms of the method of basis variation, and as “sequential” in terms of the procedure of basis variation.

From the computational point of view, the MCSM has an advantage over the conventional diagonalization method in tolerance to the increase of the model space and the particles. When it is assumed that the number of basis states needed in a MCSM wave function is fixed, the total computational cost is scaled by the cost of each Hamiltonian matrix element. This is roughly proportional to (*N*_{sp})^{α} with *α*∼3–4 in the case of *N*_{sp}=*N*_{p}=*N*_{n}, being much weaker than the exponential increase. This advantage enabled us to perform the full *pf*-shell calculation in ^{56}Ni (with ∼10^{9}*M*-scheme dimension) with good accuracy [20] several years before the exact diagonalization was carried out [19]. See the review paper [17] for more details of early achievements.

Recently, the computational environment has been changing rapidly. The number of available CPU cores is expanding to be typically in the range of tens of thousands for world-leading supercomputers, as compared to a few tens to hundreds of CPU cores used in the early MCSM calculations. The K computer will contain more than 700 000 cores upon its completion in the autumn of 2012. This situation has strongly motivated us to renew the MCSM method to be suitable for using up-to-date massively parallel computers, in addition to wider applications including ab initio shell-model calculations. Since 2009, we have developed a new MCSM package, which includes not only renewed code, but also a novel methodology and numerical algorithm [21]. Among the methodological advancements, the most important is the introduction of the energy-variance extrapolation method [22]. In the original MCSM, the energy of a many-body state is evaluated directly from the energy expectation value of the MCSM wave function, which must be higher than the exact solution. It was quite hard to know the difference between this value and the exact value. The energy-variance extrapolation method serves as a powerful tool to pin down accurately the exact solution from a series of approximated solutions. Another methodological change is incorporating a variational aspect into the MCSM by varying the basis state, which enhances the lowering of calculated energy values. Equipped with the variational-type improvement, the MCSM is now characterized as “deterministic” in terms of the method of basis variation. Thus, keeping the original idea, the present MCSM, associated with several advancements, can be regarded as a new generation.

In this paper, we review the outline of the new-generation MCSM and show some of its earliest applications that were not feasible with the original MCSM. This paper is organized as follows. In Sect. 2, the outline of the new-generation MCSM and its feasibility are presented. We also discuss the adaptability of the MCSM to massively parallel supercomputers. In Sect. 3, application to the ab initio shell model is demonstrated. Along with numerical success, the intrinsic shape and the clustering of light nuclei are also discussed with the MCSM wave function. In Sect. 4, application to the neutron-rich chromium and nickel isotopes is presented as a case of medium–heavy mass nuclei. This region is being intensively studied in radioactive isotope facilities over the world, and also challenges nuclear models because several shapes coexist and are mixed, which results in a rapid change of the dominant shape over the isotopes. In Sect. 5, we summarize this paper and give an outlook and future perspective of the MCSM toward the launch of the K computer.

We note here that many new features may appear in the structure of exotic nuclei, because unbalanced numbers of protons and neutrons may create situations where unknown or little-known aspects of nuclear forces become visible and produce large impacts. Shell evolution due to the tensor force [23, 24, 25] and the three-body force [26] are some examples. Note that shell evolution explains the basic trends of single-particle properties; one needs comprehensive calculations to obtain physical quantities and look into correlations in depth. The exploration of such unknown features needs a theoretical framework directly linked to nuclear forces, and we expect a large contribution from the new-generation MCSM.

## Outline of the new-generation Monte Carlo shell model

The new-generation MCSM can be divided roughly into two stages: a variational procedure to obtain the approximated wave function and an energy-variance extrapolation utilizing the obtained wave function. We briefly describe these two parts in Sects. 2.2 and 2.3, respectively. The shell-model Hamiltonian and the form of the variational wave function used in the MCSM is shown in Sect. 2.1. An additional improvement to make the extrapolation procedure stable is demonstrated in Sect. 2.4. The numerical aspects of the framework that are mainly for massively parallel computation are discussed in Sect. 2.5.

### Shell-model Hamiltonian and variational wave function

In conventional nuclear shell-model calculations assuming an inert core, we use a general two-body interaction:

*i*.

*H*

^{(1)}is a one-body Hamiltonian with single-particle energies

*t*

_{i}, and

*H*

^{(2)}is a two-body interaction which has parity and rotational symmetry and is represented by so-called two-body matrix elements (TBMEs) [27].

In the case of ab initio shell-model calculations, the Hamiltonian is taken as

*T*and

*T*

_{CM}are the total kinetic energy and the kinetic energy of the center-of-mass motion, respectively. Note that

*T*

_{CM}has both one-body and two-body components. The

*V*represents two-nucleon interaction, e.g., the JISP16 interaction [29, 30, 31] in Sect. 3. In this work, because we do not explicitly treat three-nucleon forces, both the Hamiltonians of these two kinds of shell-model calculations consist of one-body and two-body interactions.

If necessary, the removal of spurious center-of-mass motion can be done by utilizing the prescription of Gloeckner and Lawson [32] to suppress the contamination of the center-of-mass motion. In this prescription the Hamiltonian to be diagonalized is replaced by

with*H*

_{CM}being defined as

**R**and

**P**are the position and momentum of the center of mass. By taking

*β*

_{CM}large enough, 〈

*H*

_{CM}〉 is suppressed as a small value.

In the framework of the MCSM, the approximated wave function is written as a linear combination of angular-momentum- and parity-projected Slater determinants,

*N*

_{b}is the number of the Slater-determinant basis states. The $$P^{I\pi }_{MK}$$ operator is the angular-momentum and parity projector defined as

*Ω*≡(

*α*,

*β*,

*γ*) are the Euler angles and $$D^I_{MK}(\Omega)$$ denotes Wigner’s

*D*-function.

*Π*stands for the parity transformation. Each |

*ϕ*

_{n}〉 is a deformed Slater determinant defined in Eq. (A1). The coefficients $$f^{(N_b)}_{n,K}$$ are determined by the diagonalization of the Hamiltonian matrix in the subspace spanned by the projected Slater determinants, $$P^{I\pi }_{MK} | \phi _n \rangle $$, with −

*I*≤

*K*≤

*I*and 1≤

*n*≤

*N*

_{b}. This diagonalization also determines the energy,

*E*

_{Nb}≡〈

*Ψ*

_{Nb}|

*H*|

*Ψ*

_{Nb}〉, as a function of

*N*

_{b}. Note that the dimension of the subspace is (2

*I*+1)

*N*

_{b}, not

*N*

_{b}. The Slater determinant basis,

*D*

^{(n)}, is given by variational calculation to minimize

*E*

_{Nb=n}. We increase

*N*

_{b}until

*E*

_{Nb}converges enough, or the extrapolated energy converges. The strategy of this variational calculation will be discussed in the next subsection.

### Variational procedure

In this section, we discuss an efficient process for the determination of *D*^{(n)} in Eq. (A1), to minimize *E*_{Nb}, and the history of its development. In the original MCSM calculation, the basis states are selected from many candidates generated stochastically utilizing the auxiliary-field Monte Carlo technique, the details of which were discussed in Ref. [17]. In order to determine *D*^{(n)} in Eq. (2.5) efficiently, M. Honma et al. introduced the few-dimensional basis approximation [33, 34], in which they adopted the steepest gradient method. They succeeded in estimating the energy of *pf*-shell nuclei with a relatively small number of bases; however, the direction of the gradient is not necessarily the most efficient choice to reach the local minimum and it is difficult to control the step width in the gradient direction in practice. Mainly because of this problem, the number of basis states was limited to be rather small (*N*_{b}≃30) compared to that in the MCSM (*N*_{b}≃100). On the other hand, W. Schmid et al. adopted the quasi-Newton method for energy minimization in the VAMPIR approximation to obtain optimized quasi-particle vacua as basis states [35]. G. Puddu demonstrated that the quasi-Newton method also works in the Slater determinant bases. In the quasi-Newton method, the step width is automatically determined to minimize the energy along the modified gradient direction [36, 37]. As a result, the step width is relatively large in early steps of the quasi-Newton procedure along the sophisticated gradient direction, so that the additional basis state has sufficiently linearly independent components from the other basis states.

In the present work, we adopt the conjugate gradient (CG) method [38, 39], which also includes linear minimization in the modified gradient direction. In comparison to the quasi-Newton method, the CG method is expected to be advantageous to save memory usage and improve computational efficiency for parallel computation, because a large Hessian matrix is utilized in the quasi-Newton method, which is not used in the CG method. Note that the energy of the linear combination of projected Slater determinants is optimized, not that of the unprojected Slater determinants. In this sense, this variational procedure is “variation after projection and configuration mixing”.

We propose the following four steps to optimize the parameters of deformed Slater determinants:

Monte Carlo sampling utilizing the auxiliary field technique (the original MCSM),

Sequential optimization for each basis state (SCG),

A repeated refinement process for each basis state (refinement),

Full (simultaneous) optimization of all basis states (FCG).

In the first step, we perform the original MCSM procedure to obtain the approximated wave function, which can be used as an initial state of the CG process. The details of the original MCSM are not discussed here, but in Ref. [17]. In Refs. [17, 22], we select the best basis states from an order of 1000 candidates, which is generated stochastically. However, the necessary number of candidates can be greatly suppressed if we proceed to the next step.

In the second step, called the SCG method, we perform a variational calculation, with variational parameters being *D*^{(i)}, sequentially by minimizing the *E*_{i} without changing the other bases, *D*^{(1)},*D*^{(2)},…,*D*^{(i−1)}, which are already fixed. In other words, we perform variational calculations with a set of variational parameters, *D*^{(i)}, sequentially. We increase the number of basis states, *N*_{b}, until the energy reaches sufficient convergence.

In the third step, which is called the refinement process, we take an initial state from the result of the SCG calculation. In the SCG calculation, *D*^{(1)} is not the best optimized parameter to minimize *E*_{Nb}, since *D*^{(1)} is determined to minimize *E*_{1}. Then, we first fix the number of basis states, *N*_{m}, and restart to minimize the energy *E*_{Nm} by the CG method to optimize each basis state from *D*^{(1)} to *D*^{(Nm)}, one by one. We iterate a few times to perform a CG process routine for all basis states one by one to get a better energy calculation.

In the fourth step, called the FCG, we fix *N*_{m} at the beginning, and determine all coefficients *D*^{(Nb)} with 1≤*N*_{b}≤*N*_{m} by minimizing the *E*_{Nm} at once. In other words, we perform a variational calculation with a set of all variational parameters, *D*^{(Nb)}, simultaneously to minimize the *E*_{Nm}. Its initial state is generated by the refinement process in order to save computation time. In principle, the FCG yields the closest energy to the exact one within a fixed number of basis states, *N*_{m}, and the refinement process also reaches the energy provided by the FCG with a large number of iterations. The refinement process needs far less time than the FCG, and provides us with an approximation as good as the solution in the FCG.

Figure 1 shows the convergence pattern of the ground-state energy of ^{72}Ge in the *f*_{5}*pg*_{9} shell, which consists of the 0*f*_{5/2}, 1*p*_{3/2}, 1*p*_{1/2}, and 0*g*_{9/2} orbits. Its *M*-scheme dimension is very large and amounts to 140 050 484. However, it can be handled by the recent shell-model diagonalization code (T. Mizusaki, N. Shimizu, Y. Utsuno, and M. Honma, code MSHELL64, unpublished). The energy eigenvalue *E*_{Nb}=〈*Ψ*_{Nb}|*H*|*Ψ*_{Nb}〉 is plotted against the number of basis states, *N*_{b}. The SCG method attains faster convergence than the original MCSM: the SCG gives the same energy as the original MCSM in almost half the number of basis dimensions. The SCG method enables us to compute the energy variance in less computation time than the original MCSM, since the time to compute the energy variance is proportional to *N*_{m}(*N*_{m}+1)/2.

In the case of the FCG, the number of basis states of the variational wave function is fixed to *N*_{m}=32 and we plot energy expectation values in the subspace spanned by the *N*_{b} basis states. In the FCG, all basis states, *D*^{(1)},*D*^{(2)},…, and *D*^{(32)}, are optimized to minimize *E*_{32}, while in the SCG wave function *D*^{(1)} is optimized to minimize *E*_{1}. Thus, the calculated *E*_{1} of the FCG is much higher than the *E*_{1} of the SCG and the *E*_{1} of the original MCSM. However, *E*_{32} of the FCG is lower than that of the SCG and the original MCSM, which means the FCG provides us with the best approximation for a fixed number of basis states. Therefore, we should discuss the energy convergence as a function of *N*_{m} in the FCG, not as a function of *N*_{b}, unlike the case of the SCG.

While the FCG yields the best energy with fixed *N*_{m}, the FCG iteration needs many computational resources, since we have to replace all basis states at every FCG iteration. Which optimization level is most efficient in terms of the computation time depends on the property of the objective wave function. For applications, we adopt the original MCSM method, the SCG method in Sect. 3, and the refinement process in Sect. 4.

### Energy-variance extrapolation

Even if the variational procedure works excellently, a small gap between the exact energy eigenvalue and the energy expectation value of the approximated wave function remains. The gap can be removed by extrapolation procedures, which have been studied intensively [40, 41, 42, 43, 44, 45]. Among them, the energy-variance extrapolation method is a general framework for the supplementation of the variational calculation and is rather independent of the form of the approximated wave function. This method is known in condensed matter physics [46] and was first introduced in shell-model calculations with conventional particle-hole truncation by T. Mizusaki and M. Imada [42, 43]. In the present work, we apply this extrapolation method to estimate the exact eigenvalue precisely by utilizing a sequence of the linear combinations of the projected Slater determinants, |*Ψ*_{Nb}〉 with 1≤*N*_{b}≤*N*_{m}, which have been obtained in the variational procedure discussed in Sect. 2.2. The major obstacle to energy-variance extrapolation is the need to compute the energy variance. Its efficient computation is described in Sect. 2.5.

The energy-variance extrapolation method is based on the fact that the energy variance of the exact eigenstate is zero. The energy variance of the approximated wave function is not exactly zero, but is rather small and approaches zero as the approximation is improved. In the framework of the energy-variance extrapolation, we draw the energy *E*_{Nb}=〈*Ψ*_{Nb}|*H*|*Ψ*_{Nb}〉 against the energy variance $$\langle \Delta H^2 \rangle _{N_b} = \langle \Psi _{N_b} | H^2 |\Psi _{N_b} \rangle - E_{N_b}^2$$, the plot of which is called an “energy-variance plot”. The variance usually approaches zero as *N*_{b} increases, as the point in the energy-variance plot approaches the *y*-intercept. Following the idea of Refs. [42, 43], these points are usually fitted by a first- or second-order polynomial such as

*c*

_{0},

*c*

_{1}, and

*c*

_{2}are determined by least-squares fiting. By extrapolating the fitted curve into the

*y*-intercept we obtain the extrapolated energy, namely,

*c*

_{0}.

In the framework of the present study, the variational procedure discussed in Sect. 2.2 provides us with a sequence of approximated wave functions, which can be utilized in the energy-variance extrapolation method. The SCG procedure provides us with a successive sequence of the wave functions, |*Ψ*_{Nb}〉 with 1≤*N*_{b}≤*N*_{m}, where *N*_{m} is the maximum of *N*_{b}. For each *N*_{b}, we evaluate the energy *E*_{Nb} and energy variance 〈*ΔH*^{2}〉_{Nb}. Here, we demonstrate how the extrapolation method works with ^{56}Ni in the *pf* shell as an example. The effective interaction is the FPD6 interaction [47], which was also used in Ref. [22].

Figure 2(a) shows the energy-variance plot of the SCG method for yrast states of ^{56}Ni. In this calculation, we take the *K*=*I* state only in the angular-momentum projector, namely $$f_{n,K}^{(N_b)}=0$$ if *K*≠*I* in Eq. (2.5), for simplicity. As *N*_{b} increases, the energy-variance point moves smoothly and approaches the *y*-axis or variance zero, except for the 8^{+} state. The fitted curves for these points are shown as red solid lines, and they also show smooth behavior. The extrapolated energies, or *y*-intercepts of the fitted curves, agree quite well with the exact ones, which are shown as open circles on the *y*-axis. Especially concerning the ground-state energy, the minimum variance of the approximated wave function is 〈*ΔH*^{2}〉_{Nb=100}=0.89 MeV^{2}, which is smaller and gives a better approximation than the result of the original MCSM, 〈*ΔH*^{2}〉_{Nb=150}=1.05 MeV^{2} [22], mainly thanks to the introduction of the CG method.

On the other hand, the plot of the 8^{+} state shows anomalous behavior, in which the energy decreases when increasing the number of basis states but the variance does not decrease at 〈*ΔH*^{2}〉≃4 MeV^{2}. As a result, the simple extrapolation method apparently fails. A straightforward solution to this failure is to increase *N*_{b} until the extrapolation method works; this is shown in Ref. [28]. However, in other cases, it might be difficult due to the increase in computation time. We discuss another remedy for such cases in the following section.

### Reordering technique to improve the energy-variance extrapolation

The anomalous behavior can be removed by the reordering of the basis states [48]. In this section, we demonstrate that the reordering technique makes the energy-variance extrapolation stable and avoids the difficulty of an anomalous kink, such as the 8^{+}_{1} state of ^{56}Ni discussed in the previous subsection.

A sequence of the approximated wave functions (|*Ψ*_{1}〉, |*Ψ*_{2}〉,…,|*Ψ*_{Nb}〉,…,|*Ψ*_{Nm}〉) is specified by a set of basis states and its order (|*ϕ*_{1}〉, |*ϕ*_{2}〉, …, |*ϕ*_{Nb}〉, …, |*ϕ*_{Nm}〉). In the SCG method, the order of basis states is determined by the variational procedure. However, we can shuffle the order of basis states and make another sequence of the approximated wave functions without additional computational effort. If we assume that the extrapolated value is independent of the order of basis states, there exists an optimum order that makes the extrapolation procedure stable. In the reordering method, the order is determined so that the gradient of the curve in the energy-variance plot is as small as possible. The practical algorithm to determine the order is discussed in Ref. [48].

Figure 2(b) shows an extrapolation plot with the reordering technique using the same set of basis states as in Fig. 2(a). In this case, the gradient of the fitted curve is so stable that the points in the energy-variance plot are fitted by a first-order polynomial and the region to be used for the fit is rather broad. In Fig. 2(b), the anomalous behavior of the variance plot of the 8^{+} state disappears and the extrapolated value agrees with the exact energy quite well. ^{56}Ni is also known to have a shape coexistence feature [49], which is plausible as the origin of the anomalous behavior of the 8^{+} state in Fig. 2(a), like the case of ^{72}Ge discussed in Ref. [48].

We demonstrate another example of the energy-variance extrapolation method combined with the reordering technique in Fig. 3 using the FCG wave functions of ^{72}Ge with the JUN45 effective interaction [50]. In this figure, the energy expectation value of the FCG wave function is plotted as a function of the corresponding energy variance with the basis states being *N*_{m}=8,16,24, and 32, and their second-order fitted curves. The order of the basis states in each sequence is determined by the reordering technique[48]. Note that there is no specific order in the FCG procedure because we treat all basis states on an equal footing. The extrapolated energy apparently converges except for *N*_{m}=8. While we show the second-order fit in the figure, even the extrapolated value of the first-order fit agrees well with that of the second-order fit in the case of *N*_{m}=32.

### Computational aspects

The large-scale shell-model calculation is one of the challenging issues for nuclear physics. It is essential to develop a program which runs efficiently on recent supercomputers. Concerning a calculation using a single CPU, we proposed a sophisticated method of efficient computation of the matrix element of non-orthogonal Slater determinants in Ref. [51]. Moreover, since the main trend of recent supercomputers favors massively parallel computers, parallel efficiency is worth discussing for future studies.

In the case of the SCG method, we need to compute the Hamiltonian matrix elements and the gradient vector concerning *D*^{(n)} of two angular-momentum-, parity-projected Slater determinants. The three-dimensional integral of the Euler angles in Eq. (2.6) is performed by discretizing each range of the integral into mesh points using the Gauss–Chebyshev quadrature for the *z*-axis rotations and the Gauss–Legendre quadrature for the *y*-axis rotation, which is shown in Eq. (A3). The numbers of mesh points are taken as e.g. 26 for the *z*-axis rotation and 16 for the *y*-axis rotation. The parity projection is equivalent to two mesh points. The product of these mesh points, which is denoted by *λ* in Eq. (A3), gives rise to *N*_{mesh}=26^{2}×16×2=21 632 components in evaluating a matrix element in terms of the projected Slater determinants. Since these components can be computed independently, the program was written for massive parallel computation. When we apply the matrix-product technique discussed in Ref. [51] with the bunch size *N*_{bunch} being e.g. 30, we still have the *N*_{b}*N*_{mesh}/*N*_{bunch}≃721*N*_{b} elements to be computed in parallel.

Figure 4 shows the parallel efficiency of the benchmark calculation of the ground state of ^{64}Ge as an example. The model space consists of the *pf* shell and *g*_{9/2} orbit; the PFG9B3 effective interaction is used (M. Honma et al., unpublished). Its *M*-scheme dimension reaches 1.7×10^{14}, which is far beyond the current limitation of the Lanczos method. The MCSM result of this system has already been reported in Refs. [22, 48]. This benchmark was performed using the Intel Fortran compiler, version 11.0 [52], on the T2K open Supercomputer at the University of Tokyo [53].

Figure 4 (a) shows the performance gain with parallel computation of the SCG process to determine the 1st, 4th, 8th, 16th, and 32nd basis states respectively. Although the parallel efficiency for calculating the first basis state is not good because of the small amount of computation, the efficiency for the 32nd basis state with 2048 CPU cores reaches 82% of that with 16 cores.

We calculate the energy variance using the formula shown in Appendix A [22]. Because the two-body matrix elements in the *M*-scheme, *v*_{ijkl}, are sparse due to the symmetry that the Hamiltonian has, we store in memory only non-zero matrix elements in block-diagonal form by treating *v*_{(ij),(kl)} as a rank-2 matrix with indices (*ij*) and (*kl*). Thus, we can compute the energy variance efficiently; the details of the practical computation are given in Appendix. In a similar manner to the case of the variational process, we compute the energy variance by dividing the whole computation into matrix elements, which are moreover divided into each mesh point of Eq. (2.6), resulting in *N*_{b}(*N*_{b}+1)×*N*_{mesh} independent components to be computed in parallel. In addition, we do not need an iterative process like the CG method, and therefore a small amount of network communication appears only at the beginning and at the end of the computation. Thus the performance scaling of the parallel computation seems perfect at *N*_{b}≥4, which is shown in Fig. 4(b).

In practice, it took 807 seconds in total to obtain an SCG wave function of ^{64}Ge 0^{+}_{1} state with 32 basis states, and it took 588 seconds to compute the energy variance of this SCG wave function using 2048 CPU cores.

## Application to the ab initio shell model

In this section we focus on the latest application of the MCSM to ab initio shell-model calculations, which has recently become feasible with the aid of the major development of the MCSM algorithm discussed in Sect. 2 and also a remarkable growth in the computational power of state-of-the-art supercomputers. First, the no-core shell model (NCSM) and its variants are briefly reviewed. The limitation of the NCSM and the motivation for the application of the MCSM to the ab initio no-core full configuration interaction (FCI) approach are further discussed here. Then, the current status of the benchmarks in the no-core MCSM is referred to, mostly based on the results from Ref. [54]. Finally, our challenge to visualize the intrinsic states constructed by superpositioned non-orthogonal Slater determinants is also demonstrated.

### Ab initio shell models

One of the major challenges in nuclear theory is to understand nuclear structure and reactions from ab initio methods. Ab initio calculations for nuclear many-body systems beyond *A*=4 have recently become feasible due to the rapid evolution of computational technologies these days. In ab initio approaches for nuclear structure calculations, all the nucleons constituting the nucleus are considered as fundamental degrees of freedom and bare/effective interactions based on realistic nuclear forces are adopted. As for bare two- and three-nucleon interactions, the phase-shift equivalent family of two-nucleon interactions, derived from the meson-exchange theory and chiral effective field theory, in addition to three-nucleon interactions [55, 56, 57, 58, 59, 60] is usually used.

Ab initio NCSM has been emerging for about a decade and is now available for the study of nuclear structure and reactions in *p*-shell nuclei [9, 10, 11, 12]. Unlike the conventional shell model with a core, the NCSM does not assume an inert core, just as the name itself implies, and treats all the nucleons composing the nucleus on an equal footing. The NCSM is thus said to be one of the ab initio approaches along with Green’s function Monte Carlo [61, 62, 63, 64] and coupled cluster theories [65]. In the NCSM (in a narrow sense)[9, 10, 11, 12], the model (or basis) space is usually truncated by the so-called $$N_{\max }$$, which is the sum of the excitation quanta above the reference state. The effective interactions renormalized to that model space are used so as to obtain the faster convergences of the energy with respect to $$N_{\max }$$. Generally, the effective interactions are derived by the so-called Lee-Suzuki method [66, 67, 68, 69]. The NCSM result approaches the exact solution either by taking the larger model space with the level of the cluster approximation fixed or by improving the order of the cluster expansion with the model space fixed.

A similar but distinct approach to the NCSM is the no-core full configuration (NCFC) approach [81, 82, 83]. The NCSM result by using the effective interactions derived by the Lee-Suzuki procedure approaches the exact solution either from below or above due to the violation of the strict variational upper bound of the exact solution. Therefore, extrapolation of the NCSM result into infinite model space is obscure. The NCFC method employs bare or effective (softened) low-momentum interactions evolved from bare nuclear forces by the renormalization group transformations [70], which validates the variational upper bound of the calculated energy. The NCFC enables access to full ab initio solutions by a simple extrapolation into infinite model space in two-dimensional parameter space ($$\hbar \omega $$, $$N_{\max }$$). One of the advantages in both the NCSM and NCFC methods is the perfect factorization of the intrinsic and center-of-mass wave functions, so that the intrinsic state does not suffer from spurious center-of-mass motion.

As ab initio approaches treat all of the nucleons democratically, computational demands for the calculations explode exponentially as the number of nucleons and/or the model spaces increase. The current limitation of the direct diagonalization of the Hamiltonian matrix by the Lanczos iteration is around the order of 10^{10} as shown in Fig. 5. So far the largest calculations have been done in ^{14}N with $$N_{\max } = 8$$, which results in the *M*-scheme many-body matrix dimensions being ∼10^{9} and associated non-vanishing three-nucleon force matrix elements being ∼4×10^{13} [71]. In order to access heavier nuclei beyond the *p*-shell region with larger model spaces by ab initio shell-model methods, many efforts have been devoted for several years. One of these approaches in the $$N_{\max }$$ truncation is the importance-truncated NCSM (IT-NCSM) [72, 73]. In the IT-NCSM, the model spaces are extended by using the importance measure evaluated by the perturbation theory. Another approach is the symmetry-adapted NCSM (SA-NCSM) [74, 75, 76], where the model spaces are truncated by the selected symmetry groups.

Besides the $$N_{\max }$$ truncation of the model space in the ab initio shell models, the FCI method can give exact solutions in the fixed model space. Unlike the $$N_{\max }$$ truncation in the NCSM and NCFC methods, the FCI truncates the model space by the single-particle states, so-called *N*_{shell} or $$e_{\max } (\equiv N_{\mathrm {shell}} - 1)$$. As shown in Fig. 5, the explosion of the dimensionality prohibits full ab initio solutions of the FCI (and also the NCSM) beyond the lower *p*-shell region. Like the attempts of the IT-NCSM and SA-NCSM, the MCSM is a promising candidate to go beyond the FCI method [54, 77]. Note that there is a similar approach to the no-core MCSM, referred to as the hybrid multi-determinant method [78]. In the following subsection we will show some recent investigations by the ab initio no-core MCSM.

### Benchmarks for the MCSM to the ab initio no-core FCI

In some exploratory work, the original MCSM has been applied to no-core calculations for the structure and spectroscopy of beryllium isotopes [79]. The low-lying excited states of ^{10}Be and ^{12}Be have been investigated. The excitation energies of the first and second 2^{+} states and the B(E2; 2$$^+_{1}\rightarrow $$ 0^{+}_{g.s.}) for ^{10}Be with spurious center-of-mass motion treatment show good agreement with experimental data. The deformation properties of the 2^{+}_{1} and 2^{+}_{2} states for ^{10}Be and of the 2^{+}_{1} state for ^{12}Be are studied in terms of electric quadrupole moments, E2 transitions and single-particle occupations. The triaxial deformation of ^{10}Be is also discussed in terms of the B(E2; 2$$^+_{2}\rightarrow $$ 2^{+}_{1}) value. This work motivates a further extension of the MCSM application to ab initio FCI calculations [77]. Currently, the availability of the MCSM for no-core calculations has been tested extensively in light nuclei [54].

As a typical example, the behavior of the ground-state energies of ^{4}He (0^{+}) with respect to the number of basis states and to the energy variance in *N*_{shell}=2–5 are shown in Fig. 6. Figure 7 illustrates the comparisons of the energies for each state and model space between the MCSM and FCI methods. The FCI gives the exact energies in the fixed size of the mode space, while the MCSM gives approximated ones. Thus the comparisons between them show how well the MCSM works in no-core calculations. For this benchmark comparison, the JISP16 two-nucleon interaction is adopted and the Coulomb force is turned off. Isospin symmetry is assumed. The energies are evaluated for optimal harmonic oscillator frequencies where the calculated energies are minimized for each state and model space. Here, contributions from spurious center-of-mass motion are ignored for simplicity. In Fig. 7, the comparisons are made for the states; ^{4}He(0^{+}), ^{6}He(0^{+}), ^{6}Li(1^{+}), ^{7}Li(1/2^{−}, 3/2^{−}), ^{8}Be(0^{+}), ^{10}B(1^{+}, 3^{+}), and ^{12}C(0^{+}). The model space ranges from *N*_{shell}=2 through 5 for *A*≤6 (4 for *A*≥7). Note that the energies of ^{10}B(1^{+}, 3^{+}) and ^{12}C(0^{+}) in *N*_{shell}=4 are available only from the MCSM results. The *M*-scheme dimensions for these states (1.82×10^{10} for *M*=1 and 1.52×10^{10} for *M*=3 in ^{10}B and 5.87×10^{11} for *M*=0 in ^{12}C) are already marginal or exceed the current limitation in the FCI approach. The number of basis states is taken up to 100 in *N*_{shell}=2–4 and 50 in *N*_{shell}=5. In Fig. 7, the solid (dashed) lines indicate the MCSM (FCI) results. The shaded regions express the extrapolations in the MCSM, and the lower bound of the shaded region corresponds to the extrapolated energy. Furthermore, we also plot the NCFC results for the states of 4≤*A*≤8 as fully converged energies in the infinite model space. As seen in Fig. 7, the energies are consistent with each other where the FCI results are available to within ∼100 keV (∼500 keV) at most of the MCSM results with(out) the energy-variance extrapolation in the MCSM. The other observables besides the energies also give reasonable agreements between the MCSM and FCI results. Detailed comparisons of the MCSM, FCI, and NCFC methods can be found in Ref. [54].

By exploiting the recent development in the computation of the Hamiltonian matrix elements between non-orthogonal Slater determinants [80] and the technique of energy-variance extrapolation [22], the observables give good agreement between the MCSM and FCI results in *p*-shell nuclei. From the benchmark comparison, the no-core MCSM is now verified in its application to ab initio no-core calculations for light nuclei. Moreover, application of the no-core MCSM to heavier nuclei is expected in the near future.

### Analysis of intrinsic state

While ab initio approaches have been studied intensively in light nuclei, it is relatively difficult to study the cluster structure in an ab initio way. Among these approaches, Green’s function Monte Carlo method first provided the two-*α* structure of the ^{8}Be ground state illustratively [84]. This study has shown the possibility that the cluster structure can appear in ^{8}Be, without assuming any cluster structure in advance. Generalizing this result, it may be possible to treat cluster structure from a pure single particle picture. In this subsection, we show how to visualize the cluster state in the no-core MCSM calculation and by analyzing the calculations we discuss the appearance of the *α* cluster structure. It is also suitable for clarifying the relation between the shell-model and cluster pictures [85] from the shell-model point of view. This viewpoint has not been investigated very well yet. Recently, the density profile of lithium isotopes has been investigated by the NCFC [86]. The method has shown how to calculate the translationally-invariant density. In Li isotopes, shape distortion and cluster-like structure has been found. Thus, the study of cluster structure has become a realistic subject by using the shell-model calculation.

To extract the cluster structure from the no-core MCSM, we define the intrinsic state to visualize the cluster shape in the intrinsic framework, which is extracted from the angular-momentum-projected wave function. The wave function of the no-core MCSM, which is defined in Eq. (2.5), is represented as an angular-momentum projection of a linear combination of basis states such as

*I*is assumed to be zero and

*K*-quantum number and parity projections are omitted for simplicity. This linear combination of the unprojected basis states, |

*Φ*〉, cannot be considered an intrinsic state because the principal axis of a basis state, |

*ϕ*

_{i}〉, is not in the same direction as that of another basis state. Therefore we rotate each basis state so that it has a diagonalized quadrupole-moment;

*Q*

_{zz}>

*Q*

_{yy}>

*Q*

_{xx}and

*Q*

_{ij}=0,(

*i*≠

*j*), respectively, following the concept of Ref. [84]. As a result, these rotated basis states have a large overlap with each other and show a distinct principal axis toward the

*z*-axis. The intrinsic wave function |

*Φ*

^{intr}〉 is defined as

*R*(

*Ω*

_{n}) is a rotation operator with Euler’s angle

*Ω*

_{n}.

*Ω*

_{n}is determined so that the transformed basis state $$|\phi _n^R\rangle = R(\Omega _n) | \phi _n\rangle $$ has a diagonalized quadrupole-moment. The transformed coefficient $$D^R_n$$ (by

*R*(

*Ω*

_{n})) is derived by the relation in Ref. [17]. This state has exactly the same energy after the angular momentum projection. We calculate the one-body density of the intrinsic state such as

*r*

_{i}denotes the position of the

*i*th nucleon.

As an illustrative example, we show the ^{8}Be density in *N*_{shell}=4 and $$\hbar \omega =20$$ MeV with the JISP16 interaction for *J*=0^{+} states. The Coulomb interaction and the contamination of spurious center-of-mass motion are neglected for simplicity. We show the proton density (half of the total density) of |*Φ*〉 and the intrinsic-state density, *ρ*^{intr}, in Fig. 8.

The number of basis states is *N*_{b}=1,10, and 100 for the lower, middle, and upper rows, respectively. The energy is almost converged at *N*_{b}=100. Each density distribution is shown along the *yz* planes at *x*=0 fm and at *x*=1 fm. As shown in the *N*_{b}=1 results, clear deformation and a neck structure, to be called a dumbbell shape, appear. We can see that, as the number of basis states, *N*_{b}, increases, the density of |*Φ*〉 is much more vague and becomes ordinary prolate rather than dumbbell-like because of the mixture of different directions of the principal axes of the basis states. On the other hand, the intrinsic density has clearer dumbbell-like structure for each *N*_{b}. In addition, the density distribution of the intrinsic state is almost unchanged with respect to *N*_{b}. This result indicates the appearance of cluster structure in the no-core MCSM. We also check how the cluster shape differs between *N*_{shell}=3 and *N*_{shell}=4. We find that the neck of dumbbell shape is more enhanced in *N*_{shell}=4 than in *N*_{shell}=3. Since the weights of distribution for both sides of the principal axis are almost the same, this cluster can be considered as two *α* clusters. The stability of the *α* cluster is confirmed with respect to *N*_{b} and *N*_{shell}. This result is consistent with the result of Green’s function Monte Carlo method [84]. Using this method to draw the density, we can study the appearance of cluster structure directly, not only for *N*=*Z* nuclei but also for neutron-rich nuclei, in the ab initio approach. Research into exotic structure, including unstable nuclei in the *p*-shell region, is in progress.

## Application to neutron-rich Cr and Ni isotopes

In this section, we discuss the application of the MCSM to large-scale shell-model calculations using neutron-rich Cr and Ni isotopes as examples. We take the *pfg*_{9}*d*_{5} shell, which consists of the 0*f*1*p* shell, the 0*g*_{9/2} orbit, and the 1*d*_{5/2} orbit, as the model space. By using such a sufficiently large model space, we aim to give a unified description of medium–heavy nuclei, and to study the shell evolution [23, 24, 25, 26] and the magicity of *N*=28, 40, 50 microscopically.

### Ni isotopes and magicity of *N*=28,40,50

Nuclear shell structure evolves in neutron-rich nuclei and the magic numbers of unstable nuclei are different from those of stable nuclei. The large excitation energy of the 2^{+} yrast state and the small $$B(E2;0^+ \rightarrow 2^+)$$ value in ^{68}Ni (*Z*=28, *N*=40) might indicate that ^{68}Ni is a double-magic nucleus, although *N*=40 is a magic number of the harmonic oscillator, not a magic number of the nuclear shell model. On the other hand, the small excitation energies of the 2^{+} yrast state and the large $$B(E2;0^+ \rightarrow 2^+)$$ values in Cr (*Z*=24) isotopes of *N*∼40 suggest rather strong deformation. This change in the *N*=40 gap has been studied theoretically [87]. ^{78}Ni, which has *Z*=28 and *N*=50 doubly magic numbers, has also been investigated to discuss its magicity and the size of the *N*=50 gap [92].

In the *sd*-shell and the light *pf*-shell regions, we can describe the properties of stable nuclei in relatively small model spaces. However, we sometimes need a large model space to describe the properties of neutron-rich nuclei. In order to discuss neutron-rich Ni isotopes up to *N*=50, it is essential to include the effects of excitation across the *Z*=28 and *N*=50 gaps by adopting the *pfg*_{9}*d*_{5} model space. Concerning this model space, M. Honma et al. proposed the A3DA effective interaction (M. Honma et al., unpublished), which consists of the GXPF1A [88], JUN45 [50], and *G*-matrix effective interactions with phenomenological modifications. It has succeeded in describing neutron-rich Cr and Ni isotopes under a severe truncation of the model space utilizing the few-dimensional basis approximation [33, 34]. In this work, we use the new version of the MCSM method, which enables us to precisely evaluate the exact shell-model energy without any truncation and discuss the effective interaction.

### Effective interaction for the *pfg*_{9}*d*_{5} shell

In this section, we discuss the A3DA effective interaction (M. Honma et al., unpublished) and its improvement. The two-body matrix elements (TBMEs) of the A3DA interaction consist of three parts. The TBMEs of the *pf* shell are those of the GXPF1A interaction [88], which is successful for describing spectroscopic properties of light *pf*-shell nuclei. The TBMEs of the *f*_{5}*pg*_{9} shell related to the 0*g*_{9/2} orbit are those of the JUN45 interaction [50]. The GXPF1A and JUN45 interactions were determined by combining microscopically derived interactions (*G* matrix) with a minor empirical fit so as to reproduce experimental data. The other TBMEs are from the *G*-matrix effective interaction [89] (M. Hjorth-Jensen, private communication), which is calculated from the chiral N3LO interaction [57]. The Coulomb interaction is not included and the isospin symmetry is conserved. The *G* matrix is calculated for the *pfsdg* shell with ^{40}Ca as an inert core and the core-polarization correction is included perturbatively. The single-particle energies and the monopole interaction are adjusted to reproduce the GXPF1A and JUN45 predictions for the *pf* shell and *g*_{9/2} orbits, and the Woods–Saxon single-particle energies of stable semi-magic nuclei for the other part.

The original A3DA interaction failed to describe some nuclei around *N*∼40. We modify mainly single-particle energies and monopole components related to the 0*g*_{9/2} orbit by comparing the results of the calculations with experiments. These calculations are far beyond the current limitation of the conventional diagonalization method, and the MCSM method enables us to perform this comparison.

### MCSM results for neutron-rich Cr and Ni isotopes

We performed systematic calculations of the 0^{+} and 2^{+} yrast states of neutron-rich Cr and Ni even–even isotopes using the MCSM method and the modified A3DA interaction. We took 50 basis states for the MCSM with the refinement procedure discussed in Sect. 2.2. The energies were extrapolated by the energy-variance extrapolation method and the other values were not. The effective charges are taken as (*e*_{p},*e*_{n})=(1.5,0.5)*e*.

Figure 9 shows the 2^{+} excitation energies and the *B*(*E*2) transition probabilities of neutron-rich Cr isotopes. The MCSM results reproduce the experimental values well, while a modest overestimation remains. The Cr isotopes do not show any feature of *N*=40 magicity, while the characteristics of *N*=28 magicity can be seen, namely, a sudden increase of excitation energy and a slight decrease in the *B*(*E*2) value. On the neutron-rich side, the excitation energy decreases and the *B*(*E*2) value gradually increases as the neutron number increases, which implies gradual enhancement of the quadrupole deformation.

Figure 10 shows 2^{+} excitation energies of Ni (*Z*=28) even–even isotopes from ^{56}Ni to ^{78}Ni. The large 2^{+} excitation energy of ^{56}Ni (*N*=28) indicates *Z*=28, *N*=28 double magicity. The large value of the calculated 2^{+} excitation energy of ^{78}Ni (*N*=50) suggests *Z*=28, *N*=50 double magicity. The large 2^{+} excitation energy of ^{68}Ni (*N*=40) indicates *N*=40 magicity. The calculated values reproduce the experimental values well.

Figure 11 shows $$B(E2;0^+ \rightarrow 2^+)$$ for neutron-rich Ni isotopes. The small value of $$B(E2;0^+ \rightarrow 2^+)$$ at *N*=40 indicates *N*=40 magicity. Neither the theoretical nor the experimental value of $$B(E2;0^+ \rightarrow 2^+)$$ at *N*=28 is small, unlike that at *N*=40, and the theoretical $$B(E2;0^+ \rightarrow 2^+)$$ value at *N*=50 becomes large in comparison with those of neighboring nuclei. This suggests that, at *N*=28,50, magicity is broken to some extent for ^{56,78}Ni, respectively. Figure 12 shows the occupation number of the neutron *g*_{9/2} orbit. The occupation numbers of the 0^{+} and 2^{+} states are very close for Ni isotopes besides ^{68,78}Ni (*N*=40,50). The occupation numbers of the 2^{+} states of ^{68,78}Ni show a breakdown of the closed-shell structure. Figure 13 shows two-neutron separation energies. The calculated values of neutron-rich nuclei are smaller than experimental values. This means that the binding energies of neutron-rich nuclei are underestimated. The values of *S*_{2n} increase slightly by considering the Coulomb energy, but calculated values are still smaller than experimental values.

Figure 14 shows the total energy surface of ^{68}Ni provided by the *Q*-constrained Hartree–Fock calculation [95]. There are three minimum points for ^{68}Ni. Figure 14 also shows quadrupole deformations of the MCSM wave functions of the 0^{+}_{1} and 0^{+}_{2} states. The scattered circles correspond to the basis states in the MCSM wave function. The position of the circle indicates the quadrupole deformation of the basis state before projection. The area of the circle is proportional to the overlap probability of the projected basis and the resulting wave function. It is quite clear that the 0^{+}_{1} state of ^{68}Ni corresponds to a spherical shape and the 0^{+}_{2} state corresponds to an oblate shape. Spherical and oblate components are mixed to some extent, but the components of the prolate minimum hardly mix.

Furthermore, we consider the magicity and the energy gaps for Ni isotopes by using the effective single particle energies (ESPEs) [93]. Figure 15 shows the ESPEs of the neutron orbits for Ni isotopes. The *f*_{7/2}–*p*_{3/2} gap at *N*=28 is 7.1 MeV and gives the magicity to ^{56}Ni. The *g*_{9/2}–*d*_{5/2} gap at *N*=50 is 4.2 MeV and gives the magicity to ^{78}Ni. This is partly due to the additional lowering of the *g*_{9/2} orbit caused by the pairing correlation between two neutrons in the *g*_{9/2} orbit, and also due to the effect of the two-neutron repulsive monopole interaction originating in the three-nucleon force, like in exotic oxygen isotopes [26]. The *p*_{1/2}–*g*_{9/2} gap at *N*=40 is 2.6 MeV, which is smaller than the *N*=28,50 gaps. Figure 16 shows the ESPEs of the neutron orbits for *N*=40 isotones. As the proton number of *f*_{7/2} increases from *Z*=20 to *Z*=28, the ESPE of *f*_{5/2} decreases and the *N*=40 gap becomes larger. Because of this evolution of the *N*=40 gap, the properties of *N*∼40 nuclei depend on the proton number.

In Fig. 17, the ESPEs of the proton orbits for Ni isotopes and for *N*=40 isotones are shown. In the former, rapid lowering of the *f*_{5/2} orbit from *N*=40 to 50 is clearly seen, as suggested in [24, 25], while narrowing of the *Z*=28 gap is also visible there. Such changes are partly responsible for the origins of the structure evolution in these Ni isotopes.

## Summary and future perspectives

We have developed a new generation of the MCSM by introducing the conjugate gradient method and energy-variance extrapolation, which greatly enhance the applicability of the MCSM. We describe two major applications of this framework: ab initio shell-model calculations and conventional shell-model calculations assuming an inert core. In the former, we have compared the MCSM results with exact FCI calculations to demonstrate the validity of the MCSM framework and its feasibility beyond the limit of the FCI in Sect. 3. In addition, we have proposed a novel method to discuss the intrinsic structure and demonstrated that the cluster structure appears in shell-model-type calculations based on the harmonic-oscillator-basis wave function. In the latter, we discussed in Sect. 4 that the MCSM enables us to perform shell-model calculations of neutron-rich Cr and Ni isotopes in the *pfg*_{9}*d*_{5} model space in which isospin symmetry is conserved. We proposed a “modified A3DA” interaction which reproduces the low-lying spectra of neutron-rich Cr and Ni isotopes and guides us towards a unified description including ^{56}Ni, ^{68}Ni, and ^{78}Ni, with magic numbers 28, 40, and 50, respectively. The prediction of ^{78}Ni is especially interesting to see the evolution of shell structure. On the other hand, Cr isotopes do not show any feature of *N*=40 magicity and the collectivity enhances as the neutron number increases. The MCSM and newly proposed effective interaction are expected to provide us with a unified description of *pf*-shell nuclei.

The current status of the computer-code development was also reported in Sect. 2.5. At the present stage, we have obtained good parallel scalability of our code up to 10^{5} CPU cores via early access to the K computer at RIKEN AICS [94] as measured by benchmark testing. However, such good scalability is not always obtained and further development is in progress. This activity is strongly promoted as part of the activities of HPCI Strategic Programs for Innovative Research (SPIRE) Field 5, “The origin of matter and the universe”.

By utilizing both the developed code and the K computer, we promote further large-scale shell-model calculations as part of the SPIRE activities. We plan to perform a systematic study with ab initio calculations of light nuclei in *N*_{shell}=5 and some states in *N*_{shell}=6. Regarding medium–heavy nuclei, because it is difficult to cover the whole region of the nuclear chart, we will choose some interesting nuclides as the subjects of our investigation, and will perform shell-model calculations of these nuclides with the two-major-shell model space. For example, the shell-model calculations of ^{130}Te, ^{128}Te, and ^{150}Nd are extremely interesting to study double beta decay and the nuclear matrix element of neutrinoless decay. We also continue to study the systematic calculations of neutron-rich *pf*-shell nuclei to discuss the shell-evolution phenomenon.

## Acknowledgments

We acknowledge Professors J. P. Vary, P. Maris and Dr L. Liu for our collaboration concerning ab initio shell-model calculations. This work has been supported by Grants-in-Aids for Scientific Research (23244049), for Scientific Research on Innovative Areas (20105003), and for Young Scientists (20740127) from JSPS, the SPIRE Field 5 from MEXT, and the CNS-RIKEN joint project for large-scale nuclear structure calculations. The numerical calculations were performed mainly on the T2K Open Supercomputers at the University of Tokyo and Tsukuba University. The exact conventional shell-model calculations were performed by the code MSHELL64 (T. Mizusaki, N. Shimizu, Y. Utsuno, and M. Honma, code MSHELL64, unpublished).

### Appendix: Numeration with projected Slater determinants

In this appendix, we show some equations that are required to perform the calculation discussed in Sect. 2.

At the beginning, we define a deformed Slater determinant,

*N*

_{sp}×

*N*

_{f}matrix

*D*with the normalization condition

*D*

^{†}

*D*=1.

*N*

_{f}and

*N*

_{sp}are the numbers of fermions and single-particle states, respectively. The |−〉 denotes an inert core in conventional shell-model calculations or the vacuum in ab initio shell-model calculations. Because we do not mix the proton and neutron space in practical calculations, the wave function is written as a product of proton and neutron Slater determinants, namely, |

*ϕ*〉=|

*ϕ*

_{proton}〉⊗|

*ϕ*

_{neutron}〉. For simplicity, we do not write this isospin degree of freedom explicitly. One can easily reproduce the equations representing the explicit proton–neutron degree of freedom by taking

*D*of the proton–neutron sector as zero, such as

*D*

^{π}and

*D*

^{ν}represent Slater determinants of protons and neutrons, respectively.

The angular-momentum parity projector $$P^{I\pi }_{MK}$$ in Eq. (2.6) is performed by discretizing the integral concerning the Euler angles such as

*λ*denotes an index of the mesh point of discretization (here, a set of Euler’s angle

*Ω*=(

*α*,

*β*,

*γ*) and parity variable

*π*

^{λ}=±1). In this paper, the parity projection is described by the summation of 2 mesh points such as $$P^\pi = \frac {1 + \pi \Pi }{2} = \sum _{\lambda =1}^2 \pi ^{(\lambda)}\Pi ^{(\lambda)}$$ with $$\pi ^{(1)} = \frac {1}{2}$$, $$\pi ^{(2)} = \frac {\pi }{2}$$,

*Π*

^{(1)}=1, and

*Π*

^{(2)}=

*Π*, with

*Π*being the parity-conversion operator.

*W*

^{(λ)}is a weight of the mesh point

*λ*, and

*R*

^{(λ)}is a product of the rotation and parity-conversion operators such as

*R*

^{λ}does not change the form of a Slater determinant, i.e.,

*D*

^{n(λ)}represents the single Slater determinant $$|\phi ^{(\lambda)}_n \rangle $$, thanks to the Baker–Hausdorff theorem [17].

The norm matrix and Hamiltonian matrix spanned by *N* Slater determinants are written as

*f*

_{nK}in Eq. (2.5), is determined by solving the generalized eigenvalue problem

*Ψ*|

*Ψ*〉=1. The lowest eigenvalue of $$\mathcal {E}$$ is taken as

*E*

_{N}if you would like to obtain the yrast state.

By combining Eqs. (A3) and (A6), the norm matrix is calculated as

In the same way, the Hamiltonian matrix is obtained as

*ρ*

^{(λ)}, and the self-consistent field,

*Γ*

^{(λ)}, [95] are defined as

*m*,

*n*of

*ρ*

^{(λ)}and

*Γ*

^{(λ)}are also omitted.

The most-time-consuming part is the calculation of the $$\Gamma ^{(\lambda)}_{ik}$$, which can be rewritten following the idea of Ref. [51],

*a*=(

*i*,

*k*), and

*b*=(

*j*,

*l*). Because $$\overline {v}_{ab}$$ is a block-antidiagonal form owing to the symmetry of the Hamiltonian, Eq. (A14) is calculated as the products of the dense block matrices and the dense matrices in terms of the indices

*a*,

*b*,

*λ*, efficiently avoiding trivial zero matrix elements of

*v*

_{ab}. This method is referred to as the

*matrix–matrix method*in Ref. [51]. This

*matrix–matrix method*enables us to use a CPU utilizing the BLAS level 3 library quite efficiently, and the performance reaches 70∼80% of the theoretical peak performance [51].

This efficient computation of *Γ*^{(λ)} is also useful for the evaluation of the energy gradient, which is essential for the conjugate gradient method. The energy gradient of the Slater-determinant coefficients is written as

To evaluate the energy variance $$\langle \Delta H^2\rangle _N = \langle H^2\rangle _N - E_N^2$$, the expectation value of *H*^{2} with the wave function in Eq. (2.5) is written as

*i*,

*j*),(

*k*,

*l*),(

*α*,

*β*), and (

*γ*,

*δ*) with

*a*,

*b*,

*c*, and

*d*, respectively, this term is efficiently calculated as the products of the matrices, namely, $$\sum _{abcd} v_{ab} \Theta _{bc} v_{cd} \Lambda _{da}$$. Note that

*v*

_{ab}has a block-diagonal form, which again enables us to use the BLAS level 3 library, avoiding trivial zero matrix elements.

Another formulation to compute the expectation values in projected Slater determinants can be found in Refs. [37, 96], in which the two-body interaction is decomposed into the sum of the squares of the one-body operators.