Toward emulating nuclear reactions using eigenvector continuation

We construct an efficient emulator for two-body scattering observables using the general (complex) Kohn variational principle and trial wave functions derived from eigenvector continuation. The emulator simultaneously evaluates an array of Kohn variational principles associated with different boundary conditions, which allows for the detection and removal of spurious singularities known as Kohn anomalies. When applied to the $K$-matrix only, our emulator resembles the one constructed by Furnstahl et al. [Phys. Lett. B 809, 135719] although with reduced numerical noise. After a few applications to real potentials, we emulate differential cross sections for $^{40}$Ca$(n,n)$ scattering based on a realistic optical potential and quantify the model uncertainties using Bayesian methods. These calculations serve as a proof of principle for future studies aimed at improving optical models.


I. INTRODUCTION
There are many reasons to study rare isotopes today; e.g., they play a crucial role in obtaining a fundamental understanding of the nucleosynthesis of heavy, neutronrich nuclei and the dense matter inside neutron stars [1][2][3]. Due to their short lifetimes, rare isotopes are primarily investigated through reaction experiments conducted at radioactive beam facilities worldwide, including RIKEN, FAIR, GANIL, and soon also FRIB. For the analysis and interpretation of these experiments, reliable reaction theory is imperative. However, apart from reactions on light nuclei, reaction theory is still largely phenomenological and relies on poorly constrained effective interactions to keep calculations tractable [4].
Statistical methods such as Bayesian parameter estimation [5] and model comparison [6] can provide important insights into the issues of effective interactions. They can also help design next-generation reaction experiments (see, e.g., Ref. [7]). But in practice their applications are limited because Monte Carlo sampling of the models' parameter spaces in reaction calculations is usually computationally demanding. Hence, Bayesian studies of nuclear reactions [8][9][10] have only considered the simplest reaction theory, the optical model, which describes, e.g., nucleus-nucleus scattering as two particles interacting via a complex-valued interaction. Extending these studies to more sophisticated reaction theories (see, e.g., Refs. [11][12][13][14][15][16][17][18][19]) is a challenging yet important task.
Emulators-computationally inexpensive algorithms capable of approximating exact model calculations with high accuracy-are promising tools in this regard [20][21][22][23]. In particular, eigenvector continuation (EC) [24,25]  has been shown to be a powerful method for emulating bound-state properties such as binding energies and charge radii of atomic nuclei [26][27][28]. Furnstahl et al. [29] have recently demonstrated that EC also allows for the construction of effective trial wave functions for calculations of two-body scattering observables using the Kmatrix Kohn variational principle (KVP) [30]. Further, Melendez et al. [31] have extended the EC concept to trial K-or T -matrices in applications of Newton's variational method to two-body scattering, e.g., with a modern chiral interaction. Remarkably high accuracies and speedups relative to exact scattering calculations were obtained [29,31]. (See Ref. [32] for EC applied to Rmatrix theory calculations of fusion observables.) In this article, we improve and extend the emulator developed by Furnstahl et al. [29] in several ways. Besides the K-matrix, we emulate a variety of matrices associated with different scattering boundary conditions simultaneously via the general (complex) KVP [33]. (For pre-EC studies with this method, including nucleon-deuteron scattering, see Refs. [34][35][36].) This approach allows us to detect and remove spurious singularities known as Kohn anomalies [37,38], which can render variational calculations of scattering observables ineffective-especially when used for sampling a model's parameter space. We also propose a method for solving the emulator equations in Ref. [29] with reduced numerical noise. As a step toward emulating nuclear reactions, we apply our emulator to differential cross sections in 40 Ca(n, n) scattering using a realistic optical potential and quantify the uncertainties in the model parameters using Bayesian methods.
The remainder of this article is organized as follows. In Section II we introduce the formalism of the general KVP with EC trial wave functions. We then present several applications of our emulator to realistic potentials (including a chiral potential) in Section III. Section IV concludes the article with a summary and outlook. Additional information, e.g., on redundancies among different KVPs with EC trial wave functions, is provided in Appendices A to F. We use natural units in which = c = 1.

