Fast&accurate emulation of two-body scattering observables without wave functions

We combine Newton's variational method with ideas from eigenvector continuation to construct a fast&accurate emulator for two-body scattering observables. The emulator will facilitate the application of rigorous statistical methods for interactions that depend smoothly on a set of free parameters. Our approach begins with a trial $K$ or $T$ matrix constructed from a small number of exact solutions to the Lippmann--Schwinger equation. Subsequent emulation only requires operations on small matrices. We provide several applications to short-range potentials with and without the Coulomb interaction and partial-wave coupling. It is shown that the emulator can accurately extrapolate far from the support of the training data. When used to emulate the neutron-proton cross section with a modern chiral interaction as a function of 26 free parameters, it reproduces the exact calculation with negligible error and provides an over 300x improvement in CPU time.


I. INTRODUCTION
Nuclear scattering experiments yield invaluable data for testing, validating, and improving theoretical models such as chiral effective field theory (EFT) [1][2][3][4]-the method of choice for deriving microscopic nuclear interactions at low energies. However, there are competing formulations of chiral EFT with open questions on issues including EFT power counting, sensitivity to regulator artifacts, and differing predictions for medium-mass atomic nuclei. Low-energy nuclear scattering data combined with rigorous statistical methods such as Bayesian parameter estimation [5], model comparison [6], and sensitivity analysis [7] applied to chiral EFT predictions will provide important insights to address these issues [8].
But taking full advantage of the available data using such statistical methods requires fast & accurate predictions across a wide range of model parameters. While, in principle, the scattering equations can be solved accurately in few-body systems, doing so is prohibitively slow for statistical analyses of three-and higher-body scattering, and even for two-body scattering more efficient alternatives are appealing.
In this Letter, we study such an alternative for twobody scattering that has the potential of future extensions to higher-body systems. We introduce an efficient emulator of the Lippmann-Schwinger (LS) integral equation using Newton's variational method [9,10] combined with ideas from eigenvector continuation (EC) [11,12]. The term emulator refers here to an algorithm capable of approximating the exact solution of a scattering problem with high accuracy while requiring only a fraction of the computational resources.
The power of EC as an emulator stems from the fact that, as the Hamiltonian parameters are varied, the trajectory of each eigenvector remains within a small subspace compared to the full Hilbert space. Linear combinations of eigenvectors spanning this subspace are extremely effective trial wave functions for variational calculations (see also the reduced basis method [13,14]). Emulators based on EC have accurately approximated ground-state properties such binding energies and charge radii, and even transition matrix elements [7,[15][16][17]. Additionally, EC has recently been used to construct effective trial wave functions for applying the Kohn variational principle to emulate two-body scattering observables [18], and for R matrix theory calculations of fusion observables [19]. As we will show in this Letter, Newton's variational method has the feature that scattering observables can be predicted using trial scattering matrices (e.g., the K matrix) rather than trial wave functions. But emulated wave functions can still be obtained [9,20].
The remainder of this work is organized as follows. In Sec. II we briefly describe the formalism underlying the emulator. We then present in Sec. III several applications to short-range potentials with and without the long-range Coulomb interaction and partial-wave coupling. In addition to phase shifts, we study the neutron-proton (np) total cross section by combining multiple emulators across a set of (coupled) partial waves to assess the accuracy and speedup of the emulator in realistic scattering scenarios. We conclude this Letter in Sec. IV, and refer to the Appendices for more technical details including the emulation of gradients required by some Monte Carlo samplers and optimizers. We use natural units in which = c = 1. The self-contained set of data and codes that generates all results shown in this Letter will be made publicly available [21].

II. FORMALISM
We aim to construct an efficient emulator for the LS equation given a short-range potential V ( a) that depends arXiv:2106.15608v2 [nucl-th] 7 Sep 2021 smoothly on a set of parameters a, such as the low-energy couplings of a chiral potential. Specifically, we consider here the LS equation for the scattering K matrix, which reads in operator form 1 with the free-space Green's function operator G 0 (E q ) at the on-shell energy E q = q 2 /2µ and reduced mass µ.
The energy dependence is implicit in what follows. We stress that using the K matrix is just a convenient choice. In fact, T (±) can be emulated by imposing the associated boundary conditions on G 0 . Although the LS equation (1) has the formal solution evaluating Eq. (2) in a given basis can be prohibitively slow for large-scale Monte Carlo sampling because of the fine (quadrature) grids typically necessary to obtain highaccuracy results. Instead of solving the LS equation (1) directly for each sampling vector a, we propose a variational approach starting with a trial K matrix motivated by EC: Here, are a priori unknown coefficients. 2 To determine these coefficients at each a, we apply Newton's variational method [9,10], which provides a stationary approximation to the exact scattering K matrix using the functional given a trial matrix K such as the one in Eq. (3). The functional (4) is stationary about exact solutions of the LS equation, i.e., K[K + δK] = K + (δK) 2 .
In practice, we determine the stationary solution of the functional (4) in a chosen basis and emulate the scattering K matrix as the matrix element φ |K|φ . For example, one could choose |φ to be a plane-wave partial-wave basis |k m with momentum k and angular momentum quanta (l, m), or one could keep the angular dependence explicit via |φ = |k in a single-particle basis. We are interested in emulating K at the on-shell energy E q , so then k = k = q for |φ and φ |. Expressed in the chosen 1 All subsequent equations work for any boundary conditions imposed via G 0 , although we use its principal value formulation here. That is, using G (±) 0 and making the replacement K → T (±) will yield an emulator for T (±) . 2 The coefficients are not normalized, i.e., n t i=1 β i = 1, as opposed to the Kohn variational approach in Ref. [18]. basis, simplifying the functional (4) after inserting (3) yields with If the potential V ( a) is linear in the parameter vector a, then m and M can be efficiently reconstructed by linear combinations of matrices pre-computed during the training phase of the emulator. This results in substantial improvements in CPU time, e.g., for chiral nucleon-nucleon (NN) interactions. By imposing the stationary condition dK/d β = 0, one then finds β ( a) such that M β = m. Given that the optimal β ( a) yields a trial matrix (3) with an error δK, we insert β in Eq. (5) to obtain an error (δK) 2 . The resulting emulator K ( a) ≡ K( a, β ) is then Equations (6)-(8) are the main expressions for emulating scattering observables with short-range interactions. We extend the implementation to the long-range Coulomb potential in Sec. III B. 3 The EC-motivated trial K matrix (3) causes increasingly ill-conditioned matrices M with increasing number of training points n t . To control the numerical noise in the evaluation of Eq. (8), we follow Ref. [18] and add the regularization parameter η = 10 −12 to the diagonal elements M ii . This is a relatively simple yet effective approach compared to other regularization methods [22].
Besides numerical instabilities, Newton's variational method can also exhibit spurious singularities [23], similar to the so-called Kohn (or Schwartz) anomalies [24,25] observed in applications of the Kohn variational principle [26,27]. For instance, we expect spurious singularities to occur at energies where M is singular, i.e., when there is no (unique) stationary approximation to the K matrix due to the functional (4). Different methods to mitigate these singularities have been proposed in the literature [28,29]. Recently, Ref. [30] demonstrated that an EC-driven emulator that assesses the consistency of results obtained from a set of Kohn variational principles (with different boundary conditions) is effective in detecting Kohn anomalies. If detected, they can be mitigated at a given energy, e.g., by changing the number of training points used for emulation. A similar approach could be applied here; however, we have not encountered issues related to spurious singularities in our comprehensive proof-of-principle calculations presented in Sec. III.
In Appendix A, we show that our approach also allows for gradients with respect to model parameters to be straightforwardly propagated. This is an important feature since many optimization and sampling algorithms require gradients. In Appendix B, we discuss a simple and computationally efficient method to evaluate products involving Green's functions in the partial-wave basis.

III. RESULTS
Throughout this section, we use the convention that K = − tan δ (which is opposite to Ref. [18]) with all factors of π, the reduced mass, and momentum accounted for.

A. The Minnesota Potential
Following Ref. [18], we apply our emulator to the Minnesota potential [31] in the 1 S 0 channel as a simple test case: We now demonstrate that the emulator is also a robust tool for extrapolations. We set the Minnesota potential parameters (V 0R , κ R , κ s ) to their best fit values, train on only two purely repulsive parameter sets with V 0s = 30 and 100 MeV, and then extrapolate to purely attractive potentials (capable of supporting bound states). Figure 2 shows the resulting on-shell K matrices in the 1 S 0 channel obtained using our emulator (colored dots) in comparison to the exact solutions of the LS equation (black lines). Each set of colored dots corresponds to a specific center-of-mass energy in the range 1-70 MeV (see the legend for details). The two training points are depicted by the open circles, and the cross marks the location of the best-fit value for V 0s . As the figure illustrates, the emulator can accurately extrapolate far away from the two training points even after passing through poles in both K and K −1 .

B. Including the Coulomb Interaction
Long-range interactions, such as the Coulomb interaction, are problematic for the LS equation approach whether or not an emulator is employed. Nevertheless, we can include the Coulomb interaction via the Vincent-Phatak method [32,33]. The basic idea is to cut off the Coulomb potential at a finite radius so that Eq. (1) applies and then restore this physics using a matching procedure. Specifically, we emulate the K matrix from the potential V rc (r, r ) = V s (r, r )+V rc C (r)δ(r −r )/(rr ), where V s (r, r ) is a (non-local) short-range potential and is the Coulomb potential cut off at a radius r c large enough such that the short-range potential V s is negligible. The modified potential V rc (r, r ) is short-ranged and hence compatible with Eqs. (6)- (8); this is the potential we use to train the emulator. Suppose we want to emulate the K matrix in an uncoupled partial-wave channel with angular momentum . Solving Eq. (8) with V rc yields the associated K rc , but this is an artificial quantity representing the phase shift relative to the free radial wave functions, i.e., u rc (q, r) = j (qr) + K rc y (qr), for r r c , expressed in terms of the regular j (qr) and irregular y (qr) Riccati-Bessel function.
To obtain the phase shifts with respect to the Coulomb wave functions (Sommerfeld-parameter dependencies being implicit), i.e., u (q, r) = J (qr) + K C Y (qr), for R r r c , (12) we match the logarithmic derivatives of Eqs. (11) and (12) at r = r c . (Here R is the range of the shortrange potential.) This amounts to computing (13) and primes denote derivatives with respect to r. Now, both K C = − tan δ C and the phase shift δ C are with respect to the Coulomb wave functions such that the dependence on the choice of r c has been removed. Note that the above relies on a sign convention where, e.g., y 0 ∼ − cos(qr) and similarly for Y 0 . Each of j , y , J , Y and their derivatives can be computed once and stored for emulation purposes. Solving for K C need only be performed for the on-shell matrix and is a quick postprocessing step for the emulator. We apply this approach, with r c = 20 fm, to proton-α scattering with the non-local potential [34] in the S-wave; i.e., = 0. With the four training points V pα, = −6.5 fm −3 [34], as shown in Fig. 3. Across the energy range shown in the figure, E cm 60 MeV, the absolute residuals in the phase shift emulator (bottom panel) are negligible. The high accuracy obtained is remarkable because computing the phase shift at each energy only involves inverting a 4×4 matrix.

C. Coupled Channels
The straightforward extension to scattering in coupled channels is one of the strengths of our emulator approach. In fact, Eqs. (6)-(8) handle them as a special case. Although the term coupled channels can also refer to different reaction channels, in the following we specifically consider coupled (spin-triplet) partial-wave channels.
For potentials that are coupled across n c different channels (e.g., n c = 2 for the deuteron), solving the LS equation exactly for one on-shell point requires solving a linear system of dimension n c n k × n c n k , with n k being the size of the mesh for an uncoupled channel. The coupled-channel emulator with n t training points instead only involves operations on an n t ×n t matrix for each desired matrix element of K, where n t n k . Generally, this requires running at most n c (n c + 1)/2 of such emulations because the remaining matrix elements can be determined by symmetry.
We apply this approach to np scattering in the coupled 3 S 1 -3 D 1 channel. The potential used here is the semilocal momentum-space (SMS) regularized chiral potential at N 4 LO+ constructed by Reinert, Krebs, and Epelbaum with momentum cutoff Λ = 450 MeV [35]. At this chiral order the 3 S 1 -3 D 1 channel depends on n a = 6 nonredundant parameters, or low-energy constants (LECs), in the NN sector [35]. We choose n t = 2n a = 12 training points randomly in the range [−5, 5], where the unit of each parameter is as given in Ref. [35] and left implicit here. The emulator's predictions are then validated at the best values of the parameters found in Ref. [35]. We use a compound Gauss-Legendre quadrature mesh of 80 momentum points to exactly solve the LS equation at the training points. Figure 4 shows the resulting on-shell K matrix obtained from the exact calculation and emulator as a function of the laboratory energy. Each column corresponds to a different partial-wave component of the K matrix. The emulator accurately reproduces the exact K matrix elements across the wide range of energies shown, 1-350 MeV. Except for a spike near the energy region where the K matrix is singular, the residuals are on the order of 10 −12 . These errors are far beneath the experimental uncertainties if the K matrix were to be converted to phase shifts [36].

D. The Scattering Cross Section
We now combine multiple partial-wave emulators into an overall emulator for nuclear observables. As a simple example, we show np total cross sections using partial waves up to j max = 20-again with the Λ = 450 MeV N 4 LO+ SMS potential [35], which reproduces well the total cross sections from the partial-wave analysis [36] over a wide range of laboratory energies. This requires training partial-wave emulators across singlet and triplet channels up to j = 4, while the remaining waves are fixed with respect to a. There are a total of 26 free parameters in a. The n t training locations are again chosen randomly in [−5, 5], where n t is determined based on the n a NN LECs in each partial wave via n t = max(2n a , 4).
Upon emulating K j , the total cross section can be calculated via where S j = 1 + 2i(1 + iK j ) −1 K j and q is the centerof-mass momentum. Both S j and K j are 4 × 4 matrices that contain both the triplet-triplet and the singlettriplet channels. Figure 5 shows the emulated σ tot at the optimal values of a determined in Ref. [35]. The emulator has an error 10 −10 mb for E lab > 50 MeV. For E lab < 50 MeV, the error is 10 −8 mb except for the spike due to the singular K matrix in the 3 S 1 -3 D 1 channel, as discussed in Sec. III C. In either case, these errors are vanishingly small compared to both the size of the cross section itself and its experimental uncertainty [36].
When randomly sampling 500 values of the NN LECs in the range of [−15, 15]-an extrapolation of ±10 beyond the range of the training data (in the appropriate units)-the average absolute emulator error is less than FIG. 5. The np total cross section using the Λ = 450 MeV N 4 LO+ SMS potential [35] up to jmax = 20. The inset shows the mean absolute error between the emulator and the exact solution for 500 different samples of the NN LECs a (red line) as well as the absolute residual at the best-fit a. These errors are negligible compared to any experimental uncertainties. See the main text for more details. 10 −7 mb. Furthermore, the emulator provides a factor of > 300x improvement in terms of CPU time relative to the exact calculation. If the size of the momentum mesh used in the LS equation is increased from 80 to 160 quadrature points, then the factor becomes > 1000x. Further acceleration can be expected as finer momentum meshes are used in solving the LS equation, and as a enter into higher partial waves at higher chiral orders.
Meaningful comparisons of the speed and accuracy obtained here with the emulators in Refs. [18,30] requires, at least, that all calculations be performed in the same space. For this benchmark, we have implemented the Kohn variational principle with uncoupled channels in momentum space. Our findings so far indicate that the two variational methods are comparable in accuracy (for the same quadrature rule) for the chiral potential, although the relative speedups are implementation dependent. More work is necessary to provide more quantitative comparisons.
We also compare the accuracy of our emulator to a promising accelerator for NN scattering observables developed in Ref. [37]. Instead of a variational method, Miller et al. employed the wave-packet continuum discretization (WPCD) method to approximate scattering solutions at multiple energies at once. This method is well-suited for parallelization using Graphics Processing Units, which can lead to significant speedups compared to exact calculations via conventional matrix inversion. Depending on the laboratory energy and the number of wave packets included, Miller et al. reported averaged errors in the total cross section on the order of 1 mb at best based on the chiral interaction NNLO opt [38]. These errors suggest that our emulator motivated by EC can provide significantly higher accuracies, even when only a few training points per partial-wave channel are used. A quantitative comparison of the methods' efficiencies, however, would require a scattering scenario with matching nuclear interactions.

IV. SUMMARY AND OUTLOOK
We showed that Newton's variational method combined with ideas from eigenvector continuation allows for the construction of a fast & accurate emulator for twobody scattering observables. Our approach begins with a trial K or T matrix constructed from a small number of exact solutions to the LS equation in the parameter space of the Hamiltonian. Subsequent emulation only requires linear algebra operations on low-dimensional matrices.
We then provided several applications to short-range potentials with and without the Coulomb interaction and partial-wave coupling. In all cases studied, the emulator is capable of reproducing phase shifts and total cross sections with remarkable accuracy, even far from the support of the training data and across poles in K and K −1 . In particular, for a modern chiral interaction at N 4 LO+ the emulator reproduced the exact neutron-proton cross section with negligible error but was over 300x faster in CPU time. The code that generates all results and figures within this Letter will be made publicly available [21].
While the number of emulators applicable to boundstate observables in few-and many-body systems is growing [7,15,16,39], developing methods with similar efficacy for three-and higher-body scattering is an important avenue. Thanks to emulators, simultaneous Bayesian fits of chiral interactions to pion-nucleon, nucleon-nucleon, and three-nucleon observables [16] with theoretical uncertainties rigorously quantified [5,40] already have become feasible. Next-generation emulators have the potential to extend these studies to three-and higher-body scattering observables and to shed light on important issues inherent in chiral EFT. Our approach, together with the advances made in applying Kohn variational principles based on EC trial wave functions to three-body scattering [41], is promising in this direction. Further, the fast convergence we observed with EC-inspired trial matrices (instead of wave functions) motivates the exploration of the EC concept applied to stationary functionals in a more general context. Altogether, these are exciting prospects for rigorous Bayesian uncertainty quantification in nuclear physics and reaction theory.
The value of β = M −1 m must already be computed for the emulator itself and thus can be reused in Eq. (A1). The M matrix is symmetric, which means that m M −1 = (M −1 m) requires no further computation. If V ( a) is linear in a, then each projection of the gradient tensors ∂ m/∂a k and ∂M/∂a k can be performed once and stored. Therefore, all components of the gradient at a can either be pre-computed during training or have already been completed during the emulation step for K( a)-only matrix multiplication remains to be done to obtain ∂K/∂a k . Thus, not only is the emulator fast to compute, the gradient of the emulator is also fast because it operates in the small space of training points. This makes gradients feasible to incorporate into sampling codes with little computational overhead. As an example, Fig. 6 shows the gradients of the K matrix from the Minnesota potential considered in Sec. III A-with the same training points. The gradient is computed at the best values of V 0R and V 0s ; the residuals compared to the exact calculation are negligible for all center-of-mass energies shown.

Appendix B: A Convenient Form for the Partial-Wave Green's Function
The free-space Green's function regularly acts between operators when setting up the emulator, and possibly during each emulation if V ( a) cannot be projected and stored up front. Although a product like V G 0 K appears straightforward to evaluate, there are technicalities involving numerical instabilities and integration measures that are obscured when the LS equation (1) and the emulator equations (6) and (7) are written in operator form. Thus, it can be convenient to construct a form of G 0 that can be applied as a matrix-matrix product, such that all the aforementioned equations can be evaluated straightforwardly. This approach has the added benefit that it only requires generating the potential V on the fixed grid used for solving the LS equation [43,44], rather than appending entries for each specific on-shell solution [45]. This means that generating V and K for new parameters a becomes more efficient in both runtime and memory.
When solving the LS equation in partial waves, the following projection arises [45]: where we work in uncoupled channels for simplicity and µ is the reduced mass. The P denotes the principal value integral due to our choice of G 0 . To avoid the numerical instability (i.e., pole) at k = q, a zero integral is subtracted to yield p m|V G 0 K|p m (B2) = 2 π ∞ 0 dk V (p , k)K (k, p; q)k 2 − V (p , q)K (q, p; q)q 2 (q 2 − k 2 )/2µ .
See also Ref. [46] for a similar numerical approach. This can be compressed by writing p m|V G 0 K|p m = 2 π V (p , k) dG 0 (k; q) K (k, p; q), where, in the partial-wave basis, The left-hand term is the standard free-space Green's function with a factor of k 2 dk included. The right-hand term is the zero integral (in principal value) included for numerical stability, which has been multiplied by a factor that will get the on-shell portion of whatever matrices it acts upon. In practice, the potential-and hence K-must be evaluated on a grid using, e.g., Gauss-Legendre quadrature for k and dk. This means that the dk δ(k − q) will not be able to set the gridded matrix elements to k = q. Thus, we replace dk δ(k − q) with an interpolation vector S(q) that performs the mapping for any smooth f that has been evaluated on some grid of k, as in Ref. [44]. Because both the quadrature grid in k and the required on-shell locations q are fixed throughout the emulation process, S k (q) needs only to be computed once. Therefore, all the components of dG 0 are independent of a and we can avoid unnecessary calculations while sampling.
The resulting matrix dG 0 , which is diagonal in momentum space, is what we use as G 0 in all emulators shown here. It reduces the radial integrals over k and singularity smoothing in products like V G 0 K to a matrix product of V dG 0 K for matrices on a fixed mesh for k and dk. Another benefit of this approach is the simplicity it brings to the exact solutions of the LS equation: we can now use dG 0 in Eq. (2). Upon solving for K (p , p; q) via Eq. (2), one can again use the interpolation vector S(q) to compute the on-shell component: K (q) = S(q) K S(q). This sidesteps the requirement of creating unique (n k + 1) × (n k + 1) matrices for each desired on-shell K (q), as espoused in Ref. [45] and elsewhere. The advantage is prominent when computing the different partial waves for the total cross section calculation.