To construct an emulator for bound state calculations, we use here the Rayleigh–Ritz method ² and thus consider the energy functional E ψ ̃ = ⟨ ψ ̃ | H θ | ψ ̃ ⟩ − E ̃ θ ⟨ ψ ̃ | ψ ̃ ⟩ − 1 , (2) where the Lagrange multiplier E ̃ ( θ ) (also known as Ritz value) imposes the normalization condition ⟨ ψ ̃ | ψ ̃ ⟩ = 1 for bound states. The Generalized Ritz Theorem [47] ³ states that the functional (Eq. 2) is stationary about all (discrete) solutions of the Eq. 1, not just the ground state solution, which can be seen by imposing the stationary condition δ E ψ ̃ ≡ 0 = 2 ⟨ δ ψ ̃ | H θ − E ̃ θ | ψ ̃ ⟩ − δ E ̃ θ ⟨ ψ ̃ | ψ ̃ ⟩ − 1 , (3) and noting that Eq. 3 is only fulfilled for arbitrary variations 〈 δ ψ ̃ | if | ψ ̃ 〉 is a solution of the Schrödinger Eq. 1 with E ̃ ( θ ) = E ( θ ) .

Let us now define the trial wave function we use in conjunction with the functional (2): | ψ ̃ 〉 = ∑ i = 1 n b β i | ψ i 〉 ≡ X β ⃗ , (4a) X = | ψ 1 〉 | ψ 2 〉 ⋯ | ψ n b 〉 , (4b) where the column-vector β ⃗ contains the to-be-determined coefficients and the row-vector ⁴ X the (in principle) arbitrary basis states. Here, we use snapshots of high-fidelity solutions of the Eq. 1 at a set of given parameter values; i.e., { | ψ i 〉 ≡ | ψ ( θ i ) 〉 } i = 1 n b [2, 48–50]. No assumption has been made as to how to obtain the high-fidelity solutions.

Figure 1 motivates the efficacy of snapshot-based trial functions. Although a given eigenvector |ψ( θ )⟩ obtained from a high-fidelity solver resides in a high-dimensional (or even infinite-dimensional) space, the trajectory traced out by continuous variations in θ remains in a relatively low-dimensional subspace (as illustrated by the gray plane). Hence, linear combinations of high-fidelity eigenvectors spanning this subspace (i.e., the snapshots) make extremely effective trial wave functions for variational calculations. In nuclear physics, snapshot-based emulators already have accurately approximated ground-state properties, such as binding energies, charge radii [7, 9, 25], and transition matrix elements [9, 29], and have been explored for applications to excited states [51].

Given the trial wave function (4), we determine the coefficients β ⃗ ⋆ that render E [ ψ ̃ = X β ⃗ ] stationary under variations | δ ψ ̃ 〉 = X | δ β ⃗ 〉 of the trial wave function, as opposed to arbitrary variations. Solving for the optimal β ⃗ ⋆ occurs then in the low-dimensional space spanned by the basis elements in X (i.e., the red arrows in Figure 1) rather than in the high-dimensional space in which |ψ⟩ resides. From the stationarity condition Eq. 3, we obtain the reduced-order model [52], H ̃ θ β ⃗ ⋆ θ = E ̃ θ N ̃ β ⃗ ⋆ θ , (5a) β ⃗ ⋆ † θ N ̃ β ⃗ ⋆ θ = 1 , (5b) where H ̃ ( θ ) ≡ X † H ( θ ) X is the subspace-projected Hamiltonian and N ̃ ≡ X † X the norm matrix in the snapshot basis. As opposed to H( θ ) in Eq. 1, H ̃ ( θ ) (and likewise N ̃ ) is a n _b × n _b Hermitian matrix, H ̃ θ = ⟨ ψ 1 | H θ | ψ 1 ⟩ ⋯ ⟨ ψ 1 | H θ | ψ n b ⟩ ⋮ ⋱ ⋮ ⟨ ψ n b | H θ | ψ 1 ⟩ ⋯ ⟨ ψ n b | H θ | ψ n b ⟩ . (6)

In practice, the generalized eigenvalue problem Eq. 5 may experience numerical instabilities due to small singular values in N ̃ . (The instabilities also appear in reduced-order modeling of differential equations, see Section 3). One way to ameliorate these instabilities is to orthonormalize the snapshots in X, hence yielding N ̃ = 1 . This approach could also permit some efficiency gains, by excluding the least important vectors as measured by their singular values. Alternatively, Ref. [53] recently introduced a trimmed sampling algorithm that can substantially reduce the effects of noise in solving generalized eigenvalue problems. Finally, a well-known approach to regularize both generalized eigenvalue problems and matrix inversion is through the use of a nugget [54, 55]. Here, a regularization parameter ν ≪ 1 (called a nugget) is added to the diagonal of the ill-conditioned matrix one wishes to invert (here, N ̃ ), thus shifting its singular values.

By solving the generalized eigenvalue problem Eq. 5, one obtains n _b pairs { E ̃ ( θ ) , β ⃗ ⋆ ( θ ) } consisting of a Lagrange multiplier (i.e., an eigenvalue) and its corresponding coefficient vector. Let us index the eigenvalues of the Eq. 5 and Eq. 1 in ascending order; that is, E ̃ n ⩽ E ̃ n + 1 and E _n ⩽E _n+1, respectively, with n = 1 indicating the lowest eigenvalue. If the snapshot basis X in the trial wave function (Eq. 4) contains n _b linearly independent states, then the Min–Max Theorem [47] asserts that each Lagrange multiplier, E ̃ n θ ⩾ E n θ for 1 ⩽ n ⩽ n b , (7) provides an upper bound on its corresponding eigenvalue of the Schrödinger Eq. 1. ⁵ Furthermore, the Generalized Ritz Theorem implies that the E ̃ n ( θ ) provide not only the variational bounds (Eq. 7) but also stationary approximations for these high-fidelity eigenvalues. Adding another basis state to X can only improve these approximations, which converge to the high-fidelity eigenvalues as the projected subspace approaches the high-fidelity space [47].

Although excited states can also be emulated, especially when adding excited-state snapshots to the trial wave function to improve the emulator’s accuracy (see also Ref. [51]), we focus on ground-state properties and thus use only ground-state snapshots in the trial wave function. For brevity, we will omit the subscripts henceforth. To obtain the approximate ground-state wave function associated with E ̃ ( θ ) , one evaluates the Ritz vector | ψ ( θ ) 〉  ≈ X β ⃗ ⋆ ( θ ) . Expectations values of operators O can then be straightforwardly computed using ⟨ ψ θ | O θ | ψ θ ⟩ ≈ β ⃗ ⋆ † θ O ̃ β ⃗ ⋆ θ , (8) with the subspace-projected O ̃ ( θ ) = X † O ( θ ) X . However, these expectation values generally do not provide variational bounds unless O = H is the Hamiltonian, as discussed (see, e.g., Figure 5 in Ref. [25] for emulated ⁴He ground-state radii).