II. FORMALISM
We consider here local short-range potentials V (θ) in coordinate space that depend on a set of free parameters θ; e.g., the parameters of an optical model or lowenergy couplings of a chiral potential. Further, V (θ) is assumed to be partial-wave decomposed into an uncoupled channel with angular momentum . Following Furnstahl et al. [29], we then use EC to construct an effective trial wave function for our (nonrelativistic) variational calculations of two-body scattering observables: Here, each of the N b basis wave functions, i.e., is an exact (partial-wave) solution to the Schrödinger equation for V (θ i ) at the center-of-mass energy E > 0, and the coefficients c ,E are to be determined. The radial wave functions in Eq. (2) are normalized by imposing asymptotic boundary conditions 1 of the general form [33] φ ,E (r) ∼φ where the two independent free-space solutions are expressed in terms of a nonsingular (complex) matrix u that is associated with the generic L-matrix in Eq. (3): with η (r) = pr − π 2 and p = √ 2µE, and an arbitrary normalization constant N = 0. For instance, the familiar K-, S-, and T -matrix respectively correspond to 2 (5) But any other nonsingular parametrization (L, u) of the asymptotic limit (3) is equally valid. 3 The corresponding K-matrix can be obtained using the Möbius transformation (see also Appendix A) 1 The boundary condition (3) can be extended to the (long-range) Coulomb potential. See, e.g., Eq. (S11) in Ref. [29]. 2 The matrices are determined only up to scalar multiples. Lucchese [33] uses a different convention for the S-matrix (u → −u) and T -matrix parametrization (T → −πT ). which is related to the phase shift via K ,E = tan δ ,E . In Appendix B we give the technical details for solving the (radial) Schrödinger equation numerically to determine the basis wave functions in Eq. (1), including a generalization of Eq. (6) to transform from a given Lmatrix to any other L -matrix-not just the K-matrix.
We determine the coefficients c (i) ,E in Eq. (1) using the general (complex) KVP. Given (L, u) and a trial wave function r|ψ trial subject to the boundary condition (3), the general KVP provides a stationary approximation tõ L ,E = L ,E /N using the functional [33] with the reduced mass µ ≈ A p A t /(A p + A t ) m n , mass number of the projectile A p = 1 (here, a neutron with mass m n ) and target A t , respectively, as well as the Hamiltonian H(θ) = −∇ 2 /(2µ) + V (θ) in coordinate space. A derivation similar to the one in Refs. [29,39] for the K-matrix shows that the functional (7) is indeed stationary about exact solutions to the Schrödinger equation, i.e., β u [|ψ ,E +|δψ ,E ] =L ,E +(δL ,E ) 2 , although it does not provide an upper or lower bound in general.
We impose the normalization constraint i c (i) ,E = 1 on the EC trial wave function (1) to fulfill the boundary condition (3) required by the general KVP. Constrained optimization of the functional (7) using a Lagrange multiplier λ then leads to the system of (N b + 1) linear equations and eventually to the desired stationary solution [29] c (i) In Eq. ,E . Further, we have defined the which can be efficiently evaluated for a variety of different (L, u) at once, as discussed in Appendix C. Hence, the stationary approximation to L ,E reads ,E .
(11) Equations (8) to (11) are the main expressions of our emulator. They reduce to the ones derived in Ref. [29] for the K-matrix KVP (i.e., L = K) with N = p. Hence, the discussion of the computational complexity in Ref. [29] also applies to our emulator for each (L, u).
The EC trial wave function (1) renders the kernel matrix ∆ U ,E . Although this simple approach typically works well in practice, we find that solving instead the system of equations (8) numerically using a least-squares solver [to circumvent explicit matrix inversion] is less sensitive to numerical noise-especially at low energies. The least-squares solver uses a different regularization method, where singular values less than a given cutoff ratio times the largest singular value are considered zero. In the cases we studied, this cutoff ratio could be as small as the machine epsilon to avoid potential fine-tuning.
In addition to these numerical instabilities, the general KVP is prone to spurious singularities known as Kohn (or Schwartz) anomalies [37,40], which occur at energies where the functional (7) does not provide a (unique) stationary approximation (11). (See Section III A for several illustrations.) For the realistic potentials studied here, we find that neither real KVPs, such as the one for the K-matrix, nor complex KVPs, such as the one for the Smatrix [41,42], can guarantee anomaly-free results [33].
We therefore emulate a wide range of matrices associated with different scattering boundary conditions simultaneously using the general KVP and assess their consistency. As pointed out in Appendix D, however, not all KVPs (with EC trial wave functions) provide independent stationary approximations-we derive a simple condition to identify those. Results that do not pass the consistency checks, e.g., SS −1 = 1 [33,43], are disregarded by our anomaly detection algorithm and the remaining ones averaged over in an attempt to obtain anomaly-free results. If none of the KVPs evaluated are consistent, our algorithm iteratively adapts the size of the training set, which usually shifts the Kohn anomalies in each iteration. We refer to this approach as the "mixed approach." More details on detecting and removing Kohn anomalies are presented in Appendix E.

