Efficient emulators for scattering using eigenvector continuation

Eigenvector continuation EC has been shown to accurately and efficiently reproduce ground states for targeted sets of Hamiltonian parameters. It uses as variational basis vectors the corresponding ground-state eigensolutions from selected other sets of parameters. Here we extend the EC approach to scattering using the Kohn variational principle. We first test it using a model for S-wave nucleon-nucleon scattering and then demonstrate that it also works to give accurate predictions for non-local potentials, charged-particle scattering, complex optical potentials, and higher partial waves. These proofs-of-principle validate EC as an accurate emulator for applying Bayesian inference to parameter estimation constrained by scattering observables. The efficiency of such emulators is because the accuracy is achieved with a small number of variational basis elements and the central computations are just linear algebra calculations in the space spanned by this basis.

Eigenvector continuation (EC) has been shown to accurately and efficiently reproduce ground states for targeted sets of Hamiltonian parameters. It uses as variational basis vectors the corresponding ground-state eigensolutions from selected other sets of parameters. Here we extend the EC approach to scattering using the Kohn variational principle. We first test it using a model for S-wave nucleon-nucleon scattering and then demonstrate that it also works to give accurate predictions for non-local potentials, charged-particle scattering, complex optical potentials, and higher partial waves. These proofs-of-principle validate EC as an accurate emulator for applying Bayesian inference to parameter estimation constrained by scattering observables. The efficiency of such emulators is because the accuracy is achieved with a small number of variational basis elements and the central computations are just linear algebra calculations in the space spanned by this basis.

I. OVERVIEW
Bayesian inference is increasingly favored for uncertainty quantification in nuclear physics calculations (e.g., see [1][2][3][4][5]), but the computational requirements can be substantial. In particular, Bayesian parameter estimation generally requires Monte Carlo sampling of the parameter space, with many evaluations of the likelihood with different parameters. Each evaluation may be sufficiently expensive that a full parameter estimation is infeasible. Eigenvector continuation (EC) [6,7] has already shown that it can be used as an efficient and accurate emulator [8] to ameliorate this problem. In applying an emulator, one trains computer models of the relevant calculations using a representative set of parameters and then samples for other parameters from the model instead of full calculations. Efficient and effective EC emulators for nuclear bound-state properties and transitions have been demonstrated for many-body calculations using chiral effective field theory (χEFT) Hamiltonians [8,9].
We would also like to have fast EC emulators for scattering, e.g., for treating reactions and for few-body scattering used to constrain χEFT low-energy constants [10]. The variational method for ground-state energies is well known from elementary quantum mechanics. In addition, there are variational formulations of scattering, such as those by Schwinger and Kohn (see Refs. [11][12][13][14] and references therein). The conventional applications in scattering are for two-body scattering, usually in partial waves, but the literature contains adaptations to three-body scattering, including nucleon-deuteron scattering [15,16], a process of particular interest for χEFT [10]. Here we merge EC and the Kohn variational principle and explore how well it works using a series of model calculations, starting with two-body scattering in partial waves. As * furnstahl.1@osu.edu † garcia.823@osu.edu ‡ millican.7@osu.edu § zhang.10038@osu.edu demonstrated below, a small number of variational basis based on EC can reproduce the exact calculations with great accuracy. As a result, the main computational cost is just linear algebra in this low-dimensional space.

