## 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 40Ca 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 1s-0d (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 1p–0f (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 NpCnp×NnCnn without symmetry consideration, where Np (Nn) is the number of proton (neutron) single-particle states taken and np (nn) is the number of protons (neutrons) activated. It roughly increases exponentially with Np (or Nn) and also with np (or nn). Hence, without even considering the ab initio shell model with a huge number of Np and Nn, 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 (Np or Nn) increases for heavier nuclei. For instance, the dimension needed for A∼80 N=Z nuclei is estimated to be ∼1027 without any symmetry consideration [13], when the 1p–0f–2s–1d–0g 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 Jz as are usually taken (i.e., M-scheme calculation), it is still far from the current computational limit of 1010−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 (Nsp)α with α∼3–4 in the case of Nsp=Np=Nn, being much weaker than the exponential increase. This advantage enabled us to perform the full pf-shell calculation in 56Ni (with ∼109M-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:

(2.1)
$$\label{eq2.1} H = H^{(1)} + H^{(2)} = \sum_{i} t_{i} c^\dagger_i c_i + \sum_{i<j, k<l} v_{ijkl} c^\dagger_i c^\dagger_j c_l c_k,$$
where $$c^\dagger _i$$ denotes a creation operator of single particle state i. H(1) is a one-body Hamiltonian with single-particle energies ti, 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

(2.2)
$$\label{eq2.2} H = T -T_{\mathrm{CM}} + V = \sum_{ij} t_{ij} c^\dagger_i c_j + \sum_{i<j, k<l} v_{ijkl} c^\dagger_i c^\dagger_j c_l c_k ,$$
where T and TCM are the total kinetic energy and the kinetic energy of the center-of-mass motion, respectively. Note that TCM 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

(2.3)
$$\label{eq2.3} H^{\prime} = H + \beta_{\mathrm{cm}} H_{\rm cm}$$
with HCM being defined as
(2.4)
$$\label{eq2.4} H_{\mathrm{cm}} = \frac{\textbf{P}^{\textbf{2}}}{2AM} + \frac{1}{2} MA\omega^2 \textbf{R}^2 - \frac{3}{2} \hbar \omega,$$
where R and P are the position and momentum of the center of mass. By taking βCM large enough, 〈HCM〉 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,

(2.5)
$$\label{eq2.5} |\Psi_{N_b} \rangle = \sum_{n=1}^{N_b} \sum_{K=-I}^I f^{(N_b)}_{n,K} P^{I \pi}_{MK}| \phi_n \rangle,$$
where Nb is the number of the Slater-determinant basis states. The $$P^{I\pi }_{MK}$$ operator is the angular-momentum and parity projector defined as
(2.6)
$$\label{eq2.6} P^{I\pi}_{MK} = \frac{1 + \pi \Pi}{2} \frac{2I+1}{8\pi^2} \int d\Omega \ {D^{I}_{MK}}^{\ast}(\Omega) e^{i\alpha J_z} e^{i\beta J_y} e^{i\gamma J_z},$$
where Ω≡(α,β,γ) 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 −IKI and 1≤nNb. This diagonalization also determines the energy, ENb≡〈ΨNb|H|ΨNb〉, as a function of Nb. Note that the dimension of the subspace is (2I+1)Nb, not Nb. The Slater determinant basis, D(n), is given by variational calculation to minimize ENb=n. We increase Nb until ENb 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 ENb, 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 (Nb≃30) compared to that in the MCSM (Nb≃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:

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

2. Sequential optimization for each basis state (SCG),

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

4. 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 Ei 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, Nb, 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 ENb, since D(1) is determined to minimize E1. Then, we first fix the number of basis states, Nm, and restart to minimize the energy ENm 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 Nm at the beginning, and determine all coefficients D(Nb) with 1≤NbNm by minimizing the ENm at once. In other words, we perform a variational calculation with a set of all variational parameters, D(Nb), simultaneously to minimize the ENm. 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, Nm, 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 72Ge in the f5pg9 shell, which consists of the 0f5/2, 1p3/2, 1p1/2, and 0g9/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 ENb=〈ΨNb|H|ΨNb〉 is plotted against the number of basis states, Nb. 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 Nm(Nm+1)/2.

Fig. 1.

Convergence patterns of the ground-state energy of 72Ge in the f5pg9 shell. The circles, open triangles, and squares denote the results of the original MCSM, SCG, and FCG with Nm=32, respectively.

Fig. 1.

Convergence patterns of the ground-state energy of 72Ge in the f5pg9 shell. The circles, open triangles, and squares denote the results of the original MCSM, SCG, and FCG with Nm=32, respectively.

In the case of the FCG, the number of basis states of the variational wave function is fixed to Nm=32 and we plot energy expectation values in the subspace spanned by the Nb basis states. In the FCG, all basis states, D(1),D(2),…, and D(32), are optimized to minimize E32, while in the SCG wave function D(1) is optimized to minimize E1. Thus, the calculated E1 of the FCG is much higher than the E1 of the SCG and the E1 of the original MCSM. However, E32 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 Nm in the FCG, not as a function of Nb, unlike the case of the SCG.

While the FCG yields the best energy with fixed Nm, 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≤NbNm, 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 ENb=〈Ψ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 Nb 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

(2.7)
$$\label{eq2.7} E = c_0 + c_1 \langle \Delta H^2 \rangle + c_2 (\langle \Delta H^2 \rangle)^2,$$
where the coefficients c0, c1, and c2 are determined by least-squares fiting. By extrapolating the fitted curve into the y-intercept we obtain the extrapolated energy, namely, c0.

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≤NbNm, where Nm is the maximum of Nb. For each Nb, we evaluate the energy ENb and energy variance 〈ΔH2Nb. Here, we demonstrate how the extrapolation method works with 56Ni 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 56Ni. In this calculation, we take the K=I state only in the angular-momentum projector, namely $$f_{n,K}^{(N_b)}=0$$ if KI in Eq. (2.5), for simplicity. As Nb 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 〈ΔH2Nb=100=0.89 MeV2, which is smaller and gives a better approximation than the result of the original MCSM, 〈ΔH2Nb=150=1.05 MeV2 [22], mainly thanks to the introduction of the CG method.

Fig. 2.

Energy-variance plot for 56Ni (a) without reordering, (b) with reordering. The open circles on the y-axis denote the exact shell-model energies. See text for details.

Fig. 2.

Energy-variance plot for 56Ni (a) without reordering, (b) with reordering. The open circles on the y-axis denote the exact shell-model energies. See text for details.

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 〈ΔH2〉≃4 MeV2. As a result, the simple extrapolation method apparently fails. A straightforward solution to this failure is to increase Nb 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 56Ni 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. 56Ni 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 72Ge 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 72Ge 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 Nm=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 Nm=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 Nm=32.

Fig. 3.

Energy vs. energy-variance plot of the ground state energy of 72Ge in the f5pg9 shell. The wave functions are provided by the FCG with Nm=8 (triangles), 16 (squares), 24 (diamonds), and 32 (circles).

Fig. 3.

Energy vs. energy-variance plot of the ground state energy of 72Ge in the f5pg9 shell. The wave functions are provided by the FCG with Nm=8 (triangles), 16 (squares), 24 (diamonds), and 32 (circles).

### 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 Nmesh=262×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 Nbunch being e.g. 30, we still have the NbNmesh/Nbunch≃721Nb elements to be computed in parallel.

Figure 4 shows the parallel efficiency of the benchmark calculation of the ground state of 64Ge as an example. The model space consists of the pf shell and g9/2 orbit; the PFG9B3 effective interaction is used (M. Honma et al., unpublished). Its M-scheme dimension reaches 1.7×1014, 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].

