# BUQEYE guide to projection-based emulators in nuclear physics

^{1}Department of Physics and Astronomy & Institute of Nuclear and Particle Physics, Ohio University, Athens, OH, United States^{2}Facility for Rare Isotope Beams, Michigan State University, East Lansing, MI, United States^{3}Department of Physics, The Ohio State University, Columbus, OH, United States

The BUQEYE collaboration (Bayesian Uncertainty Quantification: Errors in Your effective field theory) presents a pedagogical introduction to projection-based, reduced-order emulators for applications in low-energy nuclear physics. The term *emulator* refers here to a fast surrogate model capable of reliably approximating high-fidelity models. As the general tools employed by these emulators are not yet well-known in the nuclear physics community, we discuss variational and Galerkin projection methods, emphasize the benefits of offline-online decompositions, and explore how these concepts lead to emulators for bound and scattering systems that enable fast and accurate calculations using many different model parameter sets. We also point to future extensions and applications of these emulators for nuclear physics, guided by the mature field of model (order) reduction. All examples discussed here and more are available as interactive, open-source Python code so that practitioners can readily adapt projection-based emulators for their own work.

## 1 Introduction

Nuclear systems are notoriously complex. But typically, our theoretical modeling of nuclear phenomena contains superfluous information for quantities of interest. Model order reduction (MOR) refers to powerful techniques that enable us to reduce a system’s complexity systematically (*e*.*g*., see Refs. [1–3] for comprehensive introductions). These techniques enable emulators, which are low-dimensional surrogate models capable of rapidly and reliably approximating high-fidelity models, making practical otherwise impractical calculations. But the nuclear physics community has barely scratched the surface of the types of emulators that could be crafted or explored their full range of applications.

A fertile area for new emulators is uncertainty quantification (UQ) [4–13] in nuclear physics, which is the general theme of this Frontiers Research Topic [14]. Quantifying theoretical uncertainties rigorously is crucial for comparing theory predictions with experimental and/or observational constraints and performing model comparison and/or mixing [15]. However, UQ has only recently drawn much attention as nuclear theory has entered the precision era. Bayesian parameter estimation for nuclear effective field theory (EFT) and optical models, UQ for nuclear structure pushing toward larger masses and for reactions across the chart of nuclides, experimental design [15–17] for the next-generation of precision experiments probing the nuclear dripline, and many other applications will all benefit from emulators. This Research Topic [14] already contains several new applications of emulators for nuclear physics. Key to the wider adoption of these tools is the evangelization of their potential and the creation of pedagogical guides for those first starting in this field [18]. This article is aimed at both goals.

To do so, the BUQEYE collaboration (Bayesian Uncertainty Quantification: Errors in Your EFT) [19] has created a rather unconventional document comprised of the article you are reading now along with a companion website [20] containing interactive supplemental material and source code that generates all the results shown, and much more. Interested individuals can dynamically generate different versions of this document based on tunable parameters. We hope that this format encourages readers to experiment and build upon the examples presented here, thereby facilitating new applications.

Various types of emulators have already been applied with success within nuclear physics. A non-exhaustive list of applications includes Refs. [7, 9, 21–42]. But as emphasized in Refs. [18, 43], there is a broad and relatively mature MOR literature outside of nuclear physics waiting to be exploited (*e*.*g*., see Ref. [44] for an overview of the Universe of MOR approaches). Our goal in this guide will be to facilitate this exploitation through a selective treatment of physics-informed, projection-based emulators relevant to a wide range of nuclear physics problems.

To this end, we organize this guide as follows. Section 2 focuses on emulators for bound-state calculations using subspace-projection methods. We then provide a more general introduction to MOR for solving differential equations in Section 3, which leads to our discussion of scattering emulators in Section 4. Section 5 concludes with a summary and outlook. Throughout, we draw connections between variational and Galerkin projection methods and illustrate these concepts with pedagogical examples, supplemented by source code on the companion website [20].

## 2 Eigen-emulators

In this section, we discuss the construction of fast and accurate emulators for bound-state calculations. Given a (Hermitian) Hamiltonian *H*(** θ**) parameterized by

*E*(

**), |**

*θ**ψ*(

**)⟩} of the Schrödinger Equation**

*θ*subject to the normalization ⟨*ψ*(** θ**)|

*ψ*(

**)⟩ = 1. The components of the vector**

*θ***may be model parameters, such as the low-energy couplings of a nuclear EFT, or other parameters describing the system of interest [33, 40]. We consider here cases in which Eq. 1 can be solved with high fidelity, but doing so requires a significant amount of compute time. This compute time is compounded when repeated solutions are required throughout the parameter space,**

*θ**e.g.*, during optimization routines or Monte Carlo sampling. In the following, we will discuss how the Ritz variational principle and the Galerkin method can be used to construct rapid and reliable

^{1}emulators that facilitate these calculations.

### 2.1 Variational approach

To construct an emulator for bound state calculations, we use here the Rayleigh–Ritz method^{2} and thus consider the energy functional

where the Lagrange multiplier ^{3} 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

and noting that Eq. 3 is only fulfilled for arbitrary variations

Let us now define the trial wave function we use in conjunction with the functional (2):

where the column-vector ^{4} *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.*,

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].

**FIGURE 1**. Illustration of a projection-based emulator using only two snapshots |*ψ*_{i}⟩ ≡ |*ψ*(*θ*_{i})⟩ (dark gray points). These snapshots are high-fidelity solutions of the Schrödinger Eq. 1, which span the subspace of the reduced-order model, as indicated by the red arrows and the gray plane. The trajectory of a high-fidelity eigenvector is denoted by the blue curve. The orange dot depicts an eigenvector |*ψ*(** θ**)⟩ along the trajectory that, when projected onto the reduced space, corresponds to the turquoise point; hence, the difference between the orange and turquoise points represents the error due to the emulator’s subspace projection (

*i.e.*, the dotted line). Inspired by Figure 2.1 in Ref. [2].

Given the trial wave function (4), we determine the coefficients *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],

where *H*(** θ**) in Eq. 1,

*n*

_{b}×

*n*

_{b}Hermitian matrix,

In practice, the generalized eigenvalue problem Eq. 5 may experience numerical instabilities due to small singular values in *X*, hence yielding *ν* ≪ 1 (called a nugget) is added to the diagonal of the ill-conditioned matrix one wishes to invert (here,

By solving the generalized eigenvalue problem Eq. 5, one obtains *n*_{b} pairs *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}, 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,

provides an upper bound on its corresponding eigenvalue of the Schrödinger Eq. 1.^{5} Furthermore, the Generalized Ritz Theorem implies that the *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 *O* can then be straightforwardly computed using

with the subspace-projected *O* = *H* is the Hamiltonian, as discussed (see, *e*.*g*., Figure 5 in Ref. [25] for emulated ^{4}He ground-state radii).

### 2.2 Galerkin approach

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

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 *ψ*⟩ 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

or likewise

But replacing *E* is in general not exactly reproduced unless *ψ*(** θ**)⟩. Hence, we also had to apply the approximation

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.*, *i* ∈ [1, *n*_{b}]. This yields a system of *n*_{b} equations with *n*_{b} unknowns *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.

### 2.3 Emulator workflow and offline-online decomposition

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,

2. the snapshot-based trial wave function Eq. 4 with the optimal coefficients (*i.e.*, the weights) determined by the Eq. 5, and

3. an efficient offline-online decomposition in which the computational heavy lifting is performed once *before* the emulator is invoked.