A. Realistic Real Potentials
We apply our emulator first to three real potentials as test cases. Specifically, we consider nucleon-nucleon (NN) scattering in the 1 S 0 channel 5 based on the Minnesota potential [44], and the local chiral potential at next-to-next-to-leading order (N 2 LO) developed by Gezerlis et al. [45] with regulator cutoff R 0 = 1.0 fm and spectral-function cutof Λ = 1000 MeV. The Minnesota potential allows for direct comparisons with the emulators constructed in Refs. [29,31], and the chiral potential is commonly used in quantum Monte Carlo calculations of atomic nuclei and nuclear matter (see, e.g., Refs. [46,47] for recent reviews). The latter potential depends on 8 parameters (i.e., NN low-energy couplings) in the 1 S 0 channel. 6 Both were constructed to reproduce 1 S 0 scattering phase shifts. We also consider the scattering states of n+ 10 Be based on the real Woods-Saxon potential with the spin-orbit (LS) term added, i.e., with the function which was fit in Ref. [48] to low-lying states in 11 Be, including the d 5/2 resonance. 7 Equation (13) is commonly used to describe the interaction of the valence nucleon(s) with the core nucleus in halo nuclei, such as 11 Be(n+ 10 Be) in reaction models [11,13,48], or 16 Be(n + n+ 14 Be) in decay studies [17,49]. We consider here the d 5/2 channel (rather than s 1/2 ) because the breakup calculations in Ref. [48] identified this channel as the dominant one for this scattering process (see Figure 1 in Ref. [48]). For the Minnesota potential (12), we follow Furnstahl et al. [29] and train our emulator on the set of points (V 0R , V 0s ) = {(0., −291.85), (100., 8.15), (300., −191.85), 5 The spectroscopic notation 1 S 0 indicates that the angular momentum = 0 ("S") and the total spin S = 0 of the two nucleons couple to the total angular momentum J = 0. 6 Only two independent (spectroscopic) low-energy couplings contribute to the 1 S 0 channel, which are given by linear combinations of the couplings mentioned in the text. For details see, e.g., Appendix A in Ref. [45]. 7 The spectroscopic notation d 5/2 indicates that the angular momentum = 2 ("d") is coupled to a total angular momentum of the valence particle j = 5/2.
(300., 8.15)} in units of MeV, while the other (nonlinear) parameters are fixed at their best fit values, i.e., κ R = 1.487 fm −2 and κ s = 0.465 fm −2 [44]. For the other two potentials, we randomly select the training points within a ±20% interval (in the appropriate units) of the parameters' best fit values, as given in Table I of Ref. [45] for the chiral potential (N b = 4) and Table I of Ref. [48] for the Woods-Saxon potential (13) with fixed . In all cases, we emulate the scattering phase shifts at the best fit values. Figure 1 shows the emulated phase shifts (a-c) and their absolute residual (d-f) relative to the exact scattering solution as a function of the center-of-mass energy. From left to right, the columns correspond to the results obtained for the Minnesota, chiral, and Woods-Saxon potential, respectively. Each panel depicts the emulated results based on the KVPs for the K-, T -, and T −1 -matrix, as well as our mixed approach as solid lines. The KVPs for the other canonical matrices (i.e., K −1 , S, and S −1 ) do not provide complementary stationary solutions (as discussed in Appendix D) and therefore are not shown.
Overall, our emulator reproduces well the exact phase shifts. The absolute residuals typically are 0.01 • , except for the Minnesota potential at the low energies where the phase shift is large. As expected, the K-matrix KVP (orange lines) reproduces the phase shifts obtained by Furnstahl et al. [50] for the Minnesota potential, including the noticeable Kohn anomaly at E ≈ 13 MeV. The T −1 -matrix KVP is anomalous at E ≈ 59 MeV. In the energy range shown, we also find such an anomaly for the chiral interaction at E ≈ 61 MeV, and for the Woods-Saxon potential at E ≈ 8 MeV. Additional Kohn anomalies, however, may be present and only noticeable when using extremely fine energy grids [29]. Figure 1 emphasizes the need for efficient anomaly removal algorithms beyond proof-of-principle calculations, where the exact scattering solution as a reference is not available.
Such an algorithm is implemented in our emulator (see Section II). Depicted by the red lines in Fig. 1, the mixed approach is capable of detecting and removing Kohn anomalies by assessing the consistency of the results obtained from a set of different KVPs and (if necessary) adaptively removing basis wave functions from the training set used for emulation. In this specific case, we have simultaneously emulated the complementary matrices L = (K, T, T −1 ), as shown in the figure, as well as the three additional matrices specified in Appendix E. No changes in the training set were necessary to mitigate these Kohn anomalies.