Fig. 4.

(a) Speedup of the parallel computation of the SCG process in units of computation time using 16 CPU cores. The squares, open circles, triangles, and filled circles represent the inverses of computation times of the variational process to obtain the 1st, 4th, 16th, and 32nd basis states respectively. The solid line shows ideal scaling to guide the eye. (b) Speedup of the parallel computation of the energy variance. The squares, triangles, and circles represent the inverses of computation times of 1, 4, and 8 basis states, respectively.

Fig. 4.

(a) Speedup of the parallel computation of the SCG process in units of computation time using 16 CPU cores. The squares, open circles, triangles, and filled circles represent the inverses of computation times of the variational process to obtain the 1st, 4th, 16th, and 32nd basis states respectively. The solid line shows ideal scaling to guide the eye. (b) Speedup of the parallel computation of the energy variance. The squares, triangles, and circles represent the inverses of computation times of 1, 4, and 8 basis states, respectively.

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, vijkl, 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 Nb(Nb+1)×Nmesh 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 Nb≥4, which is shown in Fig. 4(b).

In practice, it took 807 seconds in total to obtain an SCG wave function of 64Ge 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 1010 as shown in Fig. 5. So far the largest calculations have been done in 14N with $$N_{\max } = 8$$, which results in the M-scheme many-body matrix dimensions being ∼109 and associated non-vanishing three-nucleon force matrix elements being ∼4×1013 [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.

Fig. 5.

M-scheme dimensions as functions of basis-space size, Nshell.

Fig. 5.

M-scheme dimensions as functions of basis-space size, Nshell.

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 Nshell 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 10Be and 12Be have been investigated. The excitation energies of the first and second 2+ states and the B(E2; 2$$^+_{1}\rightarrow$$ 0+g.s.) for 10Be 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 10Be and of the 2+1 state for 12Be are studied in terms of electric quadrupole moments, E2 transitions and single-particle occupations. The triaxial deformation of 10Be 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 4He (0+) with respect to the number of basis states and to the energy variance in Nshell=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; 4He(0+), 6He(0+), 6Li(1+), 7Li(1/2, 3/2), 8Be(0+), 10B(1+, 3+), and 12C(0+). The model space ranges from Nshell=2 through 5 for A≤6 (4 for A≥7). Note that the energies of 10B(1+, 3+) and 12C(0+) in Nshell=4 are available only from the MCSM results. The M-scheme dimensions for these states (1.82×1010 for M=1 and 1.52×1010 for M=3 in 10B and 5.87×1011 for M=0 in 12C) are already marginal or exceed the current limitation in the FCI approach. The number of basis states is taken up to 100 in Nshell=2–4 and 50 in Nshell=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].