The reduced-order model (5) can be alternatively derived via a Galerkin projection, as we will also see with the variational emulators for scattering in Section 4. To this end, we construct the weak form of the Eq. 1 by left-multiplying it by an arbitrary test function ⟨ζ| and asserting that ⟨ ζ | H θ − E θ | ψ ⟩ = 0 , ∀ 〈 ζ | . (9) If the weak form (9) is satisfied for all ⟨ζ| for a given set {E, |ψ⟩}, then the set must also satisfy the Eq. 1. The proof of this statement can be obtained via a contrapositive: if Eq. 1 were not satisfied, then one could find a ⟨ζ| such that Eq. 9 is nonzero.

The weak form of the high-fidelity system is the starting point for deriving a reduced-order model. Although Eq. 9 still operates in the large space in which |ψ⟩ resides (cf. Figure 1), we can reduce its dimension by replacing | ψ 〉 → | ψ ̃ 〉 , where | ψ ̃ 〉 is defined in Eq. 4. With the degrees of freedom for |ψ⟩ reduced, we enforce a less strict orthogonality condition: we select n _b test functions ζ _i and assert that the residual due to the trial wave function (cf. Figure 1) should be orthogonal to the subspace Z spanned by these test functions Z = [ | ζ 1 〉 , … , | ζ n b 〉 ] : H θ − E ̃ θ | ψ ̃ 〉 ⊥ Z , (10) or likewise ⟨ ζ | H θ − E ̃ θ | ψ ̃ ⟩ = 0 , ∀ | ζ 〉 ∈ Z . (11) But replacing | ψ 〉 → | ψ ̃ 〉 also implies that the true eigenvalue E is in general not exactly reproduced unless X contains |ψ( θ )⟩. Hence, we also had to apply the approximation E ̃ ≈ E in Eqs. 10 and 11.

In the Galerkin method, which is also known as the “method of weighted residuals,” the test and trial function bases are chosen to be equivalent; i.e., Z = X . The so-called Galerkin condition Eq. 11 is then equivalent to imposing that ⟨ ψ i | H − E ̃ | ψ ̃ ⟩ = 0 holds for i ∈ [1, n _b ]. This yields a system of n _b equations with n _b unknowns β ⃗ and, together with the normalization condition, reduces to Eq. 5 obtained from the variational principle in Section 2.1. However, we stress that the test and trial function bases can be chosen differently (i.e., |ζ _i ⟩ ≠ |ψ _i ⟩), which makes the Galerkin method more general than the variational approach. Note that the normalization condition does not affect the Galerkin condition Eq. 11 and can be implemented by normalizing the trial function.

Figure 2 illustrates the workflow for implementing fast and accurate emulators as described in Section 2.1. The workflow involves

1. a computational framework capable of reliably solving the high-fidelity system Eq. 1,

Several computational frameworks exist in nuclear physics (and quantum chemistry) for solving the few- and many-body Schrödinger Eq. 1 [56]. For illustration, Figure 2 assumes that the high-fidelity solver performs a direct diagonalization of the N _h × N _h Hamiltonian in a chosen (truncated) model basis of length N _h . The corresponding runtime t _s per sampling point θ _i is indicated by the width of the blue bar in Figure 2. In nuclear physics, such approaches are referred to as Configuration Interaction (CI). However, the following discussion will be independent of how the high-fidelity solutions of the Eq. 1 are obtained in practice.

Using the high-fidelity solver, one constructs a set of snapshots { | ψ ( θ i ) 〉 } i = 1 n b in the truncated model basis to build the columns of the N _h × n _b matrix X. The runtime for this task is n _b × t _s . For simplicity, Figure 2 assumes n _b = 3 and depicts the basis functions schematically in different colors. This phase of the emulator needs to be completed only once before the emulator is invoked and is thus called the offline stage as opposed to the online stage of the emulator. The predictions are made quickly and with little memory footprint in the online stage.

The appearance of the full-order Hamiltonian during the offline stage, where the projected Hamiltonian H ̃ ( θ ) ≡ X † H ( θ ) X is computed (see Figure 2), implies that this class of projection-based emulators is intrusive in nature. In general, intrusive emulators apply the basis expansions and projections to the operators implemented in the high-fidelity model [57]. On the other hand, non-intrusive emulators use only outputs of the full-order solver without access to the full-order operators such as the Hamiltonian. Non-intrusive emulators include Gaussian processes [58], Dynamic Mode Decompositions [59, 60], and other machine learning methods [61–63]. More details on this classification scheme can be found in Section 8 in Ghattas and Willcox [57].

The emulator’s efficiency greatly benefits from moving all size-N _h operations into the offline stage, which can easily be achieved for Hamiltonians H( θ ) with an affine parameter dependence. These affine operators can be written as a sum of products of parameter-dependent functions h n ( θ ) and parameter-independent operators H n , H ( θ ) = ∑ n h n ( θ ) H n . (12) Note that the functions h _n ( θ ) are only required to be smooth but not necessarily linear in θ . The affine parameter dependence in Eq. 12 then allows one to store the subspace-projected operators H ̃ n = X † H n X separately up front in the offline phase, from which H ̃ θ = ∑ n h n θ H ̃ n , (13) can be efficiently constructed for each θ _i during the online stage to solve the emulator equation 5. For instance, Hamiltonians derived from chiral EFT can be cast into the form (12) due to their affine dependence on the low-energy couplings. The runtime per sample θ _i in the online phase is therefore typically just a small fraction of that of the high-fidelity solver, as depicted by the small blue box in Figure 2. Likewise, emulating expectation values of other operators with an affine parameter dependence via Eq. 8 also benefits from this offline-online decomposition. For non-affine operators, various hyperreduction methods have been developed to construct approximate affine representations [50, 64], including the empirical interpolation (EIM) [65–68] and gappy proper orthogonal decomposition [69, 70]. See also Refs [48, 71–73] for hyperreduction methods that interpolate X or H ̃ ( θ ) directly, and Refs [33, 74] for recent applications of machine learning tools for hyperreduction.

How should one choose the snapshots in the trial wave function Eq. 4 effectively? For relatively small parameter spaces, one can use Latin hypercube sampling to obtain space-filling snapshots or choose the snapshots in the proximity of the to-be-emulated parameter ranges, keeping n _b ≪ N _h in practice. A chosen set of snapshots expressed in the (truncated) model basis can be optimized by applying a singular value decomposition (SVD) or the closely related proper orthogonal decomposition (POD) [75] to the N _h × n _b matrix X. One then creates a new set of snapshots from the (orthonormal) left-singular vectors associated with the singular values greater than a chosen threshold [64]. This (optional) preprocessing step can be performed during the offline stage, as illustrated in Figure 2, thereby rendering the Eq. 5 an eigenvalue problem (i.e., N ̃ = 1 ) and less sensitive to numerical noise.