**FIGURE 2**. Illustration of the workflow for implementing fast and accurate emulators, including a high-fidelity solver (left) and an intrusive, projection-based emulator with efficient offline-online decomposition (right), for sampling the (approximate) solutions of the Schrödinger Eq. 1 in the parameter space ** θ**. For brevity, the figure assumes that the snapshots are orthonormalized during the offline stage such that

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 *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 *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

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**

*θ*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

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.*,

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].

### 2.4 Illustrative example

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:

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].

**FIGURE 3**. Basis functions for training the snapshot-based eigen-emulator. The black curves show the potential, and the blue curves show the wave functions as functions of the radial coordinate *r*. The wave functions are offset vertically by their corresponding energies for clarity. See the main text for details.

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].

**FIGURE 4**. Emulated wave functions for the RBM emulator (top panel) and HO emulator (bottom panel) as a function of the radius. The solid black lines represent the exact solution, and the dots represent the emulator result. The gray lines give the wave functions used to train the emulator. See the main text for details.

**FIGURE 5**. Absolute residuals of the emulated wave functions (in fm^{−1/2}) based on the RBM and HO emulator as a function of *r*, with colors corresponding to those in Figure 4. The emulator results are compared to the exact solutions. See the main text for details.

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

*O*can be emulated

*via*

*R*

^{2}, defined here to be

with the normalization

and then the online stage emulation can occur quickly *via* Eq. 8; *i.e.*,

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.

**FIGURE 6**. Absolute residuals in the energy (in MeV) at the 50 validation points for the RBM, HO, and GP emulators. The validation points are chosen randomly from a uniform distribution within the same range as the training points. See the main text for details.

**FIGURE 7**. Similar to Figure 6 but for the root-mean-squared radius

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.

## 3 Model reduction

In this Section, we provide a more general discussion of variational principles and the Galerkin method as the foundations for constructing highly efficient emulators for nuclear physics (see also Ref. [18]). The general methods discussed here will be used as a springboard to develop emulators for the specific case of scattering systems in Section 4.

We consider (time-independent) differential equations that depend on the parameter vector ** θ** and aim to find the solution

*ξ*of

where {*D*, *B*} are differential operators and Ω is the domain with boundary Γ. See Ref. [18] for illustrative examples. Here, we use the generic function *ξ* because different choices of *ξ* will be made in Section 4. In what follows, we will discuss two related methods for constructing emulators from Eq. 17, which states the physics problem in a *strong form* (*i.e.*, Eq. 17 holds for each point in the domain and on the boundary). The first begins by finding a variational principle whose stationary solution implies Eq. 17. The second instead constructs the corresponding *weak form* of Eq. 17.

### 3.1 Variational principles

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

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

However, VPs can also lead straightforwardly to a reduced-order model. This follows from the following trial ansatz

with the to-be-determined coefficients vector *δξ*, we instead extract the optimal coefficients, ^{6}

The general case would involve a numerical search for the solution to Eq. 20 but if *ξ*, as are all the examples considered here, then the solution can be determined exactly. In this case,

for some matrix *A*, vector *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,

which can be solved for *n*_{b}, the number of basis elements *ξ* itself. Therefore, as long as *ξ* 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

### 3.2 Galerkin emulators

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

If Eq. 23 holds for all *ζ* and *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 *residual*, and the goal is to find the *ζ* and

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

is zero for all *i* ∈ [1, *n*_{b}], from which the optimal

In a variety of cases [82], the subspace *X*; *i.e.*, *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.*,

## 4 Scattering emulators

In this Section, we describe various reduced-basis emulators one could construct for quantum scattering systems. Throughout, we note how the variational principles used to construct emulators in recent works are related [30–33, 88]. We also describe how each of the results from VPs could instead be derived from Galerkin projections.

For scattering problems, the Equation 1 is no longer an eigenvalue problem. The task is to solve the differential equation for the wave function at a given energy *E* rather than searching for discrete energies with normalizable wave functions. Differential equations are well studied in the field of MOR, where parametric reduced-order models have been constructed with great success across a multitude of fields [44, 89]. This is a relatively mature field whose formal results are quite extensive. For example, UQ for the RBM has been well studied, along with the development of effective algorithms for choosing the best training points [64, 76, 77, 90].

One can formulate the Schrödinger equation in multiple ways, including any flavor of Lippmann-Schwinger (LS) integral equation (which builds in boundary conditions) or as a differential equation in either homogeneous or inhomogeneous form. This freedom, along with the freedom of trial and test bases for the Galerkin projection, leads to multiple alternative emulators that one could construct for quantum scattering systems. For simplicity, we restrict our discussion to two-body scattering for Hermitian Hamiltonians. (See, however, Section 4.6.2 for an extension to higher-body systems.) As a concise reference, we provide Table 1 to show the connections between the fundamental differential or integral equations, variational principles, and Galerkin projections. This section thus provides multiple distinct examples of using Galerkin projections to create emulators, which may prove useful to newcomers wishing to apply model reduction to their own systems, and ends with an example for an emulator applied to a separable potential.

**TABLE 1**. ** **Description of common variational principles (VPs) in quantum scattering, and how to relate them to a Galerkin projection. The quantities are defined as the free wave function |*ϕ*⟩, the full wave function |*ψ*⟩, the scattered wave function |*χ*⟩, and the reactance matrix *K* along with its on-shell form *K*_{E}. Tildes denote trial quantities. The expressions for the Newton VP are written in operator rather than scalar form; any matrix element can be made individually stationary (see Section 4.4 for details). To compute the weak form of the Schwinger and Newton VPs, one must first left multiply by *V*(** θ**) and

*G*

_{0}, respectively, before orthogonalizing against the test basis. The rightmost column specifies whether a constraint for the trial wave function has to be imposed (

*e*.

*g*., using a Lagrange multiplier

*λ*).

### 4.1 Constrained Kohn emulators

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]:

where *K* matrix with on-shell energy *E* = *q*^{2}/2*μ*. This flavor of KVP applies when *ψ* satisfies the asymptotic normalization condition in position space

where *ϕ*(*r*) = *j*_{ℓ}(*qr*) is the (regular) free-space wave function, and *j*_{ℓ}(*qr*) and *n*_{ℓ}(*qr*) are spherical Bessel and Neumann functions, respectively. ^{7} Note that we define the on-shell *K*_{E} matrix as

which differs from the convention in Ref. [30]. The KVP is stationary about exact solutions *ψ*, such that

Eq. 26 can be cast into the form of the generic functional (Eq. 18) by noting that, in position space,

which has defined the surface functional *G* in Eq. 18 and where we have used the Wronskian

Both *rϕ*(*r*) and *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 *ψ*_{i} satisfies Eq. 27, then

whose first term implies that we must impose the constraint *∑*_{i}*β*_{i} = 1.

We are now in a position to find the

where we define *V*_{i} = *V*(*θ*_{i}) and

In the second line we have used that the |*ψ*_{j}⟩ are eigenstates with the corresponding *V*_{j}. If *V*(** θ**) is affine in

**, then**

*θ*Now we can follow the steps outlined in Section 3.1 to determine

where *n*_{b} on-shell *K*-matrices used to train the emulator, and *λ* simply returns the constraint. This system can be solved *via* the system of equations

where *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*

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:

To make explicit the boundary conditions, we make use of the relation

where we have again used the Wronskian from Eq. 30, and assigned *ζ*| instead of |*ψ*⟩. By adding Eqs. 38 and 37, we have

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

and thus, it follows after including the Lagrange multiplier that

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