Fig. 6.

4He ground-state energies as a function of the number of basis states (left) and energy variance (right). The red, green, blue, and purple solid symbols (horizontal dashed lines in the left figure and open symbols at the zero energy variance in the right figure) are the MCSM (FCI) results in Nshell=2, 3, 4, and 5, respectively. The harmonic oscillator energies are taken at optimal values for each state and model space. The Coulomb interaction and the spurious center-of-mass motion effect are not considered. Isospin symmetry is assumed.

Fig. 6.

4He ground-state energies as a function of the number of basis states (left) and energy variance (right). The red, green, blue, and purple solid symbols (horizontal dashed lines in the left figure and open symbols at the zero energy variance in the right figure) are the MCSM (FCI) results in Nshell=2, 3, 4, and 5, respectively. The harmonic oscillator energies are taken at optimal values for each state and model space. The Coulomb interaction and the spurious center-of-mass motion effect are not considered. Isospin symmetry is assumed.

Fig. 7.

Comparisons of the energies between the MCSM and FCI along with the fully converged NCFC results where available [54]. The MCSM (FCI) results are shown as solid (dashed) lines that nearly coincide where both are available. The extrapolated MCSM results are illustrated by bands. From top to bottom, the truncation of the model space is Nshell=2 (red), 3 (green), 4 (blue), and 5 (purple). Note that the MCSM results are extrapolated by the energy variance with second-order polynomials [22]. Also note that all of the results for 10B and 12C at Nshell=4 were obtained only with MCSM.

Fig. 7.

Comparisons of the energies between the MCSM and FCI along with the fully converged NCFC results where available [54]. The MCSM (FCI) results are shown as solid (dashed) lines that nearly coincide where both are available. The extrapolated MCSM results are illustrated by bands. From top to bottom, the truncation of the model space is Nshell=2 (red), 3 (green), 4 (blue), and 5 (purple). Note that the MCSM results are extrapolated by the energy variance with second-order polynomials [22]. Also note that all of the results for 10B and 12C at Nshell=4 were obtained only with MCSM.

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 8Be ground state illustratively [84]. This study has shown the possibility that the cluster structure can appear in 8Be, 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

(3.1)
$$|\Psi \rangle = P^I | \Phi\rangle | \Phi\rangle = \sum_n f_n |\phi_n\rangle,$$
where the total 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; Qzz>Qyy>Qxx and Qij=0,(ij), 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
(3.2)
$$\label{eq3.1} |\Phi^{\rm intr} \rangle \equiv \sum_{n} f_{n} R(\Omega_n)|\phi_n \rangle = \sum_{n} f_{n}|\phi^R_n\rangle,$$
where 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
(3.3)
$$\label{eq3.2} \rho^{\mathrm{intr}} (r) = \langle \Phi^{\rm intr} |\sum_i \delta(r-r_i) | \Phi^{\rm intr} \rangle,$$
where ri denotes the position of the ith nucleon.

As an illustrative example, we show the 8Be density in Nshell=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.

Fig. 8.

8Be proton density for |Φ〉 (left panels) and intrinsic (right panels) states for various Nb and sliced along the yz plane. The number of basis states is Nb=1,10, and 100 for the lower, middle, and upper figures, respectively. The slice along the yz plane is the x=0 fm plane (left) or x=1 fm plane (right) for each panel. The size of each box is 8 fm × 8 fm.

Fig. 8.

8Be proton density for |Φ〉 (left panels) and intrinsic (right panels) states for various Nb and sliced along the yz plane. The number of basis states is Nb=1,10, and 100 for the lower, middle, and upper figures, respectively. The slice along the yz plane is the x=0 fm plane (left) or x=1 fm plane (right) for each panel. The size of each box is 8 fm × 8 fm.

The number of basis states is Nb=1,10, and 100 for the lower, middle, and upper rows, respectively. The energy is almost converged at Nb=100. Each density distribution is shown along the yz planes at x=0 fm and at x=1 fm. As shown in the Nb=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, Nb, 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 Nb. In addition, the density distribution of the intrinsic state is almost unchanged with respect to Nb. This result indicates the appearance of cluster structure in the no-core MCSM. We also check how the cluster shape differs between Nshell=3 and Nshell=4. We find that the neck of dumbbell shape is more enhanced in Nshell=4 than in Nshell=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 Nb and Nshell. 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 pfg9d5 shell, which consists of the 0f1p shell, the 0g9/2 orbit, and the 1d5/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 68Ni (Z=28, N=40) might indicate that 68Ni 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]. 78Ni, 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 pfg9d5 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 pfg9d5 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 f5pg9 shell related to the 0g9/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 40Ca 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 g9/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 0g9/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 (ep,en)=(1.5,0.5)e.

Figure 9 shows the 2+ excitation energies and the B(E2) 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(E2) value. On the neutron-rich side, the excitation energy decreases and the B(E2) value gradually increases as the neutron number increases, which implies gradual enhancement of the quadrupole deformation.

Fig. 9.