The basis states of the trial wave function can also be obtained iteratively, using a greedy algorithm [64, 76, 77]. These algorithms estimate and then minimize the emulator’s overall error by adding basis states (obtained from a high-fidelity solver) in the parameter space where the error is expected to be the largest. Greedy algorithms require fast approximations of the emulator’s error and terminate when either a requested error tolerance or a maximum number of iterations has been achieved. Uncertainty quantification for reduced-order models has been studied in various contexts, including differential equations [64, 77, 78] and nuclear physics problems [24, 43].

The formal results so far in this Section can be illuminated by a simple example, which allows us to compare results from a snapshot-based emulator to more conventional approaches, such as direct diagonalization in a harmonic oscillator basis and Gaussian process emulation. Let us define the system we would like to solve as a single particle with zero angular momentum in three dimensions and trapped in an anharmonic oscillator potential. This example can be directly generalized to few- and many-body systems. The potential operator is the sum of a conventional harmonic oscillator (HO) potential and a finite-range piece: V r ; θ = V HO r + ∑ n = 1 3 θ n ⁡ exp − r 2 / σ n 2 , (14) with σ _n = [0.5, 2, 4] fm. The potential Eq. 14 has the affine structure defined in Eq. 12 for θ and hence can be emulated rapidly after projecting into the snapshot basis during the offline stage. Even the high-fidelity system considered here is still small enough to be solved quickly and accurately using a fine radial mesh on a standard laptop. However, this provides an illuminating setting within which we can observe many qualities seen in more complicated scenarios.

Following the MOR paradigm, we take snapshots of the high-fidelity wave function at various training parameters { θ _i } and collect them into our basis X. Here, we choose n _b = 6 training points randomly and uniformly distributed in the range [−5, 5] MeV for all θ _n ; 50 validation parameter sets are chosen within the same range. The snapshots and the corresponding potentials are shown in Figure 3. These snapshots are then used to construct the reduced-order system as in Eq. 5. All of this, and more, is made simple by the EigenEmulator Python class provided in the supplemental material [20].

Once the reduced system has been constructed and the affine structure of the Hamiltonian exploited to store the projected matrices during the offline stage, we can begin rapid emulation during the online stage. To help provide a baseline to a common approach in nuclear physics, we provide an emulator constructed with the first n _b = 6 HO wave functions as the trial basis X in Eq. 5. We label this approach the HO emulator and the snapshot-based approach the reduced-basis method (RBM) emulator. See Ref. [18] for a guide to the extensive literature on RBMs. One can emulate quantities with this HO approach via our OscillatorEmulator class [20].

For example, we take three of the validation parameter sets we sampled and compare the exact and emulated wave functions for both emulators. Figure 4 shows the results. The gray lines depict the n _b wave functions used to create the reduced-order models, and the colored lines show the emulated results on top of the high-fidelity solutions (black lines). Although both the reduced basis and HO basis are rich enough to capture the main effects of varying θ , the RBM emulator is much more effective at capturing the fine details of the wave function. This can be seen in more detail in Figure 5, where the absolute residuals of the RBM emulator are orders of magnitude smaller than those of the HO emulator. The sensitivity of the emulator accuracy as n _b is varied can be readily studied using the Python code provided on the companion website [20].

The quality of the emulators can be understood by noting in Figure 4 that the basis functions of the RBM emulator match much more closely with the emulated wave functions than the HO emulator, whose wave functions have nodes not seen in the ground state (see the gray lines). Thus, although the HO basis functions may be better at spanning the space of all possible wave functions, they are, in fact, a poor basis for spanning the set of all possible ground states as θ are varied. The RBM emulator constructs an extremely effective basis almost automatically, with minimal input required by the modeler. This can prove particularly effective for cases where the system’s complexity limits the quality of the basis that can be constructed from intuition or expertise alone.

Next, we discuss the emulation of bound-state observables. Straightforward to emulate are the eigen-energies E( θ ), whose emulated values E ̃ ( θ ) are the result of solving the Eq. 5. But as discussed in Section 2.1, other observables associated with the operator O can be emulated via ⟨ ψ ( θ ) | O | ψ ( θ ) ⟩ ≈ ⟨ ψ ̃ ( θ ) | O | ψ ̃ ( θ ) ⟩ using the ψ ̃ ( θ ) found from Eq. 5. We choose to show the results of emulating the radius-squared operator R ², defined here to be ⟨ ψ | R 2 | ψ ⟩ = 1 N ∫ 0 ∞ r 2 d r r 2 ψ 2 r , (15) with the normalization N = ∫ 0 ∞ r 2 d r ψ 2 ( r ) . As stated previously, the bulk of the numerical effort in the evaluation of this matrix element is handled during the offline stage, where the integration is performed once, R ̃ i j 2 ≡ ⟨ ψ i | R 2 | ψ j ⟩ = 1 N ∫ 0 ∞ r 2 d r r 2 ψ i r ψ j r , (16) and then the online stage emulation can occur quickly via Eq. 8; i.e., ⟨ ψ ( θ ) | R 2 | ψ ( θ ) ⟩ ≈ ∑ i j β ⃗ ⋆ ( i ) ( θ ) R ̃ i j 2 β ⃗ ⋆ ( j ) ( θ ) .

For illustrative purposes, we continue our example using the trained RBM and HO emulators, but add a popular emulation tool to the discussion: Gaussian processes (GPs). GPs are non-parametric, non-intrusive machine learning models for both regression and classification tasks [58, 79, 80]. Their popularity stems partly from their convenient analytical form and flexibility in effectively modeling various types of functions. GPs benefit from treating the underlying set of codes as a black box [57]; as we will soon see, this is a double-edged sword. We employ two independent GPs to emulate the ground-state energy and the corresponding radius expectation value. Each GP uses a Gaussian covariance kernel and is fit to the observable values at the same values of θ _i used to train the RBM emulator. We use the maximum likelihood values for the hyperparameters.

The absolute residuals at the validation points for each of the RBM, HO, and GP emulators are shown in Figures 6 and 7 for the energy and radius, respectively. Among these emulators, the GP emulators perform the worst, despite being trained on the values of the energies and radii themselves to perform this very emulation task. Furthermore, its ability to extrapolate beyond the support of its training data is often poor unless great care is taken in the design of its kernel and mean function (see Figures 1 and 2 in Ref. [25]). The GP suffers from what, in other contexts, could be considered its strength: because it treats the high-fidelity system as a black box, it cannot use the structure of the high-fidelity system to its advantage (although some information can be conveyed via physics-informed priors for the hyperparameters). Note that the point here is not that it is impossible to find some GP that can be competitive with other RBM emulators after using expert judgment and careful (i.e., physics-informed) hyperparameter tuning. Rather, we emphasize that with the reduced-order models, remarkably high accuracy is achieved without the need for such expertise.