which is exactly Eq. 34 found by making the KVP stationary. This simplification can be understood by noting that if *ψ*_{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.

### 4.2 Unconstrained Kohn emulators

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* *e*.*g*., explicit substitution of

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

where we used (*via* integration by parts)

We choose our trial function as

Now we can construct the set of linear equations that makes Eq. 43 stationary in *β*_{i}, we find

where

which is the set of equations used to obtain _{ij} can be evaluated with the help of

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 |*ψ*⟩ = |*ϕ*⟩ + |*χ*⟩:

We can construct the weak form by multiplying by a generic test function |*ζ*⟩, which yields

Next, we insert the trial function

We have shown that the coefficients *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*|*ϕ*⟩ = *V*|*ψ*⟩,

with the factors in brackets equating to *β*_{i} using Eq. 45. The equivalence to the KVP is demonstrated in Refs. [94, 95].

### 4.3 Schwinger emulators

The Schwinger variational principle (SVP) is given by [94].

where *G*_{0} is the Green’s operator. This too has the stationary property *ψ* is a wave function satisfying the LS equation. Following the MOR philosophy and inserting a trial function

where

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,

and create a weak form by left-multiplying by *V*(** θ**) along with the test function |

*ζ*⟩:

The weak form can then be converted to its discrete form by setting *ζ*_{i}⟩ = |*ψ*_{i}⟩ for *i* ∈ [1, *…*, *n*_{b}]. ^{8} 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 *K* through

This Equation is exactly the solution for *K* found *via* the LS equation while assuming a finite-rank approximation for *V*:

where

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.

### 4.4 Newton emulators

The Newton variational principle (NVP) for the *K* matrix is given by [31, 96].

where *T*^{(±)} by imposing the associated boundary conditions on *G*_{0}. 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.*,

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*|*ϕ*⟩. 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 *ϕ*′⟩ 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}/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]

with

If the potential *V*(** θ**) has an affine parameter dependence,

*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*, one can insert *δK*)^{2}. The resulting emulator

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

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*_{0} and a test function ⟨*ζ*|:

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

with *M* and *ϕ*′| = ⟨*ϕ*|. Thus, we find the same

Given the optimal coefficients

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.

### 4.5 Origin emulators

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

where *ζ* and *ζ*| = ⟨*ψ*|, but make a Petrov-Galerkin choice for the boundary, with

Thus, the discretized weak form, from which our emulator equations follow, is given by

where we have assumed that the trial basis is constructed such that each snapshot satisfies

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,

with *κ*_{1} = 0.5 and *κ*_{2} = 1. The parameters to be varied are ** θ** = {

*θ*

_{1},

*θ*

_{2}}. 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,

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}.

**FIGURE 8**. Results for the scattering emulator with origin boundary conditions in arbitrary units. Six basis functions are shown as gray lines, and the exact wave function as a black line. Each of the basis functions and the emulated wave function satisfy the constraints at the origin.

### 4.6 General Kohn variational principle

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

^{9}

where

with *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*-matrix in Eq. 72 [32, 88, 98],

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]:

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

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** and

**′ would correspond to**

*u*

*u*_{K}and

*u*_{T}, respectively, given by

Once *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

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.^{10} 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.

#### 4.6.1 Generalization to coupled systems

Following Ref. [88], let us now extend the generalized KVP in Section 4.6 to coupled systems, which could be coupled partial-wave or reaction channels. The stationary approximation to the high-fidelity *L*-matrix then reads

where

with *s* and *s*′ corresponding to the entrance and exit channels. The uncoupled case is retrieved by replacing *ss*′ → *ℓ*. Solving for *ss*’ channels. Note that the coefficients *independently* for each *ss*′ pair.

The coefficients are independent because *ss*’ pair. This becomes apparent when considering how one would solve for the coefficients in the case where there are two uncoupled channels, where *ss*′ → *ℓ* in Eqs. 80 and 81. Here, each partial wave is completely independent of one another, and thus each VP and their corresponding coefficients *ℓ*. Without loss of generality, let the two channels be labeled as *ℓ* = 0 and *ℓ* = 1, and let

An alternative way to understand how the *s*) than the trial functions (*s*′). Because the basis of test functions differs from the basis of the trial function, this is instead a Petrov–Galerkin approach. The linear equations to be solved are exactly what one would obtain from enforcing stationarity in Eq. 80 for each *ss*′ independently. See Ref. [88] for more information on coupled channel emulation.

#### 4.6.2 Generalizations to higher-body systems

The variational emulators for two-body scattering described so far can be generalized to higher-body scattering. In fact, the KVP, as a powerful method for solving scattering problems, has been applied in developing high-fidelity solvers (as opposed to a KVP-based emulator) for studying three- and four-nucleon systems (*e*.*g*., nucleon-deuteron elastic scattering below and above the deuteron break-up threshold) [103, 104].^{11} It is then natural to combine the KVP with the variational emulation strategy to develop fast and accurate emulators beyond just two-body scattering.

Here, we follow Ref. [33], which developed KVP-based emulators for three-body systems. We focus on systems of three identical spinless bosons, particularly the elastic scattering between boson and two-boson bound state in the channel without any relative angular momenta and below the bound state’s break-up threshold. The corresponding scattering *S*-matrix can be estimated *via* a variational functional that resembles Eq. 74 in the two-body case:

Here, *S* is the *S*-matrix associated with the trial three-body wave function *H* and *E* the full Hamiltonian and energy, respectively. The trial wave function has the following asymptotic behavior [33]:

with *R*_{1}, *r*_{1} as one of three different Jacobi coordinate sets; *v* and *P* as the relative velocity and momentum between the scattering particles, *u*_{B}(*r*_{1}) the radial wave function of the two-body bound state.

The emulation procedure is generally similar to those for two-body emulations. We first collect high-fidelity calculations of *i.e.*, training) stage and then use these snapshots as the basis to construct the trial solution (see Eq. 19b) to be used in the variational functional during the online emulation stage. A similar set of the low-dimensional linear equations as in Eq. 35 can be derived to fix the weights *β*_{i}. The variational functional with these inputs produce accurate results for the *S*-matrix at the emulation points [33]. This is all straightforward if we vary only the three-body interactions in *H*(** θ**) when exploring its parameter space. If the two-body interactions are also changed, the two-body bound states of those snapshots are different among themselves and therefore the trial wave functions based on Eq. 19b fails to satisfy the asymptotic behavior described in Eq. 83. In Ref. [33], proper modifications were applied to the constructions of the trial wave functions to satisfy the asymptotic condition. The resulting emulator is again a low-dimensional equation system, but the projected

The results in Ref. [33] are encouraging: the time cost for emulating three-boson scattering is on the order of milliseconds (on a laptop), while the emulation’s relative errors vary from 10^{–13} to 10^{–4} depending on the case. It is straightforward to generalize it to elastic scattering above the break-up threshold, but more studies need to be done for emulating the break-up processes and even higher-body systems. Of course, the Fermi statistics, spin and isospin degrees of freedom, and partial waves beyond the s-wave need to be included to realize emulation for realistic three and higher-body scatterings (*e*.*g*. for three-nucleon systems).

### 4.7 A scattering example

We have covered the reduced-order models that can be constructed from the Kohn, Schwinger, and Newton VPs, and now we put them into action. This example is given in the context of a rank-*n* separable potential where simple analytic forms are available for the snapshots. This provides a sandbox to explore many aspects of the RBM for quantum scattering without the complicating details of more realistic systems. All of the source code that generates the results shown here is available to explore on the companion website [20].