B. Realistic Optical Potential
We also apply our emulator to a realistic optical potential for 40 Ca(n, n) scattering at E = 20 MeV in the center-of-mass frame. Parametrizations of optical potentials (see, e.g., Ref. [51]) typically contain real and imag-inary terms of the Woods-Saxon form: We do not consider the spin-orbit term in Eq. (15) and assume in the following R w = R v and a w = a v , as in Ref. [51]. To train the emulator, we randomly select N b points for the remaining seven parameters, i.e., within a ±20% interval (in the appropriate units) centered around the Koning-Delaroche (KD) parameterization [51] at E = 20 MeV. This approach allows us to probe a realistic region of the parameter space. Figure 2 shows (a) the emulated differential cross sections (mixed approach) and (b) their corresponding average relative residuals as a function of the scattering angle θ-which is not to be confused with the parameter set of the interaction, V (θ). The exact scattering solutions serve as the reference for the residuals and their mean value is depicted by the black-dotted line in panel (a). We emulate the differential cross section at 500 randomly selected points in the parameter space similar to the training phase, and determine the bands shown in panel (b) as the range spanned by the 50% limit (i.e., median) and (upper) 95% limit of the residuals. The solid lines in both panels correspond to the average results for the emulators with N b = 4 (red lines), N b = 6 (orange lines), N b = 8 (green lines), and N b = 10 (blue lines), respectively. We include partial-wave channels with angular momentum 10 in the calculations. As shown in Fig. 2, the accuracy of the emulator roughly improves by an order of magnitude when increasing the size of the training set by increments of two, from N b = 4 to 10. But the accuracy can also vary by more than an order of magnitude within the 500 sampled points. Furthermore, increasing the scattering angle tends to decrease the accuracy, which is lowest at the backward angles where the differential cross section is smallest. Nevertheless, for N b 6, the emulator residual does not exceed the experimental uncertainty, typically of the order of ≈ 10% (see, e.g., Refs. [52,53]).