The HO emulator performs better than the GP emulator, but it was not “trained” per se, it was merely given a basis of the lowest six HO wave functions as a trial basis, from which a reduced-order model was derived. However, the HO emulator can still outperform the GP emulator because it takes advantage of the structure of the high-fidelity system: it is aware that the problem to be solved is an eigenvalue problem, for this is built into the emulator itself. This feature permits a single HO emulator to emulate the wave function, energy, and radius simultaneously.

Coming in first in the comparison of the emulators’ performances is the RBM emulator, which typically results in higher accuracies than the HO and GP emulators by multiple orders of magnitude. The RBM emulator combines the best ideas from the other emulators. Like the GP, the RBM emulator uses evaluations of the eigenvalue problem as training data. However, its “training data” are curves (i.e., the wave functions) rather than scalars (e.g., eigen-energies), like the GP is trained upon. Like the HO emulator, the RBM emulator takes advantage of the structure of the system when projecting the high-fidelity system to create the reduced-order model. With these strengths, the RBM emulator is highly effective in emulating bound-state systems, even with only a few snapshots and far from the support of the snapshots (see Figures 1 and 2 in Ref. [25]). As we will see in Section 4, many of these strengths carry over to systems of differential equations.

Variational principles (VPs) have a long history in physics and play a central role in a wide range of applications; e.g., for deriving equations of motion using stationary-action principles and Euler–Lagrange equations in classical mechanics (see, e.g., Ref. [81] for a historical overview). Here, we use VPs as an alternate way of solving the differential Eq. 17.

Variational principles are based on scalar functionals of the form S [ ξ ] = ∫ Ω d Ω F ξ + ∫ Γ d Γ G ξ , (18) where F and G are differential operators. Many differential Eq. 17 can be solved by finding stationary solutions of a corresponding functional Eq. 18; i.e., the solution ξ _⋆ that leads to δ S [ ξ ⋆ ] = 0 .

However, VPs can also lead straightforwardly to a reduced-order model. This follows from the following trial ansatz | ξ ̃ 〉 ≡ ∑ i = 1 n b β i | ξ i 〉 = X β ⃗ , (19a) X ≡ | ξ 1 〉 | ξ 2 〉 ⋯ | ξ n b 〉 , (19b) with the to-be-determined coefficients vector β ⃗ . Rather than stipulate that δ S = 0 for any arbitrary variation δξ, we instead extract the optimal coefficients, β ⃗ ⋆ , as those for which S is stationary under variations in β ⃗ : ⁶ δ S = ∑ i = 1 n b ∂ S ∂ β i δ β i = 0 . (20) The general case would involve a numerical search for the solution to Eq. 20 but if S is quadratic in ξ, as are all the examples considered here, then the solution can be determined exactly. In this case, S can be written as S [ β ⃗ ] = 1 2 β ⃗ ⊺ A β ⃗ + b ⃗ ⋅ β ⃗ + c (21) for some matrix A, vector b ⃗ , and scalar c. Symmetrizing the quadratic portion—if it is not already symmetric—by rewriting A ← (A + A ^⊺)/2 can be desirable for numerical purposes. It then follows that the optimal coefficients, β ⃗ ⋆ are those for which δ S = A β ⃗ ⋆ + b ⃗ = 0 , (22) which can be solved for β ⃗ ⋆ using standard linear algebra methods. Solving for β ⃗ ⋆ occurs only in a space of size n _b , the number of basis elements { ξ i } i = 1 n b , rather than in the much larger space of ξ itself. Therefore, as long as { ξ i } i = 1 n b approximately spans the space in which ξ lives, the trial function constructed by Eqs. 19a and 22 will be both a fast and accurate emulator of ξ.

Similar to the discussion in Section 2.1, the matrix A in Eq. 22 may be ill-conditioned and require regularization. A nugget ν ≪ 1 can be added to the diagonal elements of A to help stabilize the solution for β ⃗ ⋆ [54, 55].

The Galerkin approach, also more broadly called the “method of weighted residuals,” relies on the weak formulation of the differential Eq. 17. To obtain the weak form, the differential equation and boundary condition (in Eq. 17) are left-multiplied by arbitrary test functions ζ and ζ ̄ and integrated over the domain and boundary, respectively, and then their sum is set to zero: ∫ Ω d Ω ζ D ξ + ∫ Γ d Γ ζ ̄ B ξ = 0 . (23) If Eq. 23 holds for all ζ and ζ ̄ , then Eq. 17) must be satisfied as well. The form of Eq. 23 is often rewritten using integration by parts to reduce the order of derivatives and simplify the solution. Importantly, the weak form has the integral form needed for our emulator application. The weak form and its Galerkin projection are used extensively, e.g., in the finite element method; see Refs. [82–84] for an in-depth study and examples. For a discussion of the convergence properties of the Galerkin method, its relation to abstract variational problems, and other salient mathematical details, see Refs. [64, 85–87]. Here, we follow the introduction of Galerkin methods as provided in Ref. [82].

Starting with the weak form, we can construct an emulator that avoids the need for an explicit variational principle. It begins by first noting that substituting our trial function Eq. 4 into D(ξ) and B(ξ) will not, in general, satisfy Eq. 17 regardless of the choice of β ⃗ . Therefore, there will be some residual, and the goal is to find the β ⃗ ⋆ which minimizes that residual across a range of test functions ζ and ζ ̄ . This system would be over-determined in the case of truly arbitrary test functions, so instead, we propose the test bases | ζ 〉 = ∑ i = 1 n b δ β i | ζ i 〉 , | ζ ̄ 〉 = ∑ i = 1 n b δ β i | ζ ̄ i 〉 , (24) where δβ _i are arbitrary parameters, not related to the β _i in Eq. 19a. The δβ _i will play the same role as those in Eq. 20, namely as a bookkeeping method for determining the set of equations that are equivalently zero. By enforcing that the residuals against these test functions vanish for arbitrary δβ _i , the bracketed expression in δ β i ∫ Ω d Ω ζ i D X β ⃗ ⋆ + ∫ Γ d Γ ζ ̄ i B X β ⃗ ⋆ = 0 , (25) is zero for all i ∈ [1, n _b ], from which the optimal β ⃗ ⋆ are extracted. Because this approximately satisfies the weak formulation, we have found an approximate solution to Equation (17).

In a variety of cases [82], the subspace Z spanned by the test function basis is chosen to coincide with the subspace X spanned by the trial function basis X; i.e., Z = X . This particular choice is known as the Galerkin method, but it is sometimes further specified as the Ritz–Galerkin or Bubnov–Galerkin methods. The Galerkin method is more general than the variational methods described in Section 3.1 because the test space need not be equivalent to the trial space (i.e., Z ≠ X ). In these cases, the approach is described as the Petrov–Galerkin method [82]; this can result in more efficient emulators for some differential equations (84).