Separable potentials lead to simple formulas for the *K* matrix and the scattering wave function [108]. A rank-*n* separable potential in momentum space is given by

where Λ_{ij} = Λ_{ji} are the coefficients of the potential that will be varied during emulation. For simplicity, we consider here only *s*-wave scattering (*i.e.*, *ℓ* = 0). The potential (Eq. 84) leads to an affine structure that lends itself to the offline-online decomposition discussed in Section 3. From the potential (Eq. 84), simple expressions for *K* and *ψ* can be derived [109]. For instance, the *K* matrix is given in operator form by

with the identity matrix

where the Green’s function *G*_{0} implicitly depends on the on-shell energy *E*. Thus, it follows that *K* is separable if *V* is separable.

We choose to study the Yamaguchi potential [110].

with *ℓ* = 0 and assume a rank-2 potential with *b*_{i} = [2, 4] and 2*μ* = 1. In this case,

which permits all phase shifts, wave functions, and reduced-order matrices (*e*.*g*., _{00}, Λ_{01}, Λ_{11}} are sampled randomly from a uniform distribution in [−50, 50]. The companion website [20] provides the following Python classes that implement the scattering emulators: {`Newton`, `Schwinger`, `Kohn`, `UnconstrainedKohn`}`Emulator`.

Figure 9 shows the phase shifts and the absolute residuals for the Yamaguchi potential Eq. 87 for emulators constructed with *n*_{b} = 2 training points. The top panel depicts the high-fidelity solution (solid black curve) and the emulator results (dots). Here, only the constrained KVP is shown because the others would be indistinguishable. In gray we show the basis states used to train the emulator in the offline stage. The bottom panel shows the absolute residuals for each of the emulators. We can see that all but the constrained KVP are extremely accurate, with the residuals mostly governed by the choice of nugget used to regularize the matrix inversion (see also Section 2.1). For the constrained KVP, we see that the loss of a degree of freedom to implement the constraint significantly impacts its predictive power given that we only have two basis states, although it is still quite accurate in this case. Increasing the basis to *n*_{b} = 5 yields predictions that are accurate to 10^{–13} degrees, or better, for all emulators.

**FIGURE 9**. Phase shifts (top panel) and absolute residuals (bottom panel) in arbitrary units for the Yamaguchi potential Eq. 87 for each scattering emulator discussed above. The solid black lines represent the high-fidelity solution and the dots represent the emulator results. The emulators are so accurate that they are indistinguishable unless looking at residuals. The training set is given by the two gray lines.

Figure 10 shows the high-fidelity (solid black line), emulated (dots), and basis (solid gray line) wave functions for three values of *q* with their absolute residuals using the constrained KVP constructed with *n*_{b} = 5 training points. The emulator reproduces the high-fidelity solution at all three values of *q*, with *q* = 2.0 having the smallest residual. 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]. An example of how the accuracy is affected when varying *n*_{b} is also given in Ref. [88].

**FIGURE 10**. Wave functions (top panel) and absolute residuals (bottom panel) for the Yamaguchi potential Eq. 87 in arbitrary units using the constrained KVP. The top panel legend description is similar to Figure 9, but for three different values of *q*. The bottom panel shows the relative residuals of the three values previously mentioned.

While all emulators described in this Section are applicable to scattering problems in general, their efficacy will depend in practice on various factors, such as their computational complexity and the potential to be emulated. The constrained KVP has the advantage that terms constant in ** θ**, such as the (long-range) Coulomb potential, cancel in the computation of Eq. 33 but it loses one degree of freedom due to the normalization constraint of the coefficients. On the other hand, both the NVP and SVP involve the computation of Green’s functions, which makes them computationally more complex than the KVP—especially the SVP since it also depends quadratically on the potential.

## 5 Summary and outlook

We have presented a pedagogical introduction to projection-based, reduced-order emulators and general MOR concepts suitable for a wide range of applications in low-energy nuclear physics. Emulators are fast surrogate models capable of reliably approximating high-fidelity models due to their reduced content of superfluous information. By making practical otherwise impractical calculations, they can open the door to the various techniques and applications central to the overall theme of this Frontiers Research Topic [14], such as Bayesian parameter estimation for UQ, experimental design, and many more.

In particular, we have discussed variational and Galerkin methods combined with snapshot-based trial (or test) functions as the foundation for constructing fast and accurate emulators. These emulators enable repeated bound state and scattering calculations, *e*.*g*., for sampling a model’s parameter space when high-fidelity calculations are computationally expensive or prohibitively slow. A crucial element in this emulator workflow, as summarized in Figure 2, is an efficient offline-online decomposition in which the heavy computational lifting is performed only once before the emulator is invoked. Chiral Hamiltonians allow for such efficient decompositions due to their affine parameter dependence on the low-energy couplings. Furthermore, we discussed the high efficacy of projection-based emulators in extrapolating results far from the support of the snapshot data, as opposed to the GPs.

While MOR has already reached maturity in other fields, it is still in its infancy in nuclear physics—although rapidly growing—and there remains much to explore and exploit [18, 35, 36, 43]. In the following, we highlight some of the many interesting research avenues for emulator applications in nuclear physics. All of these avenues can benefit from the rich MOR literature and software tools available (*e*.*g*., see Refs. [1–3]):

• Emulator uncertainties need to be robustly quantified and treated on equal footing with other uncertainties in nuclear physics calculations, such as EFT truncation errors. This will be facilitated by the extensive literature on the uncertainties in the RBM [76, 77, 90, 111].

• The performance of competing emulators (*e*.*g*., the Newton vs. Kohn variational principle) is typically highly implementation dependent. Best practices for efficient implementation of nuclear physics emulators should be developed. This may include exploiting MOR software libraries from other fields, such as `pyMOR` [112], when possible.

• Galerkin emulators are equivalent to variational emulators for bound-state and scattering calculations if the test and trial basis are properly chosen. But Galerkin (and especially Petrov-Galerkin) emulators are more general and exploring their applications in non-linear problems may be fruitful in nuclear physics. Emulator applied to non-linear problems will have challenges in terms of both speed and accuracy: 1) the basis size will, in general, need to be large(r) resulting in lower speed-up factors and longer offline stages; 2) using hyperreduction methods will lead to additional approximations that worsen the accuracy of the emulator and whose uncertainties need to be quantified.

• Many technical aspects should be further explored, such as greedy (or active-learning) [24] and SVD-based algorithms for choosing the training points more effectively, hyperreduction methods for non-affine problems, and improved convergence analyses.

• Scattering emulators could play pivotal roles in connecting reaction models and experiments at new-generation rare-isotope facilities (*e*.*g*., the Facility for Rare Isotope Beams). In this regard, further studies on incorporating long-range Coulomb interactions and optical potentials beyond two-body systems will be valuable. Furthermore, emulators for time-dependent density functional theories could see extensive applications in interpreting fission measurements. At facilities such as Jefferson Lab and the future Electron-Ion Collider, explorations of nuclear dynamics at much higher energy scales should also benefit from emulators.

• The emulator framework can be used to extrapolate observables far away from the support of the training data, such as the discrete energy levels of a many-body system calculated in one phase to those of another, as demonstrated in Ref. [22]. Using emulators as a resummation tool to increase the convergence radius of series expansions [26] falls into this category as well, and so does using them to extrapolate finite-box simulations of quantum systems across wide ranges of box sizes [40]. Moreover, for general quantum continuum states, emulation in the complex energy plane can enable computing scattering observables with bound-state methods [113]. Extrapolation capabilities of emulators should be investigated further.