The excitation energies of 2+1 states (left) and $$B(E2;0^+_1 \rightarrow 2^+_1)$$ values (right) obtained by the MCSM for Cr isotopes. Experimental data are taken from Refs. [90, 91].

Fig. 9.

The excitation energies of 2+1 states (left) and $$B(E2;0^+_1 \rightarrow 2^+_1)$$ values (right) obtained by the MCSM for Cr isotopes. Experimental data are taken from Refs. [90, 91].

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

Fig. 10.

2+ excitation energies for Ni isotopes. Experimental data are taken from Ref. [90].

Fig. 10.

2+ excitation energies for Ni isotopes. Experimental data are taken from Ref. [90].

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,78Ni, respectively. Figure 12 shows the occupation number of the neutron g9/2 orbit. The occupation numbers of the 0+ and 2+ states are very close for Ni isotopes besides 68,78Ni (N=40,50). The occupation numbers of the 2+ states of 68,78Ni 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 S2n increase slightly by considering the Coulomb energy, but calculated values are still smaller than experimental values.

Fig. 11.

$$B(E2;0^+ \rightarrow 2^+)$$ values for Ni isotopes. Experimental data are taken from Ref. [91].

Fig. 11.

$$B(E2;0^+ \rightarrow 2^+)$$ values for Ni isotopes. Experimental data are taken from Ref. [91].

Fig. 12.

Occupation numbers of the neutron g9/2 orbit for Ni isotopes.

Fig. 12.

Occupation numbers of the neutron g9/2 orbit for Ni isotopes.

Fig. 13.

Two-neutron separation energies S2n for Ni isotopes. Experimental data are taken from Ref. [90].

Fig. 13.

Two-neutron separation energies S2n for Ni isotopes. Experimental data are taken from Ref. [90].

Figure 14 shows the total energy surface of 68Ni provided by the Q-constrained Hartree–Fock calculation [95]. There are three minimum points for 68Ni. 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 68Ni 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.

Fig. 14.

Total energy surface of the 0+1 (left) and 0+2 (right) states of 68Ni. The positions of the red circles represent quadrupole deformations of the MCSM basis states before projection. The areas of those circles represent the overlap probabilities of the basis states and the resulting wave function.

Fig. 14.

Total energy surface of the 0+1 (left) and 0+2 (right) states of 68Ni. The positions of the red circles represent quadrupole deformations of the MCSM basis states before projection. The areas of those circles represent the overlap probabilities of the basis states and the resulting wave function.

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 f7/2p3/2 gap at N=28 is 7.1 MeV and gives the magicity to 56Ni. The g9/2d5/2 gap at N=50 is 4.2 MeV and gives the magicity to 78Ni. This is partly due to the additional lowering of the g9/2 orbit caused by the pairing correlation between two neutrons in the g9/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 p1/2g9/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 f7/2 increases from Z=20 to Z=28, the ESPE of f5/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.

Fig. 15.

ESPEs of the neutron orbits for Ni (Z=28) isotopes.

Fig. 15.

ESPEs of the neutron orbits for Ni (Z=28) isotopes.

Fig. 16.

ESPEs of the neutron orbits for N=40 isotones.

Fig. 16.

ESPEs of the neutron orbits for N=40 isotones.

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 f5/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.

Fig. 17.

ESPEs of the proton orbits for Ni isotopes (left) and for N=40 isotones (right).

Fig. 17.

ESPEs of the proton orbits for Ni isotopes (left) and for N=40 isotones (right).

## 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 pfg9d5 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 56Ni, 68Ni, and 78Ni, with magic numbers 28, 40, and 50, respectively. The prediction of 78Ni 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 105 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 Nshell=5 and some states in Nshell=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 130Te, 128Te, and 150Nd 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,

(A.1)
$$|\phi \rangle = \prod_{k=1}^{N_{\mathrm{f}}} \left(\sum_{l=1}^{N_{\rm sp}} D_{lk} c^\dagger_l \right) | - \rangle , \label{eqA.1}$$
which is parametrized by the complex Nsp×Nf matrix D with the normalization condition DD=1. Nf and Nsp 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
(A.2)
$$\label{eqA.2} D = \left(\begin{matrix} D^\pi & 0 \\ 0 & D^\nu \end{matrix}\right),$$
where 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

(A.3)
$$\label{eqA.3} P^{I\pi}_{MK} = \sum_\lambda W^{I\pi(\lambda)}_{MK} R^{(\lambda)},$$
where the λ 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
(A.4)
\begin{align} W^{I\pi(\lambda)}_{MK} &= \frac{2I+1}{8\pi^2} D^{I*}_{MK}(\alpha_\lambda, \beta_\lambda, \gamma_\lambda) \pi^{(\lambda)}, \notag \\ R^{(\lambda)} &= e^{i\alpha_\lambda J_z} e^{i\beta_\lambda J_y} e^{i\gamma_\lambda J_z} \Pi^{(\lambda)}. \label{eqA.4} \end{align}
Note that the operator Rλ does not change the form of a Slater determinant, i.e.,
(A.5)
$$\label{eqA.5} |\phi^{(\lambda)}_n \rangle = R^{(\lambda)} | \phi_n \rangle,$$
where a matrix Dn(λ) 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