The Kohn variational principle (KVP) [91, 92] is one of the most well-known VPs for quantum scattering systems. Here we focus on the KVP flavor that relates a trial wave function to the reactance matrix K. However, alternative flavors exist for other matrices such at T ^± and S (see Section 4.6). Analogously to the Ritz VP for bound states, the KVP then allows us to guess effective wave functions by finding those that make the KVP stationary. This Section will discuss a style of KVP emulator that relies on the homogeneous Schrödinger equation, which requires a normalization constraint during emulation; an alternative style without such a constraint will be discussed in Section 4.2.

We start with a rescaled version of the KVP discussed in Ref. [30]: K ψ ̃ = K ̃ E + ⟨ ψ ̃ | H − E | ψ ̃ ⟩ , (26) where | ψ ̃ 〉 is the trial wave function (denoted by | ξ ̃ 〉 in Section 3) and K ̃ E ≡ ∑ i β i K E , i the associated on-shell trial K matrix with on-shell energy E = q ²/2μ. This flavor of KVP applies when ψ satisfies the asymptotic normalization condition in position space ψ ℓ r → r → ∞ j ℓ q r + n ℓ q r tan δ ℓ , (27) where ϕ(r) = j _ℓ (qr) is the (regular) free-space wave function, and j _ℓ (qr) and n _ℓ (qr) are spherical Bessel and Neumann functions, respectively. ⁷ Note that we define the on-shell K _E matrix as K E = − tan δ ℓ 2 μ q , (28) which differs from the convention in Ref. [30]. The KVP is stationary about exact solutions ψ, such that K [ ψ + δ ψ ] = K E + O ( δ K ) 2 .

Eq. 26 can be cast into the form of the generic functional (Eq. 18) by noting that, in position space, ∫ Γ d Γ G ψ ̃ → G ψ ̃ r = 0 ∞ ≡ W r ϕ , r ψ ̃ ; r 2 μ r = 0 ∞ = K ̃ E , (29) which has defined the surface functional G in Eq. 18 and where we have used the Wronskian W ϕ , ψ ; r ≡ ϕ r ψ ′ r − ϕ ′ r ψ r . (30) Both rϕ(r) and r ψ ̃ ( r ) vanish at r = 0 so only the limit of r → ∞ contributes, from which we can use Eq. 27 when evaluating Eq. 29.

Because the Schrödinger equation is a linear, homogeneous differential equation, the normalization of rψ(r) is proportional to its derivative at, say, r = 0. Therefore, a constraint on the normalization of ψ is equivalent to a boundary condition on ψ′. However, to satisfy this boundary condition we must include a constraint on Eq. 26 if we are to ensure that the trial function ψ ̃ continues to satisfy the normalization condition of Eq. 27. If we assume that each snapshot ψ _i satisfies Eq. 27, then ψ ̃ ℓ r = ∑ i β i j ℓ q r + n ℓ q r ∑ i β i ⁡ tan δ ℓ , i , (31) whose first term implies that we must impose the constraint ∑ _i β _i = 1.

We are now in a position to find the β ⃗ that make Eq. 26 stationary. If we insert the definition of ψ ̃ and K ̃ into Eq. 26, along with the Lagrange multiplier, we have (with repeated indices indicating summations) K β ⃗ = β i K E , i + 1 2 β i Δ U ̃ i j β j + λ ∑ i β i − 1 , (32) where we define V _i = V( θ _i ) and Δ U ̃ i j ≡ ⟨ ψ i | H − E | ψ j ⟩ + i ↔ j = ⟨ ψ i | V θ − V j | ψ j ⟩ + i ↔ j . (33) In the second line we have used that the |ψ _j ⟩ are eigenstates with the corresponding V _j . If V( θ ) is affine in θ , then Δ U ̃ can be projected once in the emulator’s offline stage, and reconstructed quickly during the online stage.

Now we can follow the steps outlined in Section 3.1 to determine β ⃗ ⋆ . Taking the gradient of Eq. 32 with respect to β ⃗ and setting it equal to 0 yields K ⃗ E + Δ U ̃ β ⃗ ⋆ + λ ⋆ = 0 , (34) where K ⃗ E are the n _b on-shell K-matrices used to train the emulator, and β ⃗ ⋆ are the optimal coefficients of the trial wave function. The gradient with respect to λ simply returns the constraint. This system can be solved via the system of equations Δ U ̃ 1 ⃗ 1 ⃗ ⊺ 0 β ⃗ ⋆ λ ⋆ = − K ⃗ E 1 , (35) where 1 ⃗ is an n _b × 1 vector of ones. If the n _b number of basis functions is much smaller than the size of ψ, then Eq. 35 can be a highly computationally efficient emulator for scattering systems and requires little computer memory to store. The on-shell K matrix can then be emulated via K E ≈ K β ⃗ ⋆ = β ⃗ ⋆ ⋅ K ⃗ E + 1 2 β ⃗ ⋆ ⊺ Δ U ̃ β ⃗ ⋆ , (36) whose operations all occur quickly in the size-n _b space during the online stage.

The derivation above followed closely that of Refs [30, 93], but one could instead arrive at exactly Eq. 35 from a Galerkin projection [43]. Rather than begin with the VP, we start here with (the strong form of) the homogeneous Schrödinger equation, i.e., H( θ )|ψ⟩ = E|ψ⟩. To construct the weak form, we multiply with a test function |ζ⟩ and assert that the residual vanishes: ⟨ ζ | H − E | ψ ⟩ = 0 , ∀ | ζ 〉 . (37) To make explicit the boundary conditions, we make use of the relation 0 = ⟨ ζ | H − E | ψ ⟩ = ⟨ ζ | H ⃖ − E | ψ ⟩ − W r ζ , r ψ ; r 2 μ r = 0 ∞ , (38) where we have again used the Wronskian from Eq. 30, and assigned H ⃖ as the operator acting, after integration by parts, on ⟨ζ| instead of |ψ⟩. By adding Eqs. 38 and 37, we have ⟨ ζ | H − E | ψ ⟩ + ⟨ ζ | H ⃖ − E | ψ ⟩ = W r ζ , r ψ ; r 2 μ r = 0 ∞ . (39) This is the weak form of the homogeneous Schrödinger equation that we will use to construct the emulator, although the asymptotic normalization condition Eq. 27 still needs to be enforced. This will be imposed via a Lagrange multiplier after inserting our trial basis.