• While projection-based emulators have had successes (*e*.*g*., see Refs. [7, 9, 25]), it is also important to understand their limitations and investigate potential improvements. The synergy between projection-based and machine learning methods [114] is a new direction being undertaken in the field of MOR for this purpose (*e*.*g*., see Ref. [63]). Nuclear physics problems, with and without time dependence, will provide ample opportunities for such explorations.

• Emulators run fast, often with a small memory footprint, and can be easily shared. These properties make emulators effective interfaces for large expensive calculations, through which users can access sophisticated physical models at a minimum cost of computational resources and without specialized expertise, creating a more efficient workflow for nuclear science. As such, emulators can become a collaboration tool [33, 34] that can catalyze new direct and indirect connections between different research areas and enable novel studies.

To help foster the exploration of these (and other) research directions in nuclear physics, we have created a companion website [20] containing interactive supplemental material and source code so that the interested reader can experiment with and extend the examples discussed here.

We look forward to seeing more of the MOR methodology implemented as these research directions are being pursued. But especially we look forward to the exciting applications of emulators in nuclear physics that are currently beyond our grasp.

## Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

## Funding

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under the FRIB Theory Alliance award DE-SC0013617. This work was supported in part by the National Science Foundation under award numbers PHY-1913069 and PHY-2209442 and the NSF CSSI program under award number OAC-2004601 (BAND Collaboration [115]), and the NUCLEI SciDAC Collaboration under U.S. Department of Energy MSU subcontract RC107839-OSU.

## Acknowledgments

We thank the Editors Maria Piarulli, Evgeny Epelbaum, and Christian Forssén for the kind invitation to contribute to this Frontiers Research Topic [14]. We are also grateful to our BUQEYE [19] and BAND [115] collaboration colleagues for sharing their invaluable insights with us, and especially to Pablo Giuliani and Daniel Phillips for fruitful discussions. C.D. thanks Petar Mlinarić for sharing his deep insights into the software library pyMOR [112] for building MOR applications with the Python.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Footnotes

^{1}A *reliable* emulator may not necessarily be required to be highly accurate, *e*.*g*., if the other uncertainties of the theoretical calculation dominate the overall uncertainty budget

^{2}For a critical commentary on the history of the method’s name, see, *e*.*g*., Refs. [45, 46].

^{3}Many helpful theorems relevant to the Rayleigh–Ritz method can be found in Section 3 in Ref. [47].

^{4}In a representation of *H*, the *ψ*_{i} corresponding to |*ψ*_{i}⟩ are the *n*_{b} columns of the matrix *X* in that representation

^{5}For non-Hermitian Hamiltonians, one generally does not obtain the variational bounds (Eq. 7) as can be observed in, *e*.*g*., the subspace-projected coupled-cluster method developed in Ref. [7].

^{6}For simplicity we consider *ξ* to be a real variable; for complex variables, independent variations

^{7}We focus here on examples with real-valued potentials and without long-range Coulomb interactions. Cases with complex-valued potentials and/or the Coulomb interaction may be analyzed in similar ways; relevant discussions specific to Kohn emulators can be found in Refs. [30, 32].

^{8}Note that left-multiplying by *V*(** θ**) and enforcing orthogonality against |

*ζ*

_{i}⟩ = |

*ψ*

_{i}⟩ is different than simply defining |

*ζ*⟩ =

*V*|

*ψ*⟩ and enforcing orthogonality against |

*ζ*

_{i}⟩ =

*V*

_{i}|

*ψ*

_{i}⟩ because

*V*(

**) depends on**

*θ***. Thus, this is indeed a purely Galerkin approach, rather than a Petrov-Galerkin approach**

*θ*^{9}We follow the conventions for scattering matrices in Refs. [99, 100]. Thus, the *K* matrix in this section is dimensionless and defined without the negative sign in Eq. 28.

^{10}These anomalies are not restricted to the KVP but also appear in other VPs such as the NVP and SVP [102].

^{11}The other high-fidelity solvers in this context solve problems in momentum or coordinate space based on the Faddeev formalism [105–107].

## References

1. Benner P, Ohlberger M, Patera A, Rozza G, Urban K. *Reduction of parametrized systems*. Berlin, Germany: Springer (2017).

2. Benner P, Cohen A, Ohlberger M, Willcox K. *Reduction and approximation*. Philadelphia, PA: Society for Industrial and Applied Mathematics: Computational Science & Engineering (2017).

3. Benner P, Gugercin S, Willcox K. A survey of projection-based model reduction methods for parametric dynamical systems. *SIAM Rev* (2015) 57:483–531. doi:10.1137/130932715

4. Zhang X, Nollett KM, Phillips DR. Halo effective field theory constrains the solar ^{7}Be + p → ^{8}B + γ rate. *Phys Lett B* (2015) 751:535–40. arXiv:1507.07239. doi:10.1016/j.physletb.2015.11.005

5. Neufcourt L, Cao Y, Nazarewicz W, Olsen E, Viens F. Neutron drip line in the Ca region from bayesian model averaging. *Phys Rev Lett* (2019) 122:062502. arXiv:1901.07632 [nucl-th]. doi:10.1103/physrevlett.122.062502

6. King GB, Lovell AE, Neufcourt L, Nunes FM. Direct comparison between bayesian and frequentist uncertainty quantification for nuclear reactions. *Phys Rev Lett* (2019) 122:232502. arXiv:1905.05072 [nucl-th]. doi:10.1103/physrevlett.122.232502

7. Ekström A, Hagen G. Global sensitivity analysis of bulk properties of an atomic nucleus. *Phys Rev Lett* (2019) 123:252501. arXiv:1910.02922 [nucl-th]. doi:10.1103/physrevlett.123.252501

8. Catacora-Rios M, King GB, Lovell AE, Nunes FM. Statistical tools for a better optical model. *Phys Rev C* (2021) 104:064611. arXiv:2012.06653 [nucl-th]. doi:10.1103/physrevc.104.064611

9. Wesolowski S, Svensson I, Ekström A, Forssén C, Furnstahl RJ, Melendez JA, et al. Rigorous constraints on three-nucleon forces in chiral effective field theory from fast and accurate calculations of few-body observables. *Phys Rev C* (2021) 104:064001. arXiv:2104.04441 [nucl-th]. doi:10.1103/physrevc.104.064001

10. Svensson I, Ekström A, Forssén C. Bayesian parameter estimation in chiral effective field theory using the Hamiltonian Monte Carlo method. *Phys Rev C* (2022) 105:014004. arXiv:2110.04011 [nucl-th]. doi:10.1103/physrevc.105.014004

11. Odell D, Brune CR, Phillips DR, deBoer RJ, Paneru SN. Performing bayesian analyses with AZURE2 using BRICK: An application to the 7Be system. *Phys Rev C* (2022) 105:014625. arXiv:2105.06541 [nucl-th]. doi:10.3389/fphy.2022.888476

12. Djärv T, Ekström A, Forssén C, Johansson HT. Candidate entanglement invariants for two Dirac spinors. *Phys Rev C* (2022) 105:032402. arXiv:2108.13313 [nucl-th]. doi:10.1103/physreva.105.032402

13. Alnamlah IK, Pérez EAC, Phillips DR. Analyzing rotational bands in odd-mass nuclei using effective field theory and Bayesian methods. *Front Phys* (2022) 10:901954. arXiv:2203.01972 [nucl-th]. doi:10.3389/fphy.2022.901954