II. FORMALISM
Consider a Hamiltonian H(θ) = T + V (θ) with adjustable parameters θ. For example, the vector θ could be the depth of a simple square well or the full set of low-energy constants for an effective field theory. EC is a variational method that employs a non-orthogonal basis composed of eigenvectors from different parameter sets {θ i } of the Hamiltonian. For calculating the ground state of H(θ), the trial wave function is where |ψ gs (θ i ) is the ground-state eigenvector of H(θ i ).
(The dependence of c i and |ψ trial on θ i is suppressed for notational convenience.) The N b θ i s are chosen either systematically or randomly to span a particular range of values, see below. The effectiveness of the EC basis can be understood by an analytic continuation analysis [6,17]. The variational principle for the ground-state energy states that the expectation value of H(θ) in the trial state, subject to the condition that |ψ trial is normalized, is stationary: The stationary solution given Eq. (1) is a generalized eigenvalue problem yielding Lagrange multiplier λ min , which is an upper bound to E gs , and the {c i } provide an approximation to |ψ gs (θ) through Eq. (1) [6][7][8][9].
For the extension of EC to scattering we use the Kohn variational principle (KVP) [13,18]. There are many variational methods for scattering, but the KVP is particularly straightforward to adapt to EC in a form similar arXiv:2007.03635v2 [nucl-th] 14 Jul 2020 to Eq. (2). Let us start with the goal of finding the phase shift δ (E) at energy E for nonrelativistc two-body scattering in an uncoupled partial-wave channel with angular momentum and short-range forces only. In coordinate space, T → −∇ 2 /2µ with = 1 and reduced mass µ, and we allow V (θ) to be local or nonlocal.
We take the trial wave function for the extended EC to be (we again suppress some θ i dependence) where |ψ E (θ i ) is the partial-wave solution for the Schrödinger equation with Hamiltonian H(θ i ) at energy E > 0, normalized such that for every i, Here p = √ 2µE, the scattering wave function is decomposed as and [13] for H(θ i ) at energy E. The KVP asserts that (also see the Supplementary Material (SM)) [13] β |ψ trial ≡ τ trial − 2µ ψ trial | H(θ) − E|ψ trial , subject to the radial part of r|ψ trial being normalized as in Eq. (4) but with K (i) (E)/p → τ trial , will be a stationary approximation to [K (E)] exact (i.e., it is accurate to second order in the difference of the exact and trial wave functions although not an upper bound in general). The normalization condition for |ψ trial is fulfilled if N b i=1 c i = 1, which can be imposed with a Lagrange multiplier λ. Substituting (3) into (6) with this constraint term and requiring the derivatives with respect to c i and λ to be zero yields a simple matrix inversion problem with solution where In obtaining Eq. (9) we have used that H(θ i ) − E |ψ E (θ i ) = 0 for every i. Finally, the stationary approximation to the exact partial-wave K matrix is Thus the approximation is given by a weighted average of the K matrices from the basis Hamiltonians with a correction term. Note that the validity of the KVP relies only on the cancellation of δτ trial with surface terms arising from the variation of ψ trial | H(θ) − E|ψ trial when varying β |ψ trial , which is satisfied by Coulomb, non-local, and complex potentials, as well as for coupled channels. (When the Coulomb potential is present, the asymptotic behavior of the scattering wave function is different from Eq. (4). For complex potentials, the ψ E (θ i )| factors in Eqs. (6) and (9) need to be applied with time reversal [19,20]. See the discussion in the SM.) It is worth pointing out that any long-range potential in H(θ) independent of θ, such as Coulomb, will cancel from ∆ U ij in Eq. (9) and one needs only to evaluate the matrix element within the range of the remaining potentials, which simplifies calculations. Also note that Eq. (9) can be evaluated in momentum space or any other convenient basis. More details on the derivation of Eqs. (7)-(10) are given in the SM.
As seen in Eqs. (7) to (9), the numerical effort is mainly composed of (a) constructing the ∆ U matrix and (b) linear algebra operations with it. The computational cost in (a) can be significantly reduced by saving the θ-independent pieces, which are also the most timeconsuming ones, instead of computing them while sampling the parameter space. For (b), the small dimension space-N b ∼ 10 in the following test examples-reduces both memory and time in the linear algebra calculations.
However, the matrix ∆ U to be inverted may be expected to be increasingly ill-conditioned as the basis size N b increases. Even for conventional applications of the KVP, there will be ill-conditioning issues for certain values of E, giving rise to so-called "Kohn anomalies" [22]. The often-recommended remedy is to use a complex formulation (involving the S matrix rather than the K matrix), which mostly avoids the problem [23][24][25]. Here we also have ill-conditioning, but at all E for sufficiently large N b . We find, however, that a simple regularization of the smallest singular values of ∆ U is sufficient to ameliorate the ill-conditioning [26,27]. This can be done by adding a small value to the diagonal of ∆ U (called a nugget in this context, but cf. Tikhonov regularization [26,27]) or by using the pseudo-inverse in Eq. (8). Because we can accurately calculate test results, we can verify the efficacy of the regularization. In the following calculations, the nugget is chosen to be between 10 −10 and 10 −8 to optimize-by hand-those EC estimations with an ill-conditioning problem. FIG. 2. Sampled points in the parameter space for the Minnesota potential in the 1 S0 channel. The best parameter set from Ref. [21] is a star, four values for the EC trial basis are circles, and the crosses are test point, which are either interpolations (blue) or extrapolations (red).