Now that we have a weak form, the next step to construct the reduced-order model equations is to define our trial and test bases to project the weak form into the finite space of these bases. To align with the Kohn emulator from the variational argument above, we choose the trial and test basis to be identical as snapshots ψ _i . Then we can evaluate W r ψ i , r ψ j ; r 2 μ r = 0 ∞ = K j − K i (40) and thus, it follows after including the Lagrange multiplier that λ + ⟨ ψ i | H − E | ψ j ⟩ + ⟨ ψ i | H ⃖ − E | ψ j ⟩ β j = ∑ j β j K j − K i ∑ j β j . (41) The sum in the rightmost term can be evaluated using the constraint ∑ _j β _j = 1, and we can make the redefinition λ′ ≡ λ−∑ _j β _j K _j without impacting the solution because this term does not depend on i. Thus, we have λ ′ + K ⃗ E + Δ U ̃ β ⃗ ⋆ = 0 , (42) which is exactly Eq. 34 found by making the KVP stationary. This simplification can be understood by noting that if { β ⃗ ⋆ , λ ⋆ } satisfy Eq. 41, then we know that { β ⃗ ⋆ , λ ⋆ ′ } is the unique solution to Eq. 42. Therefore, we can solve Eq. 42 to obtain β ⃗ ⋆ rather than Eq. 41. In conclusion, using the Galerkin projection of the homogeneous Schrödinger equation with trial and test bases of ψ _i , we were able to obtain the same coefficients as the KVP in Eq. 35, which yield the same on-shell K matrix when used in Eq. 36.

The Kohn emulators from Section 4.1 start with the homogeneous Schrödinger equation, which does not enforce any specific normalization of the wave function; hence this requirement needs to be enforced at the time of emulation. This effectively takes the n _b degrees of freedom {ψ _i }—which were potentially costly to obtain—and reduces the degrees of freedom to n _b −1. However, one can instead build in the normalization from the very start, thus removing the need to constrain our basis via ∑ j = 1 n b β j = 1 during emulation. The unconstrained emulator is fundamentally different from any approach that constrains the coefficients (e.g., explicit substitution of β 1 = 1 − ∑ j = 2 n b β j ), regardless of if a Lagrange multiplier is explicitly used as in Section 4.1. This is the topic of the current section.

The full wave function |ψ⟩ can be written as the sum of the free wave function |ϕ⟩ and the scattered wave |χ⟩, that is, |ψ⟩ = |ϕ⟩ + |χ⟩. Thus, we can rewrite the KVP as K = K + ⟨ ψ | H − E | ψ ⟩ = K + ⟨ χ | H − E | χ ⟩ + ⟨ ϕ | H − E | χ ⟩ + ⟨ ϕ | H − E | ϕ ⟩ + ⟨ χ | H − E | ϕ ⟩ = ⟨ χ | H − E | χ ⟩ + ⟨ ϕ | V | χ ⟩ + ⟨ ϕ | H − E | ϕ ⟩ + ⟨ χ | V | ϕ ⟩ , (43) where we used (via integration by parts) ⟨ ϕ | H − E | χ ⟩ − ⟨ ϕ | H ⃖ − E | χ ⟩ = − K . (44) We choose our trial function as | χ ̃ 〉 , which always enforces the normalization condition | ψ ̃ 〉 = | ϕ 〉 + | χ ̃ 〉 , and so no additional constraint needs to be included in the variational principle.

Now we can construct the set of linear equations that makes Eq. 43 stationary in | χ ̃ 〉 = ∑ i β i | χ i 〉 . By taking the gradient with respect to β _i , we find Ω β ⃗ ⋆ = ω ⃗ , (45) where Ω i j = ⟨ χ i | E − H | χ j ⟩ , (46a) ω i = ⟨ χ i | V | ϕ ⟩ , (46b) which is the set of equations used to obtain β ⃗ ⋆ . The matrix elements Ω _ij can be evaluated with the help of E − H | χ j 〉 = E − H j | χ j 〉 + V j − V | χ j 〉 = V j | ϕ 〉 + V j − V | χ j 〉 , (47) with H _j = H ( θ _j ) and V _j = V ( θ _j ).

An equivalent approach follows from a Galerkin orthogonalization procedure. We begin by writing the homogeneous Schrödinger equation in inhomogeneous form using |ψ⟩ = |ϕ⟩ + |χ⟩: E − H | χ 〉 = V | ϕ 〉 . (48) We can construct the weak form by multiplying by a generic test function |ζ⟩, which yields ⟨ ζ | E − H | χ ⟩ = ⟨ ζ | V | ϕ ⟩ . (49) Next, we insert the trial function | χ ̃ 〉 and choose the test basis of { | χ i 〉 } i , which is the same as the trial basis. This yields a reduced weak form that is identical to Eq. 45.

We have shown that the coefficients β ⃗ ⋆ found via the appropriate Galerkin procedure aligns exactly with the KVP. However, we can go one step further and in fact derive an estimate for the K matrix that is equivalent to K [ β ⃗ ⋆ ] . By inserting the optimal coefficients into K|ϕ⟩ = V|ψ⟩, ⟨ ϕ ′ | K | ϕ ⟩ ≈ ⟨ ϕ ′ | V | ϕ ⟩ + ⟨ ϕ ′ | V | χ ̃ ⟩ = ⟨ ϕ ′ | V | ϕ ⟩ + ∑ i j ⟨ ϕ ′ | V | χ i ⟩ Ω − 1 i j ⟨ χ j | V | ϕ ⟩ , (50) with the factors in brackets equating to β _i using Eq. 45. The equivalence to the KVP is demonstrated in Refs. [94, 95].

The Schwinger variational principle (SVP) is given by [94]. K ψ ̃ = ⟨ ψ ̃ | V | ϕ ⟩ + ⟨ ϕ | V | ψ ̃ ⟩ − ⟨ ψ ̃ | V − V G 0 V | ψ ̃ ⟩ , (51) where G ₀ is the Green’s operator. This too has the stationary property K [ ψ + δ ψ ] = K + O ( δ K ) 2 when ψ is a wave function satisfying the LS equation. Following the MOR philosophy and inserting a trial function ψ ̃ , the stationary condition becomes W β ⃗ ⋆ = w ⃗ , (52) where W i j = ⟨ ψ i | V − V G 0 V | ψ j ⟩ (53a) w i = ⟨ ψ i | V | ϕ ⟩ , (53b) for all i ∈ [1, …, n _b ].

The system of Eq. 52 can also be determined by a Galerkin projection procedure. In this case, we start with the LS equation for wave functions, | ψ 〉 = | ϕ 〉 + G 0 V | ψ 〉 , (54) and create a weak form by left-multiplying by V( θ ) along with the test function |ζ⟩: ⟨ ζ | V | ψ ⟩ = ⟨ ζ | V | ϕ ⟩ + ⟨ ζ | V G 0 V | ψ ⟩ . (55) The weak form can then be converted to its discrete form by setting ψ → ψ ̃ and enforcing orthogonality against |ζ _i ⟩ = |ψ _i ⟩ for i ∈ [1, …, n _b ]. ⁸ This yields then Eq. 52, and so the coefficients found by making Eq. 51 stationary are indeed identical to those found via the Galerkin procedure for Eq. 54.