C. Uncertainty Quantification for Optical Models
In this section we explore Bayesian parameter estimation and uncertainty quantification of an optical model using our emulator-as a step toward systematic studies in the future. For the proof-of-principle calculation we consider again 40 Ca(n, n) scattering at 20 MeV 8 and use the mixed approach with N b = 8 training points. The  [45] at N 2 LO (center column), and Woods-Saxon potential with spin-orbit term (13) (right column). Both, the Minnesota and chiral potential, are used for NN scattering, whereas the Woods-Saxon potential is used for n+ 10 Be scattering. The dotted vertical lines highlight the (approximate) locations of the detected Kohn anomalies. Notice that the algorithm proposed here (red lines) is capable of removing these anomalies. See the main text for more details.
real and imaginary volume depths and radii of the optical potential (i.e., V v , R v , W d , and R d ) are constrained based on mock data generated from the KD potential [51] (see Ref. [54] for more details), whereas the other optical model parameters are fixed at the original KD values. Each parameter's prior is taken to be a normal distribution with mean set to the KD potential value and width of 50% of the mean, similar to previous studies [8][9][10], and the likelihood is the standard exponentiated χ 2 . The uncertainty quantification is performed through Markov Chain Monte Carlo (MCMC) sampling with 20,000 accepted parameters sets from a single Markov chain. We also obtain 95% confidence intervals for the differential cross sections, defined as the smallest interval over which the posterior distribution integrates to 0.95. Figure 3 shows the results of the parameter estimation based on the mixed approach (red lines) and the exact scattering solution (black lines). Panel (a) gives the posterior distributions for the four varied parameters along the diagonal, with contour plots displaying the correlations between each pair of parameters in the off-diagonal panels (also known as corner plot). Panel (b) compares the resulting 95% confidence intervals for emulated vs. exact differential cross sections. Apart from statistical fluctuations, the emulator reproduces well the exact calculations of the parameter posterior distributions, correlations, and the confidence intervals for the angular distributions. Remarkably, our mixed approach obtained anomaly-free results without adapting the training set in all of our MCMC runs.
The mean values of the posterior distributions match the KD parameters, as expected, and the uncertainties of the parameters and the differential cross sections are similar to what has been obtained in previous studies [8,54]. Note that the same reaction has been studied in Ref. [54] at slightly lower energy but with a larger set of parameters allowed to vary in the MCMC sampling.

