Constructing approximate shell-model wavefunctions by eigenvector continuation

Shell-model calculations play a key role in elucidating various properties of nuclei. In general, those studies require a huge number of calculations to be repeated for parameter calibration and quantifying uncertainties. To reduce the computational burden, we propose a new workflow of shell-model calculations using a method called eigenvector continuation (EC). It enables us to efficiently approximate the eigenpairs under a given Hamiltonian by previously sampled eigenvectors. We demonstrate the validity of EC as an emulator of the valence shell-model, including first application of EC to electromagnetic transition matrix elements. Furthermore, we propose a new usage of EC: preprocessing, in which we start the Lanczos iterations from the approximate eigenvectors, and demonstrate that this can accelerate subsequent research cycles. With the aid of the EC, the eigenvectors obtained during the parameter optimization are not necessarily to be discarded, even if their eigenvalues are far from the experimental data. Those eigenvectors can become accumulated knowledge.

Shell-model calculations play a key role in elucidating various properties of nuclei. In general, those studies require a huge number of calculations to be repeated for parameter calibration and quantifying uncertainties. To reduce the computational burden, we propose a new workflow of shell-model calculations using a method called eigenvector continuation (EC). It enables us to efficiently approximate the eigenpairs under a given Hamiltonian by previously sampled eigenvectors. We demonstrate the validity of EC as an emulator of the valence shell-model, including first application of EC to electromagnetic transition matrix elements. Furthermore, we propose a new usage of EC: preprocessing, in which we start the Lanczos iterations from the approximate eigenvectors, and demonstrate that this can accelerate subsequent research cycles. With the aid of the EC, the eigenvectors obtained during the parameter optimization are not necessarily to be discarded, even if their eigenvalues are far from the experimental data. Those eigenvectors can become accumulated knowledge.

Introduction
The configuration interaction method plays a key role to study various properties of manyfermion systems. Nuclei are no exception; the nuclear shell model is one of the most powerful approaches to a microscopic understanding of static properties of nuclei [1][2][3]. Its application covers a wide range of the nuclear chart, and it also provides a bridge between ab initio and phenomenological studies with the development in many-body methods to derive the shell-model Hamiltonians: Over the past decades, we have witnessed much progress in this direction with a nuclear force within chiral effective field theory (chiral EFT) [4,5]. A representative and up-to-date example of the many-body method is the valence-space in-medium similarity renormalization group (VS-IMSRG) [6][7][8][9][10], which has enlarged the reach of shellmodel studies with non-empirical interactions to nearly 700 isotopes [11]. Simultaneously, there have been much progress by other approaches such as coupled-cluster effective interaction [12,13], no-core shell-model and Okubo-Lee-Suzuki transformation [14,15], shell-model coupled cluster [16,17], and so on.
Then, the original eigenpairs can be approximated as log TBTD( f i; abcd; J ab J cd ; λ) A(cd; J cd M cd ) ≡ (−1) Jcd+Mcd A( j c j d ; J cd − M cd ), TBTD ≡ 2J ab + 1 2J i + 1 TBTD( f i; abcd; J ab J ab ; 0), H i, j = ψ( c i )|H( c )|ψ( c j ) , Then, the original eigenpairs can be approximated as log OBTD k ≡ 2 j k + 1 2J i + 1 OBTD(ii; j k j k ; 0) = ψ J TBTD( f i; abcd; J ab J cd ; λ) Then, the original eigenpairs can be approximated as log TBTD( f i; abcd; J ab J cd ; λ) Then, the original eigenpairs can be approximated as log Then, the original eigenpairs can be approximated as log