Using the emulation of ψ, which is calculated by inserting the optimal coefficients obtained from Eq. 52 into the definition of ψ ̃ , we can get the associated K through ⟨ ϕ ′ | K | ϕ ⟩ = ⟨ ϕ ′ | V | ψ ⟩ ≈ ⟨ ϕ ′ | V | ψ ̃ ⟩ = ∑ i j ⟨ ϕ ′ | V | ψ i ⟩ W − 1 i j ⟨ ψ j | V | ϕ ⟩ . (56) This Equation is exactly the solution for K found via the LS equation while assuming a finite-rank approximation for V: V f = ∑ i j V | ψ i 〉 Λ i j 〈 ψ j | V , (57) where Λ − 1 i j = ⟨ ψ i | V | ψ j ⟩ . (58) It is known that the SVP yields a K matrix that is equivalent to that found via a finite-rank approximation to V [94, 95], which shows that the Galerkin projection described in this Section is identical to the SVP.

The Newton variational principle (NVP) for the K matrix is given by [31, 96]. K K ̃ = V + V G 0 K ̃ + K ̃ G 0 V − K ̃ G 0 K ̃ + K ̃ G 0 V G 0 K ̃ , (59) where K ̃ is a trial matrix. If desired, one could instead emulate T ^(±) by imposing the associated boundary conditions on G ₀. Here it is assumed that we have chosen an on-shell energy E, which will remain implicit throughout. A separate emulator can be constructed for each choice of E. The functional Eq. 59 is stationary about exact solutions of the LS equation, i.e., K [ K + δ K ] = K + ( δ K ) 2 . If we write the trial matrix as a linear combination of exact snapshots K ̃ = ∑ i = 1 n b β i K i , (60) then we can construct an emulator of the K matrix in the spirit of the RBM.

Unlike some of the VPs discussed so far, the NVP is written here in operator form, without yet projecting it into a basis. This gives us the freedom to assert that any component ⟨ ϕ ′ | K | ϕ ⟩ constructed from Eq. 59 is stationary, which yields an emulator for ⟨ϕ′|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. Coupled channels could be emulated by choosing the angular momentum quanta differently between |ϕ′⟩ and |ϕ⟩. In fact, the NVP does not even require |ϕ⟩ and |ϕ′⟩ to be free-space states; one can impose stationarity of Eq. 59 between any two states due to its operator form. Note that the independent coefficients β ⃗ are found for each choice of |ϕ′⟩ and |ϕ⟩. Thus, in the case of coupled partial waves, for example, each channel is emulated independently. To compute phase shifts, we must emulate K at the on-shell energy E = q ²/2μ and thus, k = k′ = q for |ϕ⟩ = |kℓm⟩ and ⟨ϕ′| = ⟨k′ℓ′m′|.

Expressed in the chosen basis, simplifying the functional Eq. 59 after inserting Eq. 60 yields [31] ⟨ ϕ ′ | K θ , β ⃗ | ϕ ⟩ = ⟨ ϕ ′ | V θ | ϕ ⟩ + β ⃗ ⊺ m ⃗ θ − 1 2 β ⃗ ⊺ M θ β ⃗ , (61) with m i θ = ⟨ ϕ ′ | K i G 0 V θ + V θ G 0 K i | ϕ ⟩ , (62a) M i j θ = ⟨ ϕ ′ | K i G 0 K j − K i G 0 V θ G 0 K j | ϕ ⟩ + i ↔ j . (62b) If the potential V( θ ) has an affine parameter dependence, m ⃗ and M can be efficiently constructed by linear combinations of matrices pre-computed during the emulator’s offline stage, resulting in substantial improvements in CPU time, e.g., for chiral interactions.

By imposing the stationary condition ∂ K / ∂ β ⃗ = 0 , one then finds β ⃗ ⋆ ( θ ) such that M β ⃗ ⋆ = m ⃗ . Given that the optimal β ⃗ ⋆ ( θ ) yields a trial matrix Eq. 60 with an error δK, one can insert β ⃗ ⋆ in Eq. 61 to obtain an error (δK)². The resulting emulator K ⋆ ( θ ) ≡ K ( θ , β ⃗ ⋆ ) is then [31] ⟨ ϕ ′ | K | ϕ ⟩ ≈ ⟨ ϕ ′ | K | ϕ ⟩ = ⟨ ϕ ′ | V | ϕ ⟩ + 1 2 m ⃗ ⊺ M − 1 m ⃗ . (63)

Reference [31] studied several applications of the emulator Eq. 63 to short-range potentials with and without the Coulomb interaction and partial-wave coupling. They demonstrated that the NVP emulator has remarkable extrapolation capabilities (see Figure 2 in Ref. [31]) and can quickly reproduce high-fidelity calculations of neutron-proton cross sections based on modern chiral interactions with negligible error.

We now repeat the derivation for the NVP emulator, but instead from the perspective of a Galerkin projection. Here, we will focus on the case where ⟨ϕ′| = ⟨ϕ|. We start with the LS equation K = V + V G 0 K , (64) which, in this context, constitutes the strong form of the integral equation. Although Eq. 64 is written in terms of abstract operators, it can be turned into a vector equation in a specific representation after right-multiplying by |ϕ⟩. To derive the weak form we left-multiply by G ₀ and a test function ⟨ζ|: ⟨ ζ | G 0 K − G 0 V G 0 K | ϕ ⟩ = ⟨ ζ | G 0 V | ϕ ⟩ . (65) The trial function in this case is K|ϕ⟩, which can be expanded in a discrete (snapshot) basis using Eq. 60. We further employ the Galerkin prescription, where the test basis is equivalent to the trial basis, making ⟨ζ _i | = ⟨ϕ|K _i . With these assumptions, the reduced weak form becomes M β ⃗ ⋆ = m ⃗ , (66) with M and m ⃗ defined in Equation (62), again with ⟨ϕ′| = ⟨ϕ|. Thus, we find the same β ⃗ ⋆ using either the NVP or the Galerkin projection.

Given the optimal coefficients β ⃗ ⋆ , the emulator can be derived by substituting K ̃ into the right-hand side of Eq. 64: ⟨ ϕ | K | ϕ ⟩ ≈ ⟨ ϕ | V | ϕ ⟩ + ⟨ ϕ | V G 0 K ̃ | ϕ ⟩ = ⟨ ϕ | V | ϕ ⟩ + 1 2 m ⃗ ⊺ M − 1 m ⃗ , (67) which is equivalent to Eq. 63 under the assumption that ⟨ϕ′| = ⟨ϕ|. Therefore, both the NVP and Galerkin projection lead to identical emulators for the K matrix.

The scattering emulators discussed so far are best known as VPs but are equivalent to various types of Galerkin projections of the Schrödinger or LS equation (see Table 1). However, other types of emulators can be constructed via Galerkin projections, even if they do not necessarily correspond to any well-known VP.