IV. SUMMARY AND OUTLOOK
Motivated by the recent Letter by Furnstahl et al. [29], we constructed an efficient emulator for two-body scattering observables using the general KVP [33] and trial wave functions derived from EC. Our emulator does not only consider the K-matrix KVP (as in Ref. [29]), but rather simultaneously evaluates an array of KVPs associated with different (complex) scattering boundary conditions. This approach allows us to systematically detect and remove spurious singularities known as Kohn anomalies, which can render applications of the KVP and other variational principles ineffective; especially when used for Monte Carlo sampling of a model's parameter space. If only the K-matrix KVP is evaluated, our emulator resembles the one constructed in Ref. [29] although with reduced numerical noise. We demonstrated the efficacy of the proposed algorithm for removing Kohn anomalies by emulating scattering phase shifts obtained from the Minnesota, a local chiral, and the real Woods-Saxon potential. For each potential, we found anomalies in at least one of the applied KVPs, which the algorithm reliably detected and removed-without adapting the size of the training set. This emphasizes that Kohn anomalies need to be dealt with in practice, even in proof-of-principle calculations, but doing so does not require the exact scattering solution. The basic concept of the algorithm is general and might also be applicable to other variational methods [31]. Furthermore, we showed that, although the emulator's rate of convergence can be sensitive to the details of the interaction and the size of the training set, the high accuracies obtained with our KVP-based emulator are well-suited for scattering calculations.
After these test applications to real potentials, we studied the EC convergence for emulating differential cross sections in 40 Ca(n, n) scattering at 20 MeV using the realistic KD optical potential. A training set with N b = 6−10 wave functions, typically, was enough to obtain highaccuracy results for this observable. Next, we performed Bayesian parameter estimation for the optical model by optimizing the emulated differential cross section to reproduce mock data calculated from the KD potential. The sampled distribution functions for the model parameters and the differential cross section obtained with the emulator were in excellent agreement with those calculated from the exact scattering solution.
Important future avenues include the extension of our emulator to scattering in coupled partial-wave and reaction channels, with coordinate and momentum space interactions, as well as the inclusion of the (long-range) Coulomb interaction [29,31,39]. Technically more challenging will be the extension to emulating three-and higher-body scattering observables, where the computational efficiency of emulators is vital for rigorous uncertainty quantification. Recent developments in this direction [55,56], however, are promising and will benefit from the insights into the EC-driven general KVP provided here. As the number of efficient emulators for scattering observables increases [29,31,57], it will be important to benchmark the different emulators quantitatively, e.g., in terms of accuracy, computational speedup, and susceptibility to anomalous behavior. These advances set the stage for constructing next-generation optical models using emulators for scattering observables in the FRIB era. The Möbius (or linear fractional) transformation refers to the function (for more details, see, e.g., Ref. [59]) L a (z) = a 01 + a 11 z a 00 + a 10 z , generated by the nonsingular 2 × 2 matrix a. We have chosen the order in which the coefficients a ij appear such that Eq. (6) reads K(L) ≡ L u (L). If a was singular (i.e., det a = 0), then Eq. (A1) would be just a constant, if a 10 = 0 , a01 a00 if a 10 = 0 and a 00 = 0 , undefined if a 00 = a 10 = 0 , (A2) and thus not strictly considered a Möbius transformation. L a (z) has the properties L a (z) = L λa (z) , with λ = 0 , L a (z → ∞) = a 11 a 10 , if the limit exists, (A3b) Further, it can be efficiently implemented using (mostly) linear algebra operations; e.g., (A4) Note that the vector representation of a fraction is only determined up to an arbitrary factor λ = 0. In this work, we use the Möbius transformation to relate different asymptotic limit parametrizations with one another. We write the radial Schrödinger equation for a given angular momentum and center-of-mass energy E as a system of coupled first-order differential equations, and numerically solve it for each of the N b partial-wave decomposed potentials V (r; θ 1 ), . . ., and V (r; θ N b ) using the explicit Runge-Kutta method in Scipy's integrate.solve ivp(). The relative and absolute tolerance are each set to 10 −9 or less. As initial values for the solver we set φ(ε) = 0 and (by choice) φ (ε) = 1, where the value of the derivative will be rescaled later on by imposing an asymptotic boundary condition, and ε > 0 is a numerical value close to zero. We solve the radial Schrödinger equation up to the matching radius r m ∞ located outside the range of the potential. At r m , we smoothly match the numerical solution to the free-space solution parametrized by φ Here, j (pr) and η (pr) denote the spherical Bessel function and Neumann function, respectively. Notice that the asymptotic limit of the free-space solution (B2) is defined in Eq. (4). The (arbitrary) constant N −1 = p is chosen following Ref. [29]. Given a parametrization u, we determine the value of the L-matrix in terms of the inverse logarithmic derivative with respect to r, R(r m ) = φ(r m )/φ (r m ), as follows and then rescale φ(r) by the factor φ (free) ,E (r m )/φ(r m ). More details can be found, e.g., in Ref. [60].
The numerical solution matched to any asymptotic boundary condition of the form (3) is an equally valid solution of the radial Schrödinger equation for r ε. In practice, we choose a particular boundary condition (e.g., with L = S) for the matching. To efficiently transform wave functions normalized by this asymptotic limit parametrization (L, u) to another (L , u ), we use the analytic expressions derived in the following. Notice that primes [e.g., as in φ 0 (r)] no longer indicate derivatives.
We consider the identity in the asymptotic limit  Table 3.1 of Ref. [60].
These integrals can be evaluated to a high accuracy using Gauss-Legendre quadrature rules distributed across multiple intervals. The matrix only depends on the interactions used for training and thus needs to be evaluated only once (for a given u), whereas A ij = I ij [V (r; θ)] has to be evaluated each time the emulator is invoked after the training phase. Our emulator applies a set of KVPs with different boundary conditions. Instead of constructing the kernel matrix for each KVP individually, we make use of the analytic transform for the wave functions derived in Appendix B to relate two kernel matrices associated with (L, u) and (L , u ), respectively. This amounts to the element-wise (i.e., Hadamard) matrix product: where the subscripts index the basis wave functions used for training. While the first two factors on the right-hand side of Eq. (C2) transform the wave functions in the integrals of A ij and B ij , as discussed in Appendix B, the third factor from the left corrects for the different determinants in Eq. (10). Note that a general expression for the inverse of a Hadamard product does not exist. In conclusion, by using the analytic transform (C2) combined with Eq. (B5) we need to explicitly evaluate the kernel matrix only once each time the emulator is invoked, which allows us to efficiently evaluate an array of different KVPs.