Introduction
The typical parametrization of the shell-model Hamiltonian is the following: Then, the original eigenpairs can be approximated as log  In practice, however, it is still necessary to repeat a huge number of shell-model calculations to optimize the effective interactions, quantify uncertainties, extract physical intuitions by comparison with experimental observations, test the validity of many-body methods such as VS-IMSRG, and so on. One of the aims of this work is to propose a new workflow of shellmodel calculations, i.e., how such a research procedure using shell-model calculations can be accelerated. In Fig. 1, we show a schematic picture of the workflow. As will be detailed in the followings, the key ingredient is the eigenvector continuation (EC).
While much efforts have been spent on calibrating input parameters for nuclear models (e.g., low-energy constants in chiral EFT and phenomenological shell-model interactions), one still needs expensive numerical samplings in a high-dimensional parameter space to find a reasonable range of those parameters and to quantify the associated uncertainties in the parameters and the target observables for a deeper understanding of properties of nuclei.
Recently, the importance of the uncertainty quantification (UQ) has been widely pointed out in various contexts, such as parameter calibration of the chiral EFT potentials [18][19][20] and nuclear observables [21][22][23][24][25]. This also applies to the shell model, but the UQ for the valence shell-model studies are limited up to the sd shell [26,27]. This is due to the rapid growth in the size of matrices to be diagonalized as the number of valance nucleons or the size of model space (called valence space too) increases. To reduce the computational cost for UQ, an approximation method using the principal component analysis was proposed in Ref. [27]. However, one still needs additional efforts and/or efficient methods for parameter calibration and UQ to enlarge the scope to a heavier or neutron/proton-rich region.
The paper is organized as follows. In Section 2, we explain some basics of shell-model calculations and the eigenvector continuation. We show a couple of applications of EC as an emulator and as a preprocessing in Section 3, and conclude and give an outlook in Section 4. Technical details on formulations, codes, and results are given in Supplementary Material [28] for reproducibility.
The typical parametrization of the shell-model Hamiltonian is the following: where c † a and c c denote a creation and annihilation operator on the single-particle state a and c, respectively. One can naturally extend the above to include the three-body (and higher many-body) term, but we will not go into the details of three-body forces in this work.
The eigenvectors, called wave functions too, are expressed as a superposition of many-body configurations represented by Slater determinants. We employ the so-called M -scheme basis to express the wave function. This means that the single particle states, a, b, c, d in Eq. (2), are specified by its harmonic oscillator quanta, n, l, j (and the z component of angular momentum m z and of isospin t z , if needed). Under this, the two-body part of shell-model Hamiltonian, H (2) , can be rewritten in a more explicit form: where the label of single particle states, a, b, c, d, denote the quanta {n, l, j, t z }, (j a m a j b m b |JM ), (j c m c j d m d |JM ) are Clebsch-Gordan coefficients, and J and M are coupled angular momentum and its z component. Due to the symmetries of the nuclear force, the number of non-zero V J (abcd) is typically on the order of tens or thousands, and can be further reduced by a factor of about two when assuming the isospin symmetry of the effective interaction. Usually, the diagonal part of h (1) (c = a) is called single-particle energies (SPEs), and V J (abcd) in Eq. (3) is called two-body matrix elements (TBMEs). These SPEs and TBMEs are given in the interaction file for shell-model codes. Note that, in some codes, the so-called isospin formalism V JT (abcd) is used for interactions with the isospin symmetry.
A typical recipe to determine these SPEs and TBMEs is the following: (i) making an initial guess in some way, (G-matrix, VS-IMSRG, etc.), and (ii) modifications of the parameters so as to minimize the chi-square deviation from the selected experimental data. There have been many previous works adopting the above recipe [29][30][31][32][33][34][35]. One of the most successful and well-known examples is the USDB interaction for the sd-shell nuclei [29]. The USDB gives the root-mean-square deviation 126 keV for 380 energy data in the mass region A = 16 − 40.
Once the interaction is given and the model space is specified, the main task of shell-model codes is diagonalizing the shell-model Hamiltonian under the basis states, which are now in 3/20 the M -scheme. The matrices are in general very sparse, but the total number of non-zero matrix elements is too large to be stored on memory. For this reason, the diagonalization is typically done by means of the Lanczos method [36].

Implementation: ShellModel.jl
Although many CI codes for many-nucleon systems written in Fortran are already available (e.g., NuShellX [37], BIGSTICK [38,39], ANTOINE [40], MFDn [41][42][43], MSHELL64 [44] and KSHELL [45,46]), we developed a new code ShellModel.jl [47] written in the Julia language [48,49] for efficient samplings of shell-model results and demonstrating the benefit of the methodology in the present paper. The Julia language became highly popular in the scientific computing community because it provides interactive and dynamic behavior and high-productivity like Python, high-performance, and flexibility to combine other modules (e.g., visualization, neural networks, etc.) with its own package management system.
Most of the functions in ShellModel.jl is originally implemented combining other Julia packages, but we referred to MSHELL64 and KSHELL for implementation of some methods. More precisely, we implemented the followings: the Thick-Restart Lanczos method, block family of the Lanczos methods, and the double Lanczos method [44,50] for the projection to a specified total angular momentum J. While the KSHELL code adopts the so-called onthe-fly generation of the matrix elements for the matrix-vector product, ShellModel.jl stores matrix elements rather explicitly on memory to balance the speed and moderate memory usage. The proton-neutron interaction is stored in the form of "one-body jumps" [38,39]. The current version of ShellModel.jl is optimized to run on a relatively smaller environment such as a laptop. The main codes, inputs, and sample scripts used in this study are available on the GitHub repository [47] for reproducibility. We hope that it will facilitate other future studies.