III. EC FOR A MODEL OF NN SCATTERING
We use the so-called "Minnesota potential" [21], which was developed to reproduce 1 S 0 and 3 S 1 nucleon-nucleon (NN) scattering phase shifts with a simple functional form, as a test example to explore the application of EC for scattering. The potential is a sum of local Gaussian terms, without Coulomb interaction or coupled channels. Each S-wave channel has a repulsive short-range term and an attractive term with longer range: The best values from Ref. [21] Figure 1(a) shows the scattering wave functions at E = 50 MeV (in the centerof-mass (CM) frame) from the four basis potentials (blue dot-dashed lines), the exact wave function corresponding to the best value parameters at the same energy (red dashed), and the wave function from EC based on the four-potential basis (black solid line). It is evident that with four basis elements the EC wave function agrees very well with the exact wave function. Figure 1(b) shows the corresponding phase shifts. Again, the exact result is very well reproduced by the EC prediction, even though the wave functions and phase shifts of individual basis elements are significantly different.
Next we make a more global study with the same potential. For each channel, we vary the two potential strengths by ±100 MeV about the best values, and scan the 2-dimensional parameter space by comparing the EC -phase shift with the exact phase shift. The values for the θ i = {V 0R , V 0s } parameters are randomly drawn using Latin-hypercube sampling [28], as used in EC boundstate studies [8,9]. A range of basis sizes N b have been explored. For those parameter values, we compute scattering phase shifts and wave functions by directly solving the Schrödinger equation using an R-matrix package [29], which serves as input for the subsequent EC calculations using Eqs. (7)- (9). To explore the predictive power of the EC, we randomly sampled 200 points from the twodimensional space, and for each made EC predictions as well as direct calculations using the R-matrix package, whose phase-shift calculation, as we checked, has precision-i.e., relative error-better than 10 −8 with the order of 10 2 mesh points used therein. Comparing these results indicates the accuracy of the EC emulator.
An example of the parameter sets for this comparison protocol is shown in Fig. 2, where the sampled points in the V 0R -V 0s parameter space (for the 1 S 0 channel) are shown. The trial basis points (N b = 4) are blue circles, the tested sample points are blue crosses if within the convex hull of the basis points (for these the EC calculations are considered to be interpolations) and otherwise are red crosses (these EC calculations are extrapolations), and finally the best-value point is a red star.
In Fig. 3, the mean values of the relative error (in absolute value) of the EC calculations for the interpolated sample points (left panel) and the extrapolated points (right panel) are plotted against the scattering energy E (in the CM frame) for three calculations using N b = 4, 6, and 8 basis elements. (The errors are in the value of p/K (E) = p cot δ(E). Since the relative error is tiny here, other functions of δ have almost the same relative errors.) With a basis size of 4, the EC calculation can reproduce the phase shift to better than 0.1 percent at almost all energies. The accuracy improves to be better than 10 −4 , and for most energies it reaches 10 −6 , with N b = 6. For N b = 8, the ∆ U matrix becomes ill-conditioned, but after regularizing the small singular values by adding a nugget (10 −10 ) to the diagonal of this matrix when computing the matrix inversion in Eqs. (7) and (8), the accuracy of these calculations is comparable to the N b = 6 case. We also computed the standard deviations of the absolute value of the relative errors, and found them to be similar in size to the mean values. It is interesting to note in this case that EC works equally well for interpolated and extrapolated points.
To explore a higher-dimensional parameter set, we vary both the potential strength (±100 MeV about the best values) and the two Gaussian widths κ R and κ s within a ±50% range about their best values. So now θ i = {V 0R , κ R , V 0s , κ s }. 1 For this demonstration, we uniformly sample 1000 test points within the fourdimensional parameter space. Figure 4 shows the parallel error information to Fig. 3, with a nugget of 10 −9 size. For the interpolated points, the accuracy improves from 10 −3 -10 −4 to 10 −6 or better as N b increases from 6 to 10. For N b = 14, the ill-conditioning issues require the use of a nugget but its accuracy is similar to N b = 10. The results for the extrapolated parameter sets are worse than the interpolated results for N b = 6, but become as accurate with N b = 10 and larger. Again, the standard deviations of the relative errors are comparable to their mean values. The parallel results for the 3 S 1 channel are similar for a large enough trial basis (see the SM). To explore the effectiveness of EC for non-local potentials, the inclusion of a Coulomb potential, and for higher-partial waves, we use proton-α scattering in the S 1/2 and P 3/2 channels as examples, with the non-local potential [30]: The best values are: V pα, around its best values ±100 MeV and the width parameters β around its best values ±50%. As a representative case, the relative errors for interpolated points in the P 3/2 channel are plotted in Fig. 5 for several basis sizes (additional plots for the S 1/2 channel are given in the SM). The nugget for both channels is set to be 10 −8 . The performance of EC is again excellent except at some isolated energies, and these exceptions are not at the same energies for different basis sizes.
The Kohn variational approach also applies to complex potentials, which are extensively used in optical potentials for nuclear scattering and reactions. To test the EC for these applications, we use a Wood-Saxon optical potential constructed for describing α-208 Pb low-energy scattering [31]: with f (r, R, a) ≡ (1 + exp (r − R)/a) −1 . We take V 0 = −100 MeV, W 0 = −10 MeV, R R = R I = 8.36 fm, and a R = a I = 0.58 fm as the best value [29]. The Coulomb potential is simplified as for point charges [29]. We vary the V 0 and W 0 parameters in different partial waves (so θ i = {V 0 , W 0 }) by ±50% around their best values [29]. In Fig. 6, the size of relative errors for the = 20 channel is shown as a representative example. The results for = 0 are shown in the SM. (Note that the scattering phase shift is complex here. The vertical axis of the plot is for the modulus of the relative error.) The lower end of the energy range is chosen such that the Sommerfeld parameter η is less about 10, because the numerical calculation of Coulomb functions in the R-matrix package becomes unreliable for larger values [29,32]. The upper end of the energy range is chosen to match Ref. [31]. The nugget used in ∆ U 's inversion is set to 10 −10 in both = 20 and = 0 calculations. With 10 basis elements, the relative accuracy for interpolated points is no worse than 10 −4 , while increasing N b further improves it to 10 −5 or better. Again, the errors for interpolated and extrapolated points are similar, and the standard deviations are similar in size to the mean values.
Based on these results for p-α and α-Pb scattering, we expect EC could play an important role in fitting potential parameters for nuclear scattering and coupledchannel reactions.