(A.6)
$$\mathcal{N}_{mM,nK} = \langle \phi_m| P^{I\pi}_{MK} | \phi_n \rangle \label{eqA.6}$$

(A.7)
$$\mathcal{H}_{mM,nK} = \langle \phi_m| H P^{I\pi}_{MK} | \phi_n \rangle. \label{eqA.7}$$
The coefficient, fnK in Eq. (2.5), is determined by solving the generalized eigenvalue problem
(A.8)
$$\label{eqA.8} \sum_{nK} \mathcal{H}_{mM,nK}f_{nK} = \mathcal{E} \sum_{nK} \mathcal{N}_{mM,nK} f_{nK},$$
and the normalization condition 〈Ψ|Ψ〉=1. The lowest eigenvalue of $$\mathcal {E}$$ is taken as EN if you would like to obtain the yrast state.

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

(A.9)
$$\label{eqA.9} \mathcal{N}_{mM,nK} = \sum_\lambda W^{I\pi(\lambda)}_{MK} \langle \phi_m | R^{(\lambda)} | \phi_n \rangle = \sum_\lambda W^{I\pi(\lambda)}_{MK} \langle \phi_m | \phi_n^{(\lambda)} \rangle ,$$
with
(A.10)
$$\label{eqA.10} \langle \phi_m | \phi_n^{(\lambda)} \rangle = {\mathrm{det}}(D^{m\dagger} D^{n(\lambda)}).$$

In the same way, the Hamiltonian matrix is obtained as

(A.11)
\begin{align} \mathcal{H}_{mM,nK} &= \sum_\lambda W^{I\pi(\lambda)}_{MK} \langle \phi_m | H R^{(\lambda)} | \phi_n \rangle \notag\\ &= \sum_\lambda W^{I\pi(\lambda)}_{MK} \langle \phi_m| \phi_n^{(\lambda)} \rangle {\rm Tr} \left(\rho^{(\lambda)} \left(t + \frac{1}{2} \Gamma^{(\lambda)}\right) \right), \label{eqA.11} \end{align}
where the generalized density matrix, ρ(λ), and the self-consistent field, Γ(λ), [95] are defined as
(A.12)
$$\rho^{(\lambda)}_{ij} = \frac{\langle \phi_m| c^\dagger_j c_i| \phi^{(\lambda)}_n \rangle} {\langle \phi_m| \phi_n^{(\lambda)} \rangle} = (D^{n(\lambda)}(D^{m\dagger} D^{n(\lambda)})^{-1} D^{m\dagger})_{ij}\label{eqA.12}$$

(A.13)
$$\Gamma^{(\lambda)}_{ik} = \sum_{jl} \overline{v}_{ijkl} \rho^{(\lambda)}_{lj} \label{eqA.13}$$
with $$\overline {v}_{ijkl} = v_{ijkl}-v_{ijlk}$$. The trivial summations and their indices for the matrix products are omitted for readability. The indices 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.14)
$$\label{eqA.14} \Gamma^{(\lambda)}_{a} = \sum_{b} \overline{v}_{ab} \rho^{(\lambda)}_{b},$$
where 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 vab. 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

x1-18024r15
\begin{align} \frac{\partial E_N}{\partial D^{m*}} &= (1-D^mD^{m\dagger}) \sum_{M,n,K,\lambda} f_{mM}^* f_{nK} W^{I\pi(\lambda)}_{MK} \langle \phi_m |\phi^{(\lambda)}_n \rangle \notag \\ & \quad \times \left((1 - \rho^{(\lambda)}) (t + \Gamma^{(\lambda)}) + \left({\mathrm{Tr}} \left(\left(t + \frac{1}{2} \Gamma^{(\lambda)}\right) \rho^{(\lambda)}\right) - E_N\right) \right) \rho^{(\lambda)} D^m. \label{eqA.15} \end{align}

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

x1-18026rA16
$$\langle \Psi_N | H^2 |\Psi_N \rangle = \sum_{m,M,n,K,\lambda} f_{mM}^* f_{nK} W^{I\pi(\lambda)}_{MK} \langle \phi_m| H^2 | \phi_n^{(\lambda)}\rangle. \label{eqA.16}$$
From Ref. [22], the matrix element of the Hamiltonian squared is computed such as
(A.17)
\begin{align} \frac{\langle \phi_m | H^2 | \phi^{(\lambda)}_n \rangle}{\langle \phi_m | \phi^{(\lambda)}_n \rangle} & = \sum_{i<j, k<l, \alpha<\beta, \gamma<\delta} v_{ijkl} \Theta^{(\lambda)}_{kl\alpha\beta} v_{\alpha\beta\gamma\delta} \Lambda^{(\lambda)}_{\gamma\delta ij} \notag\\ & \quad + {\mathrm{Tr}}((t+\Gamma^{(\lambda)}) (1-\rho^{(\lambda)})(t+\Gamma^{(\lambda)}) \rho^{(\lambda)}) + \left({\rm Tr}\left(\rho^{(\lambda)}\left(t+\frac{1}{2} \Gamma^{(\lambda)}\right)\right)\right)^2 \label{eqA.17} \end{align}