14.Frontiers in Physics. Research topic: Uncertainty quantification in nuclear physics (2022) (Online; accessed November 2, 2022)

15. Phillips DR, Furnstahl RJ, Heinz U, Maiti T, Nazarewicz W, Nunes FM, et al. Get on the BAND wagon: A bayesian framework for quantifying model uncertainties in nuclear dynamics. *J Phys G* (2021) 48:072001. arXiv:2012.07704 [nucl-th]. doi:10.1088/1361-6471/abf1df

16. Melendez JA, Furnstahl RJ, Grießhammer HW, McGovern JA, Phillips DR, Pratola MT. Designing optimal experiments: An application to proton compton scattering. *Eur Phys J A* (2021) 57:81. arXiv:2004.11307 [nucl-th]. doi:10.1140/epja/s10050-021-00382-2

17. Farr JN, Meisel Z, Steiner AW. Decision theory for the mass measurements at the facility for rare isotope beams. arXiv:2111.11536 [nucl-th]. doi:10.48550/arXiv.2111.11536

18. Melendez JA, Drischler C, Furnstahl RJ, Garcia AJ, Zhang X. Model reduction methods for nuclear emulators. *J Phys G* (2022) 49:102001. arXiv:2203.05528 [nucl-th]. doi:10.1088/1361-6471/ac83dd

19.BUQEYE collaboration. Buqeye (2021). Available at: https://buqeye.github.io/ (Accessed December 15, 2022).

20.BUQEYE collaboration. Frontiers emulator review (2022). Available at: https://github.com/buqeye/frontiers-emulator-review (Accessed December 15, 2022).

21. Higdon D, McDonnell JD, Schunck N, Sarich J, Wild SM. A Bayesian approach for parameter estimation and prediction using a computationally intensive model. *J Phys G* (2015) 42:034009. arXiv:1407.3017. doi:10.1088/0954-3899/42/3/034009

22. Frame D, He R, Ipsen I, Lee D, Lee D, Rrapaj E. Eigenvector continuation with subspace learning. *Phys Rev Lett* (2018) 121:032501. arXiv:1711.07090. doi:10.1103/physrevlett.121.032501

23. Sarkar A, Lee D. Convergence of eigenvector continuation. *Phys Rev Lett* (2021) 126:032501. arXiv:2004.07651 [nucl-th]. doi:10.1103/physrevlett.126.032501

24. Sarkar A, Lee D. Self-learning emulators and eigenvector continuation. *Phys Rev Res* (2022) 4:023214. arXiv:2107.13449 [nucl-th]. doi:10.1103/physrevresearch.4.023214

25. König S, Ekström A, Hebeler K, Lee D, Schwenk A. Eigenvector continuation as an efficient and accurate emulator for uncertainty quantification. *Phys Lett B* (2020) 810:135814. arXiv:1909.08446 [nucl-th]. doi:10.1016/j.physletb.2020.135814

26. Demol P, Duguet T, Ekström A, Frosini M, Hebeler K, König S, et al. Improved many-body expansions from eigenvector continuation. *Phys Rev C* (2020) 101:041302. arXiv:1911.12578. doi:10.1103/physrevc.101.041302

27. Bai D, Ren Z. Generalizing the calculable R-matrix theory and eigenvector continuation to the incoming-wave boundary condition. *Phys Rev C* (2021) 103:014612. arXiv:2101.06336 [nucl-th]. doi:10.1103/physrevc.103.014612

28. Demol P, Frosini M, Tichai A, Somà V, Duguet T. Bogoliubov many-body perturbation theory under constraint. *Ann Phys* (2021) 424:168358. arXiv:2002.02724 [nucl-th]. doi:10.1016/j.aop.2020.168358

29. Yoshida S, Shimizu N. Constructing approximate shell-model wavefunctions by eigenvector continuation. arXiv:2105.08256 [nucl-th]. doi:10.48550/arXiv.2105.08256

30. Furnstahl RJ, Garcia AJ, Millican PJ, Zhang X. Efficient emulators for scattering using eigenvector continuation. *Phys Lett B* (2020) 809:135719. arXiv:2007.03635 [nucl-th]. doi:10.1016/j.physletb.2020.135719

31. Melendez J, Drischler C, Garcia A, Furnstahl R, Zhang X. Fast and accurate emulation of two-body scattering observables without wave functions. *Phys Lett B* (2021) 821:136608. doi:10.1016/j.physletb.2021.136608

32. Drischler C, Quinonez M, Giuliani PG, Lovell AE, Nunes FM. Toward emulating nuclear reactions using eigenvector continuation. *Phys Lett B* (2021) 823:136777. arXiv:2108.08269 [nucl-th]. doi:10.1016/j.physletb.2021.136777

33. Zhang X, Furnstahl RJ. Fast emulation of quantum three-body scattering. *Phys Rev C* (2021) 105:064004. arXiv:2110.04269 [nucl-th]. doi:10.1103/physrevc.105.064004

34. Tews I, Davoudi Z, Ekström A, Holt JD, Becker K, Briceño R, et al. Nuclear forces for precision nuclear physics – A collection of perspectives. *Few-Body Syst* (2022) 63:67. arXiv:2202.01105. doi:10.1007/s00601-022-01749-x

35. Anderson AL, O’Donnell GL, Piekarewicz J. Applications of reduced-basis methods to the nuclear single-particle spectrum. *Phys Rev C* (2022) 106:L031302. arXiv:2206.14889 [nucl-th]. doi:10.1103/physrevc.106.l031302

36. Giuliani P, Godbey K, Bonilla E, Viens F, Piekarewicz J. Bayes goes fast: Uncertainty quantification for a covariant energy density functional emulated by the reduced basis method. *Front. Phys.* (2022) 13039. arXiv:2209[nucl-th]. doi:10.48550/arXiv.2209.13039

37. Sürer O, Nunes FM, Plumlee M, Wild SM. Uncertainty quantification in breakup reactions. *Phys Rev C* (2022) 106:024607. arXiv:2205.07119 [nucl-th]. doi:10.1103/physrevc.106.024607

38. Bai D, Ren Z. Entanglement generation in few-nucleon scattering. *Phys Rev C* (2022) 106:064005. doi:10.1103/physrevc.106.064005

39. Kravvaris K, Quinlan KR, Quaglioni S, Wendt KA, Navratil P. Quantifying uncertainties in neutron-α scattering with chiral nucleon-nucleon and three-nucleon forces. *Phys Rev C* (2020) 102:024616. arXiv:2004.08474 [nucl-th]. doi:10.1103/physrevc.102.024616

40. Yapa N, König S. Volume extrapolation via eigenvector continuation. *Phys Rev C* (2022) 106:014309. arXiv:2201.08313 [nucl-th]. doi:10.1103/physrevc.106.014309

41. Francis A, Agrawal AA, Howard JH, Kökcü E, Kemper AF. Subspace diagonalization on quantum computers using eigenvector continuation. 10571. arXiv:2209[quant-ph]. doi:10.48550/arXiv.2209.10571

42. Zare A, Wirth R, Haselby CA, Hergert H, Iwen M. Modewise johnson-lindenstrauss embeddings for nuclear many-body theory. 01305. arXiv:2211. doi:10.48550/arXiv.2211.01305

43. Bonilla E, Giuliani P, Godbey K, Lee D. Training and projecting: A reduced basis method emulator for many-body physics. *Phys Rev C* (2022) 106:054322. arXiv:2203.05284 [nucl-th]. doi:10.1103/physrevc.106.054322