Eigenvector continuation: an efficient emulator and a preprocessing for the nuclear shell model
The eigenvector continuation (EC), which was proposed in Ref. [51], has been widely used as an efficient emulator of nuclear many-body methods [52][53][54][55][56][57][58] and as a resummation method [59,60]. In the following, we give a brief overview of the eigenvector continuation, and then explain how to apply it to shell-model calculations.
Suppose that the eigenvectors for Eq. (1) have been already obtained at N s different parameters |ψ( c 1 ) , |ψ( c 2 ) , . . . , |ψ( c Ns ) . These eigenvectors will be hereafter referred to as samples. In such a case, the eigenpairs under a given Hamiltonian H( c ) are approximated by solving the following generalized eigenvalue problem in the subspace spanned by the samples: 4/20 Then, the original eigenpairs can be approximated as The problem is reduced from the diagonalization of a sparse Hamiltonian matrix H with size of M -scheme dimension to theH of a dense matrix of size N s , which is typically on the order of tens or hundreds. This might significantly reduce the computational cost compared to the original problem. Note that it is straightforward to extend Eqs. (7)-(11) to include excited states. With the approximated wave functions, one can approximate the expectation values for a target observable: where the operatorÔ can be e.g., electromagnetic transition operators.
The eigenvector continuation has been well known as the Rayleigh-Ritz method in (applied) mathematics, and introduced to or re-evaluated in the nuclear physics community recently. However, its properties are still under consideration [61], and the previous works in nuclear physics have proved for the first time its efficiency for large-scale many-body problems and enlightened the possibility to enlarge the scope of microscopic studies to a heavier and/or more exotic region of the nuclear landscape.
For solving the EC problems, the most time-consuming part is, in general, to evaluate the expectation value of Hamiltonians, i.e.,H in Eq. (8). However, the computational cost for evaluatingH can be somewhat alleviated in shell-model calculations. This is because of the fact that the expectation values of shell-model Hamiltonians can be factorized out by the SPEs and TBMEs:H where k denotes the labels of one-and two-body interactions in a given shell-model Hamiltonian H( c ), and the OBTD and TBTD are the abbreviations of one-and two-body transition densities, respectively. The concrete expressions for the OBTD k and TBTD k are summarized in Sec. A of Supplementary Material [28]. Once we evaluate the transition densities in Eq. (13) for arbitrary two sample eigenvectors, one can significantly reduce the computational cost to evaluateH for different parameters, because these can be evaluated simply by the dot product of the parameters and the transition densities, which are independent on the parameter values. Although the number of required samples depends on the problem and desired accuracy, the cost for solving generalized eigenvalue problem in Eq. (7) is typically negligible compared to that of the original eigenvalue problem. Hence, the eigenvector continuation provides an efficient emulator of the shell-model calculation as well as other nuclear models discussed in the previous works [52,53,[53][54][55][56][57][58] In addition to its role as an emulator, EC can also be used as a preprocessing method. In the Lanczos method, one usually starts with the Lanczos iterations from a random vector, when one has no prior information about the starting vector. The convergence of the Lanczos iterations should be accelerated by starting from a good initial guess of the eigenvectors. 5/20 The approximate eigenvectors constructed by EC would make it, and it will be also helpful to increase the sample size for further improvement of the accuracy the EC emulator.

Results
In this section, we present the main results of this study. For demonstration purposes, we restrict ourselves to mainly consider the sd shell on top of 16 O core, i.e., the model space consists of 1s 1/2 , 0d 3/2 , 0d 5/2 orbits (12 single-particle states in total) for both protons and neutrons. Since we consider only the positive parity states, the plus sign on the total angular momentum J, e.g., 0 + , 2 + , will be omitted where appropriate.
In Sec. 3.1, we briefly describe our procedures to prepare the sample eigenvectors. Next, the efficiency of the EC emulator is demonstrated in Sec. 3.2 and Sec. 3.3. We also make some remarks on the accuracy of approximated eigenvectors in terms of magnetic dipole moments and electric quadrupole moments. In Sec. 3.5, we explain the preprocessing for the Lanczos method and discuss its possible accelerating effect on the convergence. In Sec. 3.6, we reexamine the accuracy of EC approximate wavefunctions utilizing the preprocessing. Finally, we will mention the extensibility of the proposed method for larger systems in Sec. 3.7.

Calculations of sd-shell nuclei with ShellModel.jl
Our new shell-model code, ShellModel.jl, is designed to make sample eigenvectors efficiently and to demonstrate the efficiency of the emulator and the preprocessing with EC. With ShellModel.jl, the execution time to evaluate 10 lowest states of 28 Si, which has the largest M -scheme dimension in the full sd-shell space (93,710 for M = 0), is about 3 seconds on a MacBook Air (2020, Apple M1). We provide the sample codes in the GitHub repository [47], and further instructions for the sample codes are given in Sec. B of Supplementary Material [28].
To make the sample eigenvectors, we generated 50 random sd-shell interactions by adding gaussian random values with the standard deviation σ int. = 1.0 MeV to the USDB interaction [29]. Hereafter, the terms "random interaction" and "sample" denote those by σ int. = 1 unless otherwise mentioned. In the next section, we also discuss the results with σ int. = 3 and samples utilizing Latin Hypercube Sampling [62]. The USDB is a phenomenological sd-shell interaction constructed by G-matrix [63] and the chi-square fit using 608 data in 77 nuclei. In accordance with the USDB interaction having the isospin symmetry, the Hamiltonians are defined in the 66-dimensional parameter space (3 for SPEs and 63 for TBMEs). Next, we calculated the five lowest eigenvectors with J = 0, 1, 2, 3, 4 (for even nuclei) and J = 1/2, 3/2, 5/2, 7/2, 9/2 (for odd) under the 50 random interactions, i.e. sampling 250 eigenvectors in total for each J.

Approximate eigenpairs by the eigenvector continuation
In what follows, we show the results of approximate eigenpairs by EC. Using the samples discussed in the previous subsection, we solved Eq. (7) to estimate the approximate eigenvalues for other 100 random interactions, which were made in the same way as above, and compare them to the exact values. In the followings, the 100 random interactions are referred to as the validation set.
In Fig. 2 In Table 1, we summarize the typical size of two types of errors by EC estimates compared to the exact results. One is the relative error of absolute energy values, and the other one is error in terms of the excitation energies: ex. error (MeV) ≡ |E ex. exact − E ex. EC |.
Since excitation energies themselves are dependent on the level ordering, we restrict ourselves to measure excitation energies from the 0 + (1/2 + ) state for even (odd) nuclei. As shown in Table 1, the typical size of the relative errors is less than 1%, when all the N s = 250 samples for each J is used. Here we chose the four nuclei in the middle of sd shell ( 28 Si, 26 Al, 25 Mg, and 24 Mg) as examples. We found that the relative errors for the odd-odd nucleus, 26 Al, are worse than those for the others.
In the context of nuclear structure, the excitation energy may be of interest rather than the absolute value of the energy. As shown in Table 1, the typical size of the error lies on the order of a few hundred keV or 1 MeV. The errors for excitation energies show relatively slower convergence to the exact ones compared to their absolute energies. This can be general tendency because the excitation energies are sensitive to the relative convergence speed of the EC wavefunctions. The origin of this tendency will be examined in the Sec. 3.7.
To show the dependence on the sampling procedure, we summarized the results for five different N s in Table 1. The sample size N s is the product of the number of random interactions and the number of excited states used as the sample eigenvectors for EC, and some cases are marked with an asterisk to indicate that the detailed conditions differ from the others. For example, the N s = 50 * case means that only the first and second lowest states under the first half of the 50 random interactions are used, while the N s = 250 * (σ int. = 3) case corresponds to the result using the samples, where the size of Gaussian noise σ int. increased from 1 to 3. The eigenvalues of the five yrast states of 28 Si under the N s = 250 * (σ int. = 3) spread over a wide range of −350 to −50 MeV, but the relative errors are still ∼ 3%. This indicates that Table 1 Average size of two errors by EC estimates for the five yrast states of the four sd-shell nuclei: One is the relative error (%) of absolute energies, and the other one is error of excitation energies. The sample size N s means the product of the number of random interactions and the number of excited states used as the sample eigenvectors in Eqs. (8)(9). The N s = 250 * with σ int. = 3 means that the standard deviation to generate the random interactions is increased from the default value σ int. = 1, and the N s = 250 * (LHS, L = 2) corresponds to the result using Latin Hypercube Sampling (LHS).  one can make use of the eigenvector continuation for rough parameter optimization even if the prior knowledge of the parameter range is poor. Besides the random samples generated with Gaussian noises, we tested the dependence on the sampling procedure by utilizing Latin Hypercube Sampling (LHS) [62]. To compare with gaussian random samples, we fix the lattice size L for LHS to be 2, i.e. to cover 2σ (95%) domain with σ int. = 1 case. As shown in Table 1, the result of N s = 250 * (LHS) samples show almost the same accuracy with the N s = 250 case. One possible advantage of using LHS for EC is that it provides numerical stability for generalized eigenvalue problems, since the distance between each LHS sample tends to be larger than that for the Gaussian case. In this study, the generalized eigenvalue problem, Eq. (7) is not ill-conditioned with the help of relatively high dimension, D = 66. Hence, we mainly focus on the results with random samples generated with gaussian noise for brevity of the following discussions.
From the results above, one can see that there are two strategies to improve the accuracy of the EC estimates. One is to sample more excited states; the N s = 250 case using the five lowest states as the samples improves the accuracy by a factor of two in terms of the relative errors from the N s = 50 case. The other one is to increase the sample points in the parameter space. Looking at the results of the N s = 50 and N s = 50 * cases, we expect that the latter one is more beneficial unless we change the number of states of interest. On the other hand, the former one is relatively easy in terms of computational costs. Since the efficiency may depend on the sampling procedure and the distribution of eigenvalues, how to select the strategy and how many states to be sampled will be determined by the trade-off between the costs and desired accuracy.
Note that one can receive benefit from the factorization of transition densities in Eq. (13) during the sampling procedure; to increase the samples, all one needs is evaluating transition densities among already sampled eigenvectors and the new sample. 8/20 The execution time to calculate approximate eigenvalues for a random sd-shell interaction is about a few tens of milliseconds in our settings. It significantly facilitates the sampling procedures for the parameter optimization or uncertainty quantification.

Uncertainty quantification
In this subsection, we consider the Monte Carlo sampling over the parameter space utilizing the EC emulator. The following analysis is by no means an exhaustive study of quantifying uncertainties in the sd-shell effective interactions, but it is meant to illustrate the usefulness of the emulator.
We simplify the problem by picking up twelve energy data of 28 Si from those used to construct the USDB [29]: 0 + 1−2 , 1 + 1−2 , 2 + 1 , 3 + 1−4 , 4 + 1−3 . This is apparently too small data set to constrain the 66 parameters in the sd-shell interaction (in the isospin formalism), but enough for the demonstration purpose.
For this small data set, we consider to get samples obeying a certain posterior distribution. Let the log-likelihood be where E EC,i ( c), E Exp.,i , and σ err,i are the EC estimate of energy, its corresponding experimental value, and the expected typical error of the calculation. Since the errors stemming from EC are dominant, we fixed the σ err,i as follows where the first term is fixed as 1.0 MeV, which is typical error in order of 1% of the energy eigenvalues, and the second term is small correction with i dependence, which is taken from the difference of E EC,i with N s = 250 under the USDB interaction from the corresponding exact value. Since the number of data is less than the number of parameters, we introduced the independent Gaussian prior, which is equivalent to the L2 regularization, to avoid highly multi-modal structure of the log-likelihood in the parameter space. The log-prior is defined as Here || · || 2 is the L2 norm, i.e., the sum of squared deviation between a parameter c and the reference value c ref. , which is now the USDB. The choice of variance of the Gaussian prior Λ is non-trivial, but we fixed Λ = 10 to simplify the discussions. To sample from the posterior ∝ L( c) × P r( c), we adopted a Markov chain Monte Carlo (MCMC) method with the Metropolis-Hastings algorithm. We set a random interaction around the USDB parameters as an initial state of the MCMC, and generated 100,000 MCMC samples (after 10,000 burn-in) with the Metropolis-Hastings updates. The elapsed time to achieve it with the EC emulator on a laptop was about 30 minutes. In Fig. 3, we show one realization of parameter-distributions for some selected channels. Supplementary Material [28], and one can see that the independent runs give the converged energy distributions.
The posterior is similar to the prior, but some are slightly modified through the loglikelihood, Eq. (16). This indicates that the USDB values for these channels were affected and thereby constrained by a variety of data, not just the 28 Si data with lower J. In this way, the efficient emulator allows a detailed sensitivity analysis on the parameter and adopted data.
There still remain a lot of investigations undone such as dependence on the choices of the loss function (log-likelihood), experimental data, and σ err,i , and the more efficient sampling schemes for MCMC. In addition to these, it is also desired to develop a method to quantify the uncertainties taking into account of the variational property of the EC estimates for the energy eigenvalues. However, the above analysis manifests that one can achieve the samplings, which were computationally too heavy to carry out on a laptop, in tens of minutes with the help of the EC emulator. It is encouraging to proceed towards full evaluation of uncertainties for deeper understanding of the relation among nuclear structure, effective interactions, the choice of the experimental data, and the underlying nuclear force.

Accuracy of the approximate wave functions
Next, we test the accuracy of the approximate wave functions by evaluating other observables using Eq. (12). In Fig. 4, we show the approximated values of magnetic dipole moments µ and electric quadrupole moments Q, which are calculated for the validation set, of the yrast states of the four nuclei in the middle of sd shell. The symbols denote the yrast states with J ranging 1/2, 1, ..., 9/2 which are allowed by the selection rules, i.e., J = 0 for µ and J = 0, 1/2 for Q. The values are given in units of nuclear magneton µ N = 0.105 e fm for µ and e fm 2 for Q, and we used the free-nucleon values for the g-factors (g p = 1.0, g n = 0, g sp = 5.586, and g sn = −3.826) and effective charges (e p = 1.5, e n = 0.5), while these choices do not have much effect on the correlation plot.
One can see that the agreement between the EC estimates and the exact ones for both µ and Q become worse than that for energy eigenvalues. This is simply because the electromagnetic transitions are much sensitive to the accuracy of the wave functions. Note, however, that 10 the correspondence between the exact and approximated eigenstates are deduced simply from their energy eigenvalues. As described below, it is highly non-trivial to find one by one correspondence among them.
Most of the approximate values of the magnetic dipole moment are nearly on the diagonal line in Fig. 4 (a). The points for even nuclei with J = 1, 2, 3, 4 concentrate on µ ∼ 0.5, 1.0, 1.5, 2.0, respectively. On the other hand, results of electric quadrupole moments are rather scattered, while most of the results are still on the diagonal line. As an extreme example, let us look at the lower right-most triangle in Fig. 4 (b), which is J = 4 + 1 state of 24 Mg with Q Exact = +28.340 (e fm 2 ). Diagonalizing the corresponding random interaction, we found that the first and second eigenstates have small energy difference ∼ 400 keV and, the Q-moments with opposite sign and similar amplitude. While the approximate energy eigenvalues are close to the exact ones with ∼ 2% accuracy, the Q moment by EC has the opposite sign to the exact ones: This indicates that the approximate eigenvectors for the nearly degenerate states can be in the wrong order in terms of the energy, or contaminated by other states. In addition to the energy eigenvalues, the moments and electromagnetic transitions are much useful for sensitivity analyses on the input parameters and on adopted approximations and truncations (if there were). When discussing these observables, one needs additional 11/20 efforts to improve the accuracy of the approximate wave functions. We revisit this issue in Sec. 3.6.

Acceleration of the Lanczos iterations by the preprocessing
From the results shown in the previous subsections, it is expected that the convergence of the EC estimates to the exact values is not so fast as a function of the sample size. While the eigenvector continuation provides an efficient emulator, one still needs additional efforts to discuss the observables which demand high accuracy of the wave functions.
In general, the matrix-vector product of the shell-model Hamiltonian matrix takes up most of the calculation time, and one usually starts the Lanczos iterations from a normalized random vector. In the following, we consider the preprocessing, which is to start the Lanczos iterations from the EC approximate eigenvectors, to accelerate the shell-model calculations. For simplicity, we discuss the numbers of Hamiltonian operations in the Lanczos iterations with and without the preprocessing, although the elapsed time of the computation depends not only on the number of the operations but also on various conditions such as computing environment and the block size for a block algorithm.
As an example, we consider the four sd-shell nuclei discussed above, and fix the target shell-model Hamiltonian to the USDB interaction. In Table 2, we summarized the number of iterations, which is equivalent to the number of Hamiltonian operations, required to get converged results of the lowest states with the specified J. The number of iterations for the block Lanczos method is summarized in Table 3. When considering the excited states too, we use the block algorithm for the Lanczos method and both the block size q and the number of states of interest n are fixed as four to simplify the discussions. It means that we prepare the four lowest approximate eigenvectors from the samples and then use them as the starting block vectors with the size of q = 4.
As a general trend, we can see that the number of iterations becomes smaller as the number of samples increases. This tendency is prominent in the case of the single target state, Table 2. One exception is seen in the 26 Al in Table 3; the N s = 50 * and N s = 150 cases are worse than the N s = 50 case. This suggests that in order to speed up the convergence, the initial guesses of the nearly degenerate states must be made simultaneously with high accuracy. Indeed, the exact eigenvalues of third and fourth J = 0 states of 26 Al are close to each other with a difference of about 50 keV.
It should be noted that the preprocessing can deteriorate the convergence with some naive usage. If one calculates the ten lowest J = 0 states of 24 Mg starting from four approximate eigenvectors which are constructed from the N s = 250 (50 × 5) samples, one may need additional several iterations compared to a random starting vector.
The convergence patterns of the block Lanczos iterations are shown in Fig. 5. The four rows from the top to the bottom corresponds to the case of (q, n) = (4, 4), (10,4), (4, 10), (10, 10), respectively. Here, the q is the block size and the n is the number of states of interest. The right panels are difference in a logarithmic scale from the exact value. Note that the tolerance for the Lanczos method is set as 10 −6 , and the difference smaller than the tolerance is regarded as log 10 (E Ritz − E Exact ) = −6 in the log-plot. The solid line with circles and the dashed lines with diamonds correspond to the cases of preprocessed and random initial vectors, respectively. The symbols represent the points at which the results converged. 12/20   As already discussed above, the (q, n) = (4, 4) case exhibits the effect of the preprocessing. Looking at (q, n) = (10, 4), we can see that the number of Hamiltonian operation is much reduced. On the other hand, the (q, n) = (4, 10) case shows that insufficient preprocessing can deteriorate the convergence, and the effect of the preprocessing is limited for the (q, n) = (10, 10) case due to the ninth and tenth lowest states. This is because the N s = 250 (50 × 5) samples have less information about the higher excited states.
As shown in Fig. 5 and Supplementary Figure D1-D3 [28], the efficiency of the preprocessing depends on the problem. A rule of thumb to receive the benefit of the preprocessing is starting with a larger number of states than that of interest. It does not always reduce the computation time, because the number of manipulations in each Lanczos iteration becomes larger as the block size increases. In such a case, it would be better to start with a large initial block vector and, after several iterations, compress it to a smaller block vector by the Thick-Restart method.
One can consider another preprocessing by taking eigenstates within a truncated subspace (e.g., particle-hole truncation). In general, however, the overlap between the two eigenvectors obtained within a truncated subspace and the full model space strongly depends on the states and the choice of the truncation scheme for those target states. For example, third 0 + state of 56 Ni is known to show rather slow convergence as a function of the number of allowed 13 The four rows, respectively, correspond to (q, n) = (4, 4), (10,4), (4,10), (10,10) cases, where q is the block size and n is the number of states of interest during the Lanczos iterations. The left panels are for Ritz values, and the right panels show the difference on a logarithmic scale from the exact value with the Lanczos tolerance 10 −6 . The solid lines with filled circle and dashed lines with diamond denote the cases starting from preprocessed and random initial vectors, respectively. The symbols represent the points at which the results get converged.
particle-hole excitations in the truncation scheme than yrast states [64]. The EC emulator gives a more general way to make starting vectors for nuclei and/or states without such a special consideration. 14/20 3.6. Improvement of approximate wavefunctions by the preprocessing In Sec. 3.4, we showed that the EC emulator can fail to approximate the magnetic dipole moments and electric quadrupole moments due to the level inversion for nearly degenerate levels. Here, we reconsider this issue and its improvement utilizing the preprocessing discussed above. For this purpose, taking EC emulated wave functions as the starting Lanczos vectors, we perform only two additional Lanczos iterations. We will call it EC+2, and the reevaluation of µ and Q moments with EC+2 wavefunctions are shown in Fig. 6. It is seen that the scatter of the symbols has been considerably improved with just a few additional Lanczos iterations. For some points that are still far from the diagonal line, these are, in general, corresponding to the random interactions that induce nearly degenerate states within a few tens of keV, so that their moments are very sensitive to the number of Lanczos iterations.
This result manifests the efficacy of our workflow using EC as an emulator and a preprocessor to estimate any observables of interest.

Extensibility of EC emulator: to go beyond the sd-shell
Regarding the extensibility of the EC emulator for shell-model calculations for larger systems, we explore the behavior of the accuracy of the emulator against the model space size in this subsection. In Fig. 7, we plot the relative errors defined in Eq. (14) as a function of J-scheme dimension, where d J (J) and d M (M ) are J-and M -scheme dimensions under given total angular momentum J and its projection M , respectively. It makes dependence of the accuracy on the dimension and total J clearer. As a whole, the relative errors are increasing functions of J-scheme dimension, and the lower J states show relatively faster convergence to the exact ones. Since the smaller J generally gives smaller J-scheme dimension, the results indicate that the number of effective degrees of freedom governs the accuracy. From Fig. 7   In addition to the four sd-shell nuclei discussed above, we show the results with lower pf -shell nuclei ( 46 V, and 47,48 Ti) to see the model-space dependence. For the pf -shell model space, there are 199 parameters (4 SPEs and 195 TBMEs) in isospin formalism, which is three-times larger than the sd-shell case, and we used the GXPF1A interaction [33] as the reference to prepare the random samples. The other steps such as generating samples are the same as in the case of the sd shell. As a whole, the relative errors become worse than the sd-shell case, when we use σ int. = 1 for 199 parameters. For the target nuclei ( 46 V, and 47,48 Ti), it seems that the SPEs and TBMEs involved in f 7/2 and p3/2 mainly affect the wavefunctions, so random samples made in such a way that all the parameters are independent can be less informative to span the exact eigenvectors, and thereby gives larger errors than the sd-shell case. Indeed, if we restrict the sampling space to the parameters related only to 0f 7/2 and 1p3/2 (i.e. the other parameters are fixed as the reference values), the relative errors against the validation set are improved especially for higher J states. The corresponding results are shown by the open symbols in Fig. 7.
For larger model spaces, it is expected that a strategy such as principal component analysis to sample mainly in low-dimensional subspaces that are sensitive to wavefunctions will be necessary to achieve the same level of accuracy as in the sd-shell case. Additionally, if one goes beyond the sd shell and lower pf shell, we need code development with massive parallelization and more efficient designs of the sampling scheme. Even in such cases, the 16/20 workflow re-using already sampled wave functions gives acceleration cycles for the sampling procedure, as demonstrated in the previous subsections.

Summary and Outlook
We demonstrated that the eigenvector continuation (EC) can be used as an efficient emulator for shell-model calculations, as in previous works for other models [52][53][54][55][56][57][58]. As an example, we considered sd-shell nuclei and the randomly sampled effective interactions around the USDB interaction.
With the sample eigenvectors under the 50 random interactions, we demonstrated that the exact eigenenergies under other 100 random interactions in a similar range can be estimated with an accuracy of a few percent. The EC emulator is found to be helpful for rough parameter optimization and quantifying uncertainty. The latter was partially demonstrated by Markov chain Monte Carlo sampling on a simplified problem.
It should be noted that one should pay attention to the accuracy of the approximate wave functions by EC when discussing excitation energies and observables such as electromagnetic transitions, which are sensitive to the relative convergence speed among different states. In such cases, one needs additional efforts to improve the accuracy of the approximate wave functions. As one solution to make it feasible, we proposed a new usage of the eigenvector continuation as the preprocessing for the shell-model calculations. The exact diagonalization and subsequent manipulations can be accelerated by starting the Lanczos iterations from the EC approximate eigenvectors.
As is schematically shown in Fig. 1, the eigenvector continuation can be used (a) as an efficient emulator, (b) as a preprocessing to accelerate the exact calculations, and (c) to generate another sample eigenvector efficiently with the help of the preprocessing. The eigenvectors in the optimization process of effective interactions are no longer necessarily to be discarded, even if their eigenvalues are far from the corresponding experimental data. These eigenvectors are not wasted, but can become accumulated knowledge and accelerate the workflow of researchers.
In addition to the methodology, we developed the new open-source shell-model code, ShellModel.jl, written in the Julia language. The code enables us to sample shell-model results efficiently, and is highly flexible to add any extensions to suit the users' purposes such as MCMC samplings. The main codes and sample scripts are provided in the GitHub repository [47].
There are many possible future works along the line presented in this work. While the microscopically derived effective interactions give reasonable results, it is still indispensable to construct better phenomenological interactions for more quantitative discussions in comparison with experimental data, designing new experiments, and testing the validity of many-body methods (and the approximations in them) to derive effective interactions. Of particular importance for this purpose is the extension of the proposed method to larger model spaces. This requires a more efficient sampling strategy that can work even when the computational cost for a single sample is high. Such a development is expected to facilitate the construction of phenomenological interactions and the validation of microscopically derived effective interactions, even in mass regions where gold-standard phenomenological interactions like USDB do not exist. 17/20 For parameter calibration using the EC emulators, we still have interesting topics to be further studied. Unlike phenomenological interactions used in this study, the dependence of the wavefunctions on the controlling parameters such as low-energy constants (LECs) in chiral EFT will not be factorized out as simple as Eq. (13) when a shell model effective interaction derived from a chiral potential is considered. In such a case, sensitivity analyses to the parameters will become more beneficial, but more non-trivial. If one could extend the EC emulator to treat the LECs as the controlling parameters for shell-model effective interactions, uncertainties in the EC emulator would be reduced. This is because the number of parameters is reduced from e.g., 199 (pf -shell case above) to the typical number of LECs, 30-40. Those sensitivity analyses on the LECs in medium mass regions will make the relationship between nuclear structure and nuclear force clearer.
For uncertainty quantification, efficient emulators allow to evaluate full posteriors for the parameters and/or observables. It would also be interesting to explore such a direction using more efficient MCMC methods.
It is also desirable to open a database of transition densities and wave functions so that one can easily obtain approximate eigenpairs and benefit from the preprocessing without repeating massive computations. This will further facilitate the cycle of cross-validation between theory and experiment. These would contribute to the long journey towards understanding of how nuclei emerge from the underlying nuclear force.