V. SUMMARY AND OUTLOOK
We have extended the eigenvector continuation method to scattering using the Kohn variational principle. The EC enables accurate calculations of observables for any parameter set θ given calculations of scattering wave functions and K-matrix elements from a limited number N b of parameter sets θ i . Unlike the bound-state application of EC, for scattering the KVP does not give an upper bound to observables but is only guaranteed to give stationary results. Nevertheless, for good trial functions the KVP has been demonstrated in the literature to give accurate results for a wide range of applications [14]. An EC basis provides a very effective trial function and its application to the KVP is simple, involving only the inversion of the matrix defined in Eq. (9). Issues of illconditioning with increasing basis size are successfully treated with simple regularizations.
Here we have provided representative results from a wide range of tests of the EC for scattering using model problems. These include multi-dimensional parameter sets, both local and non-local potentials, charged-particle scattering, and complex optical potentials. In all cases shown here and in all our other tests to date, the EC is found to be effective with moderate basis sizes both for interpolated and extrapolated parameter sets. We are working to formulate a robust uncertainty quantification and to develop a procedure for determining the optimal regularization parameter for ill-conditioning, which has thus far been fixed empirically.
The success of the EC enables the development of efficient emulators for scattering. In subsequent work we will demonstrate the application to coupled channels in both coordinate-space and momentum-space (which is a straightforward generalization of the presentation here) and set up the application to Nd scattering [15,18]. It would be also interesting to apply our method to fit an NN potential to the NN energy spectra from Lattice QCD calculations, since the eigenenergies and phase shifts are directly connected. To have a self-contained presentation, we briefly review the KVP [13,18] for single-channel scattering. The generalization to coupled-channel processes is straightforward [18,19]. For a specified partial wave with local potential V (r), the KVP stationary functional (6) with the trial radial wave function u t (we suppress the dependence on and energy E) reduces to [13] β[u t ] = τ trial − ∞ 0 dr u t (r)Du t (r). (S1) Here, p is the asymptotic momentum (p = √ 2µE) and with U (r) ≡ 2µV (r). The value of τ trial is extracted from the asymptotic behavior of the trial wave function, where the sine term sets the normalization. This normalization convention goes hand-in-hand with the particular form of the β[u t ] functional, which needs to be modified for other normalizations. For the exact radial function u exact , we have and therefore from Eqs. (S3) and (4) with τ trial → τ exact , To see that β is a stationary functional, we write We want to act D to the left to take advantage of Eq. (S4), which requires partially integrating the deriatives in D. This yields only surface terms, to which we can use the asymptotic forms: The lower limit does not contribute because while the upper limit, after using Eq. (S3) for u exact and u t , yields −δτ [sin 2 (pr − 1 2 π) + cos 2 (pr − 1 2 π)] = −δτ (with two other terms canceling), so in the end (S10) Thus, the functional β[u t ] approximates τ exact at its stationary point up to O δu 2 .
If the long-range Coulomb potential is present, the asymptotic behavior of the radial basis functions and of u t differs by the argument in the sine and cosine functions, namely, and similarly for the cos function. Here η is the so-called Sommerfeld parameter and σ the pure Coulomb phase shift [13]. Then the phase shift δ would be the shorterrange interaction-induced phase shift, i.e., the total phase shift with σ subtracted. (Note that the total phase shift here is measured with respect to the incoming and outgoing spherical waves in the form of exp [±i(pr − η ln 2pr)].) For a general non-local potential, the above derivation still holds, except that the integral involving V (r) in Eq. (S1) needs to be changed to a double integration, i.e., drdr u t (r)V (r, r )u t (r ). It should be emphasized that the potential in the function may be complex in general, but in this case no complex conjugation is to be applied to u t in the integral [19]. The applicability of KVP for complex potentials is important for this approach to be used in optical potential model fitting.
The generalization to coupled channels can be found e.g., in Ref. [19,20]. The central step involves finding that the variation of the integral in the definition of the β functional comes from the end points/surface terms of the integral (using Green's theorem), which exactly cancel the variation of the first term in Eq. (S1), which now turns into the K-matrix K [20].

II. ADAPTING EC TO SCATTERING
To adapt EC to the scattering problem, we do not solve an energy eigenvalue problem but instead we find solutions at a specified scattering energy E = p 2 /2µ. We still have a set of N b Hamiltonians specified by {θ i } but now for each one we find the scattering solution for energy E and the corresponding phaseshift. This will be our basis. This identifies the τ i = K (i) /p value in Eq. (4) for each basis wave function. A trial wave function as in Eq. (3), which we write as with each u i normalized according to Eq. (4), will have the asymptotic form because the u i all have the same p. Matching, we find the constraint: So we substitute (S12) into (S1) and require β to be stationary with respect to the c i s, subject to the constraint (S14), which we incorporate using a Lagrange multiplier λ. The functional is To get the second line we have added and subtracted V k (r) inside the integral and used that each of the u j (r) wave functions are eigenstates with their corresponding V j but all the same E. We define (S16) Note that this is not a symmetric matrix. Then the functional to make stationary is: We want to be this stationary under the constraint that the sum of the c j coefficients is one. So (S19) The second equation just gives us back the constraint. The first line, for each i, gives which defines ∆ U ij as in Eq. (9). Finally, in agreement with Eq. (7). We can identify λ by summing this last equation over j and setting it equal to 1, then solving for λ: in agreement with Eq. (8). To get our estimate of τ exact = [K (E)] exact /p as in (10), we substitute from (S23) and (S24) into (S17). Note that the first term on the right side of Eq. (10) is zeroth order in δu while the correction term is first order.

III. ADDITIONAL RESULTS
Here we provide additional results in the partial wave channels not already presented. The setups for these calculations are the same as for their counterparts in the main text. Figures S1 and S2 provide results for NN scattering in the 3 S 1 channel with the Minnesota potential by varying θ i = {V 0R , V 0s } and θ i = {V 0R , V 0s , κ R , κ R }, respectively. Figure S3 shows the results for p-α scattering in the S-wave channel, in parallel to the P-wave results in the main text. Figure S4 plots results for S-wave α-208 Pb scattering. In these plots, both mean values and the standard deviations (std) of the relative errors (absolute values) are plotted against E. As mentioned in the main text and shown in these plots, the std values are similar to the mean values in general. FIG. S1. Relative errors between EC predictions and direct calculations of p/K (E) = p cot θ(E) for the Minnesota potential in the 3 S1 channel. To regularize the ∆ U matrix for its inversion, the nugget is set to 10 −9 here. The top plots show the means of errors for sampled parameter sets in Fig. 2 while the bottom plots show the standard deviations of these errors.