44. Benner P, Schilders W, Grivet-Talocia S, Quarteroni A, Rozza G, Miguel Silveira L. *System- and data-driven methods and algorithms*. Berlin, Germany: De Gruyter (2021).

45. Leissa A. The historical bases of the Rayleigh and Ritz methods. *J Sound Vibration* (2005) 287:961–78. doi:10.1016/j.jsv.2004.12.021

46. Ilanko S. Comments on the historical bases of the Rayleigh and Ritz methods. *J Sound Vibration* (2009) 319:731–3. doi:10.1016/j.jsv.2008.06.001

47. Suzuki Y, Varga K. *Stochastic variational approach to quantum-mechanical few-body problems*. Berlin, Germany: Springer Berlin Heidelberg (1998).

48. Benner P, Schilders W, Grivet-Talocia S, Quarteroni A, Rozza G, Miguel Silveira L. *Order reduction: Volume 2: Snapshot-based methods and algorithms*. Berlin, Germany: De Gruyter (2020). p. 1–348.

49. Buchan AG, Pain CC, Fang F, Navon IM. A POD reduced-order model for eigenvalue problems with application to reactor physics. *Int J Numer Methods Eng* (2013) 95:1011–32. doi:10.1002/nme.4533

50. Quarteroni A, Manzoni A, Negri F. Reduced basis methods for partial differential equations. In: *An Introduction, La Matematica per il 3+2*. Berlin, Germany: Springer International Publishing (2016). p. 92.

51. Franzke MC, Tichai A, Hebeler K, Schwenk A. Excited states from eigenvector continuation: The anharmonic oscillator. *Phys Lett B* (2022) 830:137101. arXiv:2108.02824. doi:10.1016/j.physletb.2022.137101

52. Melendez J. *Effective field theory truncation errors and why they matter*. Ph.D. thesis. Columbus, Ohio: Ohio State U. (2020).

53. Hicks C, Lee D. Trimmed sampling algorithm for the noisy generalized eigenvalue problem. 02083. arXiv:2209. doi:10.48550/arXiv.2209.02083

54. Neumaier A. Solving ill-conditioned and singular linear systems: A tutorial on regularization. *SIAM Rev* (1998) 40:636–66. doi:10.1137/s0036144597321909

55. Engl H, Hanke M, Neubauer A. Regularization of inverse problems. In: *Mathematics and its applications*. Netherlands: Springer (1996).

56. Hergert H. A guided tour of *ab initio* nuclear many-body theory. *Front Phys* (2020) 8:37905061. arXiv:2008. doi:10.3389/fphy.2020.00379

57. Ghattas O, Willcox K. Learning physics-based models from data: Perspectives from inverse problems and model reduction. *Acta Numerica* (2021) 30:445–554. doi:10.1017/s0962492921000064

58. Rasmussen CE, Williams CKI. *Gaussian processes for machine learning, adaptive computation and machine learning series*. Cambridge, MA: University Press Group Limited (2006).

59. Mezić I. Analysis of fluid flows via spectral properties of the koopman operator. *Annu Rev Fluid Mech* (2013) 45:357–78. doi:10.1146/annurev-fluid-011212-140652

60. Kutz JN, Brunton SL, Brunton BW, Proctor JL. Dynamic mode decomposition. In: *Other titles in applied mathematics*. Pennsylvania, United States: SIAM-Society for Industrial and Applied Mathematics (2016).

61. Raissi M, Perdikaris P, Karniadakis GE. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *J Comput Phys* (2019) 378:686–707. doi:10.1016/j.jcp.2018.10.045

62. Chen W, Wang Q, Hesthaven JS, Zhang C. Physics-informed machine learning for reduced-order modeling of nonlinear problems. *J Comput Phys* (2021) 446:110666. doi:10.1016/j.jcp.2021.110666

63. Fresca S, Manzoni A. POD-DL-ROM: Enhancing deep learning-based reduced order models for nonlinear parametrized PDEs by proper orthogonal decomposition. *Comp Methods Appl Mech Eng* (2022) 388:114181. doi:10.1016/j.cma.2021.114181

64. Hesthaven J, Rozza G, Stamm B. Certified reduced basis methods for parametrized partial differential equations. In: *SpringerBriefs in mathematics*. Berlin, Germany: Springer International Publishing (2015).

65. Barrault M, Maday Y, Nguyen NC, Patera AT. An ‘empirical interpolation’ method: Application to efficient reduced-basis discretization of partial differential equations. *Comptes Rendus Mathematique* (2004) 339:667–72. doi:10.1016/j.crma.2004.08.006

66. Grepl MA, Maday Y, Nguyen NC, Patera AT. Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations. *ESAIM: M2AN* (2007) 41:575–605. doi:10.1051/m2an:2007031

67. Horacio LC, Solis-Daun J. Synthesis of positive controls for the global CLF stabilization of systems. In: Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference; 15-18 December 2009; Shanghai, China (2009). p. 4316–21.

68. Chaturantabut S, Sorensen DC. Nonlinear model reduction via discrete empirical interpolation. *SIAM J Scientific Comput* (2010) 32:2737–64. doi:10.1137/090766498

69. Bui-Thanh T, Damodaran M, Willcox K. Proper orthogonal decomposition extensions for parametric applications in compressible aerodynamics. In: 21st AIAA Applied Aerodynamics Conference; 23 June 2003 - 26 June 2003; Orlando, Florida (2003).

70. Carlberg K, Bou-Mosleh C, Farhat C. Efficient non-linear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations. *Int J Numer Methods Eng* (2011) 86:155–81. doi:10.1002/nme.3050

71. Amsallem D, Cortial J, Farhat C. Towards real-time computational-fluid-dynamics-based aeroelastic computations using a database of reduced-order information. *AIAA J* (2010) 48:2029–37. doi:10.2514/1.j050233

72. An SS, Kim T, James DL. Session details: Character animation II. *ACM Trans Graph* (2008) 27:3262975. doi:10.1145/3262975

73. Farhat C, Avery P, Chapman T, Cortial J. Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency. *Int J Numer Methods Eng* (2014) 98:625–62. doi:10.1002/nme.4668

74. Guo M, Hesthaven JS. Reduced order modeling for nonlinear structural analysis using Gaussian process regression. *Comp Methods Appl Mech Eng* (2018) 341:807–26. doi:10.1016/j.cma.2018.07.017

75. Gubisch M, Volkwein S. Chapter 1: Proper orthogonal decomposition for linear-quadratic optimal control. In: *Reduction and approximation* (2017). p. 3–63.

76. Rozza G, Huynh DBP, Patera AT. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations. *Arch Comput Methods Eng* (2008) 15:229–75. doi:10.1007/s11831-008-9019-9

77. Chen P, Quarteroni A, Rozza G. Reduced basis methods for uncertainty quantification. *SIAM/ASA J Uncertainty Quantification* (2017) 5:813–69. doi:10.1137/151004550

78. Horger T, Wohlmuth B, Thomas D. Simultaneous reduced basis approximation of parameterized elliptic eigenvalue problems. *ESAIM: M2AN* (2017) 51:443–65. doi:10.1051/m2an/2016025

79. MacKay DJC. Neural networks and machine learning. In: CM Bishop, editor. *NATO ASI series*. Berlin: Springer (1998). p. 133–66.

80. MacKay DJC. *Information theory, inference, and learning algorithms*. Cambridge, UK: Cambridge University Press (2003).