Starting from the Schrödinger equation—a second-order differential equation—we must impose two boundary conditions. The first is that rψ(r) vanishes at r = 0; this constraint has been automatically satisfied by our choice of trial bases in all VPs considered above. But the second constraint is yet to be chosen. In the KVP, for example, the second constraint was obtained via the normalization of rψ(r) as r → ∞, which, in the constrained KVP, led to a normalization constraint for the coefficients β _i . Because the Schrödinger equation is linear and homogeneous, this normalization condition is equivalent to imposing a constraint on the derivative of rψ(r), e.g., evaluated at the origin.

Thus, an alternative weak form for the Schrödinger equation could be constructed using only constraints at the origin. Let us construct a coordinate-space emulator with (rψ)′(0) = 1. Starting from the generic weak form (23), we obtain ⟨ ζ | H − E | ψ ⟩ + ζ ̄ r ψ ′ − 1 r = 0 = 0 , (68) where ζ and ζ ̄ are the (independent) test functions in the domain and on the boundary, respectively, and the boundary condition is only evaluated at the origin. Here, we can make the Galerkin choice of (domain) test functions, where ⟨ζ| = ⟨ψ|, but make a Petrov-Galerkin choice for the boundary, with ζ ̄ ( 0 ) = 1 .

Thus, the discretized weak form, from which our emulator equations follow, is given by ⟨ ψ i | H − E | ψ j ⟩ β j + ∑ j β j − 1 = 0 , (69) where we have assumed that the trial basis is constructed such that each snapshot satisfies ( r ψ j ) ′ ( 0 ) = 1 .

As an example, we show the output of such an s-wave (ℓ = 0) emulator in Figure 8. Here, the potential is given by a sum of two Gaussians, V r , θ = θ 1 ⁡ exp − κ 1 r 2 + θ 2 ⁡ exp − κ 2 r 2 , (70) with κ ₁ = 0.5 and κ ₂ = 1. The parameters to be varied are θ = {θ ₁, θ ₂}. The six training and one validation parameters are selected randomly from a uniform distribution in the range of [−5, 5]. To obtain the snapshots, the partial-wave decomposed radial Schrödinger equation for ℓ = 0 can be expressed as the system of coupled first-order differential equations, y ′ 0 ( r ) y ′ 1 ( r ) = ( r ψ ) ′ ( r ) 2 μ V ( r ; θ i ) − E ( r ψ ) ( r ) , (71) and numerically solved with Runge-Kutta methods. For more details on solving the radial Schrödinger equation and matching the solutions to the asymptotic boundary condition (27), see, e.g., Ref. [97]. As we can see from Figure 8, each training and emulated wave function has matching boundary conditions at the origin, and the discrepancy from the true wave function is less than 10^–3.

The KVP functional given in Eq. 26 can be extended to include arbitrary boundary conditions [32, 98]. For simplicity, let us consider short-range potentials V( θ ) that have been partial-wave decomposed into an uncoupled channel with angular momentum ℓ. The general asymptotic form of the (coordinate-space) radial wave functions will be linear combinations of free-space solutions ⁹ ψ ℓ , E ( r ) → r → ∞ ϕ ̄ ℓ , E ( 0 ) ( r ) + L ℓ , E ϕ ̄ ℓ , E ( 1 ) ( r ) , (72) where ϕ ̄ ℓ , E ( 0 ) r ϕ ̄ ℓ , E ( 1 ) r = N − 1 u 00 u 01 u 10 u 11 j ℓ q r η ℓ q r , (73) with q = 2 μ E and an arbitrary normalization constant N ≠ 0 . Here, L _ℓ,E is a generic scattering matrix that is determined by the boundary condition, as parametrized by the non-singular matrix u .

We now define L as a general functional for the generic L-matrix in Eq. 72 [32, 88, 98], L [ ψ ̃ ] = L ℓ , E + 2 μ N 2 q d e t u 〈 ψ ̃ u | H − E | ψ ̃ u 〉 . (74) Note that Eq. 74 is not restricted to coordinate space, e.g., it also holds for scattering in momentum space (see Ref. [88]). With Eq. 26, one can follow the process described in Section 4.1 to emulate any asympototic boundary condition. Obtaining an emulator prediction for different boundary conditions does not mean that Eq. 74 has to be solve multiple times. In fact, it only needs to be solved once and each term in the functional rescaled using the relations derived in Ref. [32]: Δ U ̃ u ′ = C ′ − 1 L i C ′ − 1 L j det u det u ′ Δ U ̃ u , (75) C ′ L = det u det u ′ u 11 ′ − u 10 ′ K L u 11 − u 10 K L . (76)

The non-primed terms refer to the initial state and primed terms refer to the final state (explained below). The snapshots used to train the emulator in the offline stage are transformed using the Möbius (or linear fractional) transform L ′ L = − u 01 ′ + u 00 ′ K L u 11 ′ − u 10 ′ K L . (77)

Let us consider solving Eq. 74 using the K-matrix boundary condition, but then wanting a prediction for the T-matrix. We would first rescale Δ U ̃ using Eq. 75. Here, u and u ′ would correspond to u _K and u _T , respectively, given by u K = 1 0 0 1 , u T = 1 0 i 1 . (78)

Once Δ U ̃ ( u ′ ) is calculated and the snapshots transformed from the K- to the T-matrix according to Eq. 77, we apply Eq. 35 to obtain the emulator prediction for the T-matrix. One can also inverse transform the new emulated solution back to its K-matrix equivalent by using K L = u 01 + u 11 L u 00 + u 10 L . (79)

Variational principles may not always provide a (unique) stationary approximation, causing the appearance of spurious singularities known as Kohn (or Schwartz) anomalies [32, 98, 101], which can render applications of VPs ineffective; especially for sampling of a model’s parameter space. ¹⁰ The appearance of these anomalies depends on the parameters θ used to train the emulator in the offline stage, the scattering energy, and the evaluation set used in the online stage. However, Ref. [32] demonstrated that a KVP-based emulator that simultaneously emulates an array of KVPs with different boundary conditions can be used to systematically detect and remove these anomalies. The anomalies can be detected by assessing the (relative) consistency of the different emulated results for, e.g., the scattering S-matrix. The results that do not pass the consistency check are discarded and the remaining ones averaged to obtain an anomaly-free estimate of the S matrix (or any other matrix). If all possible consistency checks fail, one can change the basis size of the trial wave function iteratively, which typically shifts the locations of the Kohn anomalies in the parameter space in each iteration. The basic idea for removing Kohn anomalies is general and can be applied to other emulators, including the NVP-based emulator discussed in Section 4.4, as long as multiple scattering boundary conditions can be emulated independently and efficiently. Alternatively, one might also consider comparing the consistency of emulated results obtained from different VPs such as the ones summarized in Table 1.