Appendix D: Relationships between Kohn Variational Principles and Kohn anomalies
In this Appendix we inspect the relationship between two arbitrary KVPs associated with (L, u) and (L , u ), respectively, and show that their stationary solutions are identical (up to numerical noise) if the cross matrix c defined in Eq. (B7), and thus the generating matrix (B8), is singular (i.e., det c = det ucd = 0). For instance, this applies to (L, L ) = (T, S) and (K −1 , T −1 ), as well as combinations drawn from the generalized T -matrix KVP, and generalized S-matrix KVP, which reduce to the matrices given in Eq. (5) for the Tmatrix and S-matrix, respectively, when τ = 0. (See also Refs. [61,62] and the overview of the KVPs in Ref. [63].) In Appendix C we have shown that C (L) ≡ C is a constant if det c = 0. Hence, in that case, the kernel matrices for the two KVPs are directly (not just elementwisely) proportional to one another and Eq. (C2) reads Kohn (or Schwartz) anomalies occur at energies E > 0 where det ∆ U (u) = 0 or ij (∆ U (u) ) −1 ij = 0. 10 In either case, no (unique) stationary approximation for the L-matrix can be obtained from the functional (7). Generally, Kohn anomalies reduce the KVP's accuracy over a finite range of the phase space (e.g., the energy) since the stationary solution becomes unstable as det ∆ U (u) or ij (∆ U (u) ) −1 ij becomes vanishingly small (in absolute value). This can be seen in Fig. 1.
Our findings imply that the KVPs for (L, L ) with singular cross matrices c [and, for real potentials, also (L, L ) = (S, S −1 )] are equally subject to these Kohn anomalies and that deviations from Eq. (D5) are due to numerical noise, e.g., from inverting ill-conditioned kernel matrices. Apart from these cases, however, we find for real potentials that complex KVPs (e.g., the S-matrix KVP) typically are less prone to Kohn anomalies than real KVPs (e.g., the K-matrix KVP) because their kernel matrices are complex, which means that both the real and imaginary part of det ∆ U (u) or ij (∆ U (u) ) −1 ij need to approach zero simultaneously. (See also the discussion in Ref. [33].) For optical potentials, on the other hand, the kernel matrices are complex whether a real or a complex KVP is used.

Appendix E: Diagnostic tools for Kohn anomalies
As illustrated in Fig. 1, Kohn anomalies can readily be spotted when the emulated scattering observable of interest is plotted as a function of a continuous variable such as the energy and scattering angle-or likewise when the exact scattering solution is known. Since this is not the case in practice and Monte Carlo sampling can be performed at a fixed energy or scattering angle, a different (i.e., automated) method to detect and remove Kohn anomalies is required for such applications.
Our method works as follows. For a given ( , E, N b ), we simultaneously evaluate the complementary KVPs for L = (K, T, T −1 ) and L τ = (L 30 • , L 60 • , L 90 • ) associated with asymptotic boundary conditions drawn from The construction and usage of the parameter matrix (E1) is motivated by the generalized T -and S-matrix KVP (see Appendix D for details), which are redundant for EC trial wave functions. We then compute the relative residuals of, e.g., the estimated S-matrices, defined as for all considered KVPs without repetitions to avoid trivial cases where L 1 = L 2 . 11 For instance, given the KVP estimates for L = (K, T, T −1 ) we would determine δ(K, T ), δ(K, T −1 ), and δ(T, T −1 ). The expressions for the transformations S(L) are discussed in Appendix B. Let P be the set of pairs (L 1 , L 2 ) that fulfill the relative consistency check δ(L 1 , L 2 ) < ε rel , with ε rel = 10 −1 . Since such a consistency check alone does not allow one to disentangle whether L 1 , L 2 , or both are anomalous, we estimate the S-matrix by the weighted sum of averages if at least one consistency check passes (see also Ref. [65]). A small regularization parameter is added to δ(L 1 , L 2 ) to prevent a potential division by zero. If all checks fail, we partition the training set with N b wave functions in batches of size N p N b , 12 remove one batch at a time from the training set (i.e., then of size N b − N p ), and repeat the process iteratively. This usually shifts the Kohn anomalies in each iteration. The iterative process is computationally efficient because ∆ U (u) only needs to be sliced, not recomputed. Different ways of partitioning the training set can be straightforwardly implemented. For instance, suppose N b = 6 and N p = 3, we would first remove the basis wave functions with indices i = (1, 2, 3) and run the consistency checks; if all checks fail again, we add the basis wave functions with i = (1, 2, 3) back to the training sets and remove the ones with i = (4, 5, 6) next. If all checks fail once more, our algorithm signals that Kohn anomalies could not be mitigated for the given training set and terminates. This can happen in practice, e.g., if all KVPs considered are anomalous at overlapping regions in the phase space.