(A.18)
$$\Lambda^{(\lambda)}_{ijkl} = \rho^{(\lambda)}_{i k} \rho^{(\lambda)}_{j l} - \rho^{(\lambda)}_{i l} \rho^{(\lambda)}_{j k} \label{eqA.18}$$

(A.19)
$$\Theta^{(\lambda)}_{ijkl} = (1- \rho^{(\lambda)})_{i k}(1-\rho^{(\lambda)})_{j l} - (1- \rho^{(\lambda)})_{i l}(1-\rho^{(\lambda)})_{j k}. \label{eqA.19}$$
The most-time-consuming part in the evaluation of the energy variance is the calculation of the first term of the right-hand side of Eq. (A17). By substituting (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 vab 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.

## References

1
H.
, et al.  .
Phys. Rev. C
,
2001
, vol.
64
pg.
044001

()
2
Mayer
M. G.
Phys. Rev.
,
1949
, vol.
75
pg.
1969

()
3
Hazel
O.
Jensen
J. H. D.
Suess
H. E.
Phys. Rev.
,
1949
, vol.
75
pg.
1766

()
4
Brown
B. A.
Wildenthal
B. H.
Ann. Rev. Nucl. Part. Sci.
,
1988
, vol.
38
pg.
29

()
5
Brown
B. A.
Richter
W. A.
Phys. Rev. C
,
2006
, vol.
74
pg.
034315

()
6
Poves
A.
Sánchez-Solano
J.
Caurier
E.
Nowacki
F.
Nucl. Phys. A
,
2001
, vol.
694
pg.
157

()
7
Honma
M.
Otsuka
T.
Brown
B. A.
Mizusaki
T.
Phys. Rev. C
,
2002

65, 061301(R).
8
Honma
M.
Otsuka
T.
Brown
B. A.
Mizusaki
T.
Phys. Rev. C
,
2004
, vol.
69
pg.
034335

()
9
Navrátil
P.
Vary
J. P.
Barrett
B. R.
Phys. Rev. Lett.
,
2000
, vol.
84
pg.
5728

()
10
Navrátil
P.
Vary
J. P.
Barrett
B. R.
Phys. Rev. C
,
2000
, vol.
62
pg.
054311

()
11
Quaglioni
S.
Navrátil
P.
Phys. Rev. Lett.
,
2008
, vol.
101
pg.
092501

()
12
Quaglioni
S.
Navrátil
P.
Phys. Rev. C
,
2009
, vol.
79
pg.
044606

()
13
SciDAC Review, Issue 6
2007
(pg.
42
-
51
)
14
Honma
M.
Mizusaki
T.
Otsuka
T.
Phys. Rev. Lett.
,
1996
, vol.
77
pg.
3315

()
15
Honma
M.
Mizusaki
T.
Otsuka
T.
Phys. Rev. Lett.
,
1995
, vol.
75
pg.
1284

()
16
Koonin
S. E.
Dean
D. J.
Langanke
K.
Phys. Rep.
,
1997
, vol.
278
pg.
1

()
17
Otsuka
T.
Honma
M.
Mizusaki
T.
Shimizu
N.
Utsuno
Y.
Prog. Part. Nucl. Phys.
,
2001
, vol.
47
pg.
319

()
18
Shimizu
N.
Otsuka
T.
Mizusaki
T.
Honma
M.
Phys. Rev. Lett.
,
2001
, vol.
86
pg.
1171

()
19
Schmid
K. W.
Prog. Part. Nucl. Phys.
,
2004
, vol.
52
pg.
565

()
20
Otsuka
T.
Honma
M.
Mizusaki
T.
Phys. Rev. Lett.
,
1998
, vol.
81
pg.
1588

()
21
Shimizu
N.
Utsuno
Y.
Abe
T.
Otsuka
T.
RIKEN Accel. Prog. Rep.
,
2010
, vol.
43
pg.
46

22
Shimizu
N.
Utsuno
Y.
Mizusaki
T.
Otsuka
T.
Abe
T.
Honma
M.
Phys. Rev. C
,
2010

82, 061305(R).
23
Otsuka
T.
Fujimoto
R.
Utsuno
Y.
Brown
B. A.
Honma
M.
Mizusaki
T.
Phys. Rev. Lett.
,
2001
, vol.
87
pg.
082502

()
24
Otsuka
T.
Suzuki
T.
Fujimoto
R.
Grawe
H.
Akaishi
Y.
Phys. Rev. Lett.
,
2005
, vol.
95
pg.
232502

()
25
Otsuka
T.
Suzuki
T.
Honma
M.
Utsuno
Y.
Tsunoda
N.
Tsukiyama
K.
Hjorth-Jensen
M.
Phys. Rev. Lett.
,
2010
, vol.
104
pg.
012501

()
26
Otsuka
T.
Suzuki
T.
Holt
J. D.
Schwenk
A.
Akaishi
Y.
Phys. Rev. Lett.
,
2010
, vol.
105
pg.
032501

()
27
Brown
B. A.
Prog. Part. Nucl. Phys.
,
2001
, vol.
47
pg.
517

()
28
Shimizu
N.
Utsuno
Y.
Mizusaki
T.
Otsuka
T.
Abe
T.
Honma
M.
AIP Conf, Proc.
,
2011
, vol.
1355
pg.
138

29
Shirokov
A. M.
Vary
J. P.
Mazur
A. I.
Weber
T. A.
Phys. Lett. B
,
2007
, vol.
644
pg.
33

()
30
Shirokov
A. M.
Vary
J. P.
Mazur
A. I.
Zaytsev
S. A.
Weber
T. A.
Phys. Lett. B
,
2005
, vol.
621
pg.
96

()
31
Shirokov
A. M.
Vary
J. P.
Mazur
A. I.
Zaytsev
S. A.
Weber
T. A.
Subroutines to generate this interaction in the relative-center-of-mass HO basis are available at nuclear.physics.iastate.edu
2005
32
Gloeckner
D. H.
Lawson
R. D.
Phys. Lett. B
,
1974
, vol.
53
pg.
313

()
33
Honma
M.
Brown
B. A.
Mizusaki
T.
Otsuka
T.
Nucl. Phys. A
,
2002

704, 134c.
34
Honma
M.
Otsuka
T.
Brown
B. A.
Mizusaki
T.
Phys. Rev. C
,
2002

65, 061301(R).
35
Schmid
K. W.
Prog. Part. Nucl. Phys.
,
2004
, vol.
52
pg.
565

()
36
Puddu
G.
Eur. Phys. J. A
,
2007
, vol.
34
pg.
413

()
37
Puddu
G.
J. Phys. G
,
2006
, vol.
32
pg.
321

()
38
Numerical Recipes in Fortran 77, the Art of Scientific Computing
,
1992

(Cambridge University Press, Cambridge, UK), 2nd ed
39
Egido
J. L.
Lessing
J.
Martin
V.
Robledo
L. M.
Nucl. Phys. A
1995
, vol.
594
pg.
70

40
Horoi
M.
Volya
A.
Zelevinsky
V.
Phys, Rev. Lett.
,
1999
, vol.
82
pg.
2064

()
41
Horoi
M.
Brown
B. A.
Zelevinsky
V.
Phys. Rev. C
,
2003
, vol.
67
pg.
034303

()
42
Mizusaki
T.
M.
Phys, note = Rev. C
,
2002
, vol.
65
pg.
064319

()
43
Mizusaki
T.
M.
Phys, Rev. C
,
2003
, vol.
67
pg.
041301

()
44
Papenbrock
T.
Dean
D. J.
Phys, Rev. C
,
2003
, vol.
67
pg.
051303(R)

()
45
Shen
J. J.
Zhao
Y. M.
Arima
A.
Yoshinaga
N.
Phys, Rev. C
,
2011
, vol.
83
pg.
044322

()
46
M.
Kashima
T.
J, Phys. Soc. Jpn.
,
2000
, vol.
69
pg.
2723

()
47
Richter
W. A.
van der Merwe
M. G.
Julies
R. E.
Brown
B. A.
Nucl, Phys. A
,
1991
, vol.
523
pg.
325

()
48
Shimizu
N.
Utsuno
Y.
Mizusaki
T.
Honma
M.
Tsunoda
Y.
Otsuka
T.
Phys. Rev. C
,
2012
, vol.
85
pg.
054301

()
49
Mizusaki
T.
Otsuka
T.
Utsuno
Y.
Honma
M.
Sebe
T.
Phys, Rev. C
,
1999
, vol.
59
pg.
R1846

()
50
Honma
M.
Otsuka
T.
Mizusaki
T.
Hjorth-Jensen
M.
Phys, Rev. C
,
2009
, vol.
80
pg.
064323

()
51
Utsuno
Y.
Shimizu
N.
Otsuka
T.
Abe
T.
arXiv:1202.2957 [nucl-th]
54
Abe
T.
Maris
P.
Otsuka
T.
Shimizu
N.
Utsuno
Y.
Vary
J. P.

arXiv:1204.1755 [nucl-th]
55
Epelbaum
E.
Glöckle
W.
Meissner
U.-G.
Nucl, Phys. A
,
1998
, vol.
637
pg.
107

()
56
Epelbaum
E.
Glöckle
W.
Meissner
U.-G.
Nucl, Phys. A
,
2000
, vol.
671
pg.
295

()
57
Entem
D. R.
Machleidt
R.
Phys. Rev. C
,
2003
, vol.
68
pg.
041001(R)

()
58
Wiringa
R. B.
Stoks
V. G. J.
Schiavilla
R.
Phys. Rev. C
,
1995
, vol.
51
pg.
38

()
59
Pieper
S. C.
Pandharipande
V. R.
Wiringa
R. B.
Carlson
J.
Phys. Rev. C
,
2001
, vol.
64
pg.
014001

()
60
Pieper
S. C.
AIP Conf. Proc.
,
2008
, vol.
1011
pg.
143

61
Pieper
S. C.
Wiringa
R. B.
Carlson
J.
Phys. Rev. C
,
2004
, vol.
70
pg.
054325

()
62
Nollett
K. M.
Pieper
S. C.
Wiringa
R. B.
Carlson
J.
Hale
G. M.
Phys. Rev. Lett.
,
2007
, vol.
99
pg.
022502

()
63
Pieper
S. C.
proceedings of the International School of Physics “Enrico Fermi”
,
2008

Course CLXIX, in eds. A. Covello, F. Iachello, R. A. Ricci (Società Italiana di Fisica, Bologna), p. 111
64
Pieper
S. C.
Nuovo Cimento
2008
, vol.
31
pg.
709

and references therein
65
Hagen
G.
Papenbrock
T.
Hjorth-Jensen
M.
Phys. Rev. Lett.
,
2010
, vol.
104
pg.
182501

and references therein. ()
66
Okubo
S.
Phys. Theor. Phys.
,
1954
, vol.
12
pg.
603

()
67
Lee
S. Y.
Suzuki
K.
Phys. Lett. B
,
1980
, vol.
91
pg.
173

()
68
Suzuki
K.
Lee
S. Y.
Prog. Theor. Phys.
,
1980
, vol.
64
pg.
2091

()
69
Suzuki
K.
Okamoto
R.
Phys. Theor. Phys.
,
1983
, vol.
70
pg.
439

()
70
Bogner
S. K.
Furnstahl
R. J.
Schwenk
A.
Prog. Part. Nucl. Phys.
,
2010
, vol.
65
pg.
94

()
71
Maris
P.
Vary
J. P.
Navratil
P.
Ormand
W. E.
Nam
H.
Dean
D. J.
Phys. Rev. Lett.
,
2011
, vol.
106
pg.
202502

()
72
Roth
R.
Phys. Rev. C
,
2009
, vol.
79
pg.
064324

()
73
Roth
R.
Binder
S.
Vobig
K.
Calci
A.
Langhammer
J.
Navratil
P.

arXiv:1112.0287 [nucl-th]
74
Dytrych
T.
Sviratcheva
K. D.
Bahri
C.
Draayer
J. P.
Vary
J. P.
Phys. Rev. Lett.
,
2007
, vol.
98
pg.
162503

()
75
Dytrych
T.
Sviratcheva
K. D.
Bahri
C.
Draayer
J. P.
Vary
J. P.
J. Phys. G
,
2008
, vol.
35
pg.
095101

()
76
Dytrych
T.
Sviratcheva
K. D.
Draayer
J. P.
Bahri
C.
Vary
J. P.
J. Phys. G
,
2008
, vol.
35
pg.
123101

()
77
Abe
T.
Maris
P.
Otsuka
T.
Shimizu
N.
Utsuno
Y.
Vary
J. P.
Conf. Proc.
,
2011
, vol.
1355
pg.
173

78
Puddu
G.

arXiv:1201.0600 [nucl-th]
79
Liu
L.
Otsuka
T.
Shimizu
N.
Utsuno
Y.
Roth
R.
Phys. Rev. C.
,
2012
, vol.
86
pg.
014302

()
80
Utsuno
Y.
Shimizu
N.
Otsuka
T.
Abe
T.

arXiv:1202.2957 [nucl-th]
81
Maris
P.
Vary
J. P.
Shirokov
A. M.
Phys. Rev. C
,
2009
, vol.
79
pg.
014308

()
82
Maris
P.
Shirokov
A. M.
Vary
J. P.
Phys. Rev. C
,
2010
, vol.
81
pg.
021301(R)

()
83
Cockrell
C.
Vary
J. P.
Maris
P.

arXiv:1201.0724
84
Wiringa
R. B.
Pieper
S. C.
Carlson
J.
Pandharipande
V. R.
Phys. Rev. C.
,
2000
, vol.
62
pg.
014001

()
85
Itagaki
N.
Masui
H.
Ito
M.
Aoyama
S.
Phys. Rev. C
,
2005
, vol.
71
pg.
064307

()
86
Cockrell
C.
Vary
J. P.
Maris
P.

arXiv:1201.0724 [nucl-th]
87
Lenzi
S. M.
Nowacki
F.
Poves
A.
Sieja
K.
Phys. Rev. C
,
2010
, vol.
82
pg.
054301

()
88
Honma
M.
Otsuka
T.
Brown
B. A.
Mizusaki
T.
Eur. Phys. J. A
,
2005
, vol.
25
pg.
499

()
89
Hjorth-Jensen
M.
Kuo
T. T. S.
Osnes
E.
Phys. Rev.
,
1995
, vol.
261
pg.
125

90
Hjorth-Jensen
M.
Kuo
T. T. S.
Osnes
E.
National Nuclear Data Center, information extracted from the NuDat 2 database
1995

91
Pritychenko
B.
Choquette
J.
Horoi
M.
Karamy
B.
Singh
B.

arXiv:1102.3365v2
92
Porquet
M.-G.
Sorlin
O.
Phys. Rev. C.
,
2012
, vol.
85
pg.
014307

()
93
Utsuno
Y.
Otsuka
T.
Mizusaki
T.
Honma
M.
Phys. Rev. C.
,
1999
, vol.
60
pg.
054315

()
95
Ring
P.
Schuck
P.
The Nuclear Many-Body Problem
,
1980

(Springer, New York) ()
96
Puddu
G.

arXiv:1201.0600 [nucl-th]
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.