81. Gander MJ, Wanner G. From euler, Ritz, and Galerkin to modern computing. *SIAM Rev* (2012) 54:627–66. doi:10.1137/100804036

82. Zienkiewicz O, Taylor R, Zhu J. *The finite element method: Its basis and fundamentals*. 7th ed. Oxford: Butterworth-Heinemann (Butterworth-Heinemann (2013).

83. Zienkiewicz O, Taylor R, Fox D. *The finite element method for solid and structural mechanics*. 7th ed. Oxford: Butterworth-Heinemann (2014).

84. Zienkiewicz O, Taylor R, Nithiarasu P. *The finite element method for fluid dynamics*. 7th ed. Oxford: Butterworth-Heinemann (2014).

85. Mikhlin SG. *Variational methods in mathematical physics*. Oxford: Cambridge University Press (1967). [Translated by T. Boddington for Pergamon Press, Oxford, 1964].

86. Evans JD. *Straightforward statistics for the behavioral sciences*. Pacific Grove, CA: Brooks/Cole Publishing (1996).

87. Brenner SC, Scott LR. The mathematical theory of finite element methods. In: *Texts in applied mathematics*. Berlin, Germany: Springer (2008).

88. Garcia AJ, Drischler C, Furnstahl RJ, Melendez JA, Zhang X. Wave function-based emulation for nucleon-nucleon scattering in momentum space. arXiv:2301.05093.

89. Benner P, Schilders W, Grivet-Talocia S, Quarteroni A, Rozza G, Miguel Silveira L. *Order reduction: Volume 3: Applications*. Berlin, Germany: De Gruyter (2020).

90. Huynh D, Rozza G, Sen S, Patera A. A successive constraint linear optimization method for lower bounds of parametric coercivity and inf–sup stability constants. *Comptes Rendus Mathematique* (2007) 345:473–8. doi:10.1016/j.crma.2007.09.019

91. Kohn W. Variational methods in nuclear collision problems. *Phys Rev* (1948) 74:1763–72. doi:10.1103/physrev.74.1763

92. Kohn W. Variational scattering theory in momentum space I. Central field problems. *Phys Rev* (1951) 84:495–501. doi:10.1103/physrev.84.495

93. Drischler C, Holt JW, Wellenhofer C. Chiral effective field theory and the high-density nuclear equation of state. *Annu Rev Nucl Part Sci* (2021) 71:403–32. arXiv:2101.01709. doi:10.1146/annurev-nucl-102419-041903

94. Takatsuka K, Lucchese RR, McKoy V. Relationship between the schwinger and kohn-type variational principles in scattering theory. *Phys Rev A* (1981) 24:1812–6. doi:10.1103/physreva.24.1812

95. Takatsuka K, McKoy V. Variational scattering theory using a functional of fractional form. I. General theory. *Phys Rev A* (1981) 23:2352–8. doi:10.1103/physreva.23.2352

97. Thompson IJ, Nunes FM. *Nuclear reactions for astrophysics: Principles, calculation and applications of low-energy reactions*. Cambridge, United Kingdom: Cambridge University Press (2009).

98. Lucchese RR. Anomalous singularities in the complex kohn variational principle of quantum scattering theory. *Phys Rev A* (1989) 40:6879–85. doi:10.1103/physreva.40.6879

99. Taylor JR. *Scattering theory: The quantum theory of nonrelativistic collisions*. Illinois, United States: Dover (2006).

100. Morrison MA, Feldt AN. Through scattering theory with gun and camera: Coping with conventions in collision theory. *Am J Phys* (2007) 75:67–80. doi:10.1119/1.2358156

101. Nesbet R. *Variational methods in electron-atom scattering theory, Physics of atoms and molecules*. Plenum Press (1980).

102. Adhikari SK. Anomalies of variational phase shifts. *Chem Phys Lett* (1991) 181:435–40. doi:10.1016/0009-2614(91)90376-k

103. Marcucci LE, Dohet-Eraly J, Girlanda L, Gnech A, Kievsky A, Viviani M. The hyperspherical harmonics method: A tool for testing and improving nuclear interaction models. *Front Phys* (2020) 8:69. arXiv:1912.09751 [nucl-th]. doi:10.3389/fphy.2020.00069

104. Kievsky A, Rosati S, Viviani M, Marcucci LE, Girlanda L. A high-precision variational approach to three- and four-nucleon bound and zero-energy scattering states. *J Phys G* (2008) 35:063101. arXiv:0805.4688. doi:10.1088/0954-3899/35/6/063101

105. Gloeckle W, Witala H, Huber D, Kamada H, Golak J. A New look into the partial wave decomposition of three nucleon forces. *Phys Rept* (1996) 274:107. doi:10.1007/s006010050057

106. Deltuva A, Fonseca AC, Lazauskas R. Faddeev equation approach for three-cluster nuclear reactions. In: *Lecture notes in physics*. Berlin, Germany: Springer (2014). p. 1. arXiv:1201.4979 [nucl-th].

107. Lazauskas R, Carbonell J. The faddeev–yakubovsky symphony. *Few Body Syst* (2019) 60:62. arXiv:1908.04861 [quant-ph]. doi:10.1007/s00601-019-1529-5

108. Tabakin F. Inverse scattering problem for separable potentials. *Phys Rev* (1969) 177:1443–51. doi:10.1103/physrev.177.1443

109. Kwong NH, Köhler HS. SeparableNNpotentials from inverse scattering for nuclear matter studies. *Phys Rev C* (1997) 55:1650–64. doi:10.1103/physrevc.55.1650

110. Göbel M, Hammer HW, Ji C, Phillips DR. Momentum-space probability density of {}^6He in halo effective field theory. *Few Body Syst* (2019) 60:61. arXiv:1904.07182 [nucl-th]. doi:10.1007/s00601-019-1528-6

111. Haasdonk B. Chapter 2: Reduced basis methods for parametrized pdes—A tutorial introduction for stationary and instationary problems. In: *Reduction and approximation*. Philadelphia, PA: Computational Science & Engineering (2017). p. 65–136.

112. Milk R, Rave S, Schindler F. pyMOR -- generic algorithms and interfaces for model order reduction. *SIAM J Scientific Comput* (2016) 38:S194–216. doi:10.1137/15m1026614

113.Zhang X, Low Energy Community Meeting. *Recent developments in the emulations of quantum continuum states*. USA: Argonne National Laboratory (2022).

114. Boehnlein A, Diefenthaler M, Sato N, Schram M, Ziegler V, Fanelli C, et al. *Colloquium*: Machine learning in nuclear physics. *Rev Mod Phys* (2022) 94:031003. arXiv:2112.02309 [nucl-th]. doi:10.1103/revmodphys.94.031003

115.BAND. Bayesian analysis of nuclear dynamics (BAND) framework project (2020). Available at: https://bandframework.github.io/ (Accessed December 15, 2022).

Keywords: emulators, reduced-order models, model order reduction, nuclear scattering, uncertainty quantification, effective field theory, variational principles, Galerkin projection

Citation: Drischler C, Melendez JA, Furnstahl RJ, Garcia AJ and Zhang X (2023) BUQEYE guide to projection-based emulators in nuclear physics. *Front. Phys.* 10:1092931. doi: 10.3389/fphy.2022.1092931

Received: 08 November 2022; Accepted: 15 December 2022;

Published: 17 February 2023.

Edited by:

Maria Piarulli, Washington University in St. Louis, United StatesCopyright © 2023 Drischler, Melendez, Furnstahl, Garcia and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: C. Drischler, drischler@ohio.edu