Variations on the Maiani-Testa approach and the inverse problem

We discuss a method to construct hadronic scattering and decay amplitudes from Euclidean correlators, by combining the approach of a regulated inverse Laplace transform with the work of Maiani and Testa. Revisiting the original result, we observe that the key observation, i.e. that only threshold scattering information can be extracted at large separations, can be understood by interpreting the correlator as a spectral function, $\rho(\omega)$, convoluted with the Euclidean kernel, $e^{- \omega t}$, which is sharply peaked at threshold. We therefore consider a modification in which a smooth step function, equal to one above a target energy, is inserted in the spectral decomposition. This can be achieved either through Backus-Gilbert-like methods or more directly using the variational approach. The result is a shifted resolution function, such that the large $t$ limit projects onto scattering or decay amplitudes above threshold. The utility of this method is highlighted through large $t$ expansions of both three- and four-point functions that include leading terms proportional to the real and imaginary parts (separately) of the target observable. This work also presents new results relevant for the un-modified correlator at threshold, including expressions for extracting the $N \pi$ scattering length from four-point functions and a new strategy to organize the large $t$ expansion that exhibits better convergence than the expansion in powers of $1/t$.


Introduction
Over the last decades, numerical lattice QCD has become a high-precision tool for predicting several non-perturbative strong-force observables, including hadronic masses, decay constants and form factors. Looking beyond these quantities, each defined in terms of single-hadron states, lattice QCD has also shown outstanding progress in calculating multi-hadron observables including 2 → 2 scattering and 1 → 2 decay amplitudes. The successful determination of such amplitudes is especially impressive due to the presence of the Euclidean metric, required in order to apply Monte Carlo methods in estimating the lattice QCD path integral. While the analytic continuation to Minkowski correlators is formally guaranteed by the Osterwalder-Schrader theorem [2], in practice the limited knowledge of correlation functions on a finite set of points, together with the presence of statistical uncertainties, implies that direct extraction is a numerically ill-posed inverse problem, see e.g. ref. [3].
Offering another perspective on this challenge, Maiani and Testa [1] showed that (for energies above production threshold) asymptotically separating individual pion fields in Euclidean time leads to correlators dominated by off-shell contributions. As we review in more detail in section 2, this means that matrix elements of the form π|π(0)|ππ contribute, where π| represents a single pion state and |ππ a two-pion asymptotic in-or out-state. The infamous off-shellness refers to the fact that the four momentum associated with the operatorπ does not satisfy q 2 = m 2 π , where m π is the physical particle mass. As a result the matrix element depends on the details of the operator used and gives no useful information.
This limitation was circumvented by Lüscher [4,5], who showed that the values of the low-lying two-pion energies in a finite periodic spatial volume are determined by the on-shell 2 → 2 scattering amplitude (up to corrections falling faster than any power of 1/L where L denotes the box length). Thus, the scattering amplitude can be extracted indirectly, from discrete spectra determined in lattice calculations. In more recent years, due especially to efficient methods for evaluating quarkfield Wick contractions [6] and an improved understanding of the importance of a large operator basis [7], this method has proven to be extremely successful in the determination of elastic twohadron scattering amplitudes. 1 Formal extensions of the relations between amplitudes and energies  have allowed the same basic approach to be applied in calculations involving particles with spin, coupled-channel two-particle systems, and most recently to the scattering of three-pion states.
The main limitation of the finite-volume formalism, besides the technical challenge of measuring several excited-state energy levels and matrix elements, is the proliferation of multi-particle channels as the scattering energy increases. Here it is important to note that any finite-volume method must formally treat all open multi-hadron channels (or argue that they are irrelevant) in order to reach a prediction about any individual scattering, decay, or transition amplitude. This is because finitevolume energies and matrix elements generically depend on a mixture of all physically allowed scattering processes.
For these reasons, recent work has revisited prospects for the direct analytic continuation of numerical correlation functions in the context of inverting the Laplace transform where G(t, L) represents a general Euclidean correlator and the right-hand side defines the spectral function, ρ(ω, L). A regulated inverse of this simple relation could potentially unlock an alternative method in extracting inclusive quantities, such as heavy-particle lifetimes or the hadronic tensor [59], as well as scattering amplitudes [60], directly from spectral functions convoluted with a known resolution function where δ ∆ (ω, ω) is peaked at ω = ω with characteristic width ∆. Encouraging progress has recently been made in developing improved strategies to regulate the inverse problem and systematically target a resolution function [61,62] to extract eq. (1.2) from eq. (1.1). To reach a physical prediction, these methods formally require that the infinite-volume limit (L → ∞) is taken before the resolution width is sent to zero (∆ → 0). Such ideas might prove useful also in the context of QED corrections to semi-leptonic decays, and similar processes, 2 where intermediate on-shell states prevent the analytic continuation to Minkowski signature and approximate numerical solutions to the inverse problem could play a significant role.
In this manuscript, we present a strategy for combining the ideas summarized by eqs. (1.1) and (1.2) with the work of Maiani and Testa [1], in order to extend the reach of the latter to energies above scattering threshold, without suffering the dominance of off-shell terms. We present the idea both in the context of three-and four-point functions and, as a side benefit, we also reach new results for extracting threshold information from standard correlators.
Although we focus in this work on three-and four-point functions, the basic idea can already be expressed with a two point function of a scalar current, J(t, x), where we assume the L → ∞ limit has been taken. We then define the modified correlator where Θ(z, ∆) is a smoothened Heaviside function, interpolating from zero for z < 0 to one for z > 0. A specific definition of Θ(z, ∆) is given in eq. (2.5) below, but any function can be used provided it is smooth and becomes the usual Heaviside step function for ∆ → 0. Defining the spectral function as note that the following relations hold: These two simple results form the basis of this work. In the expression for G(t), the kernel e −ωt becomes sharply peaked at threshold for large t, leading to the threshold dominance famously identified by Maiani and Testa [1]. For G Θ (t|s), by contrast, the peak forms at ω = √ s and one can extract time-like observables away from threshold. 3 See also figure 1 below.
To give an impression of the relations derived in the following, we close this introduction by highlighting two of our key results.
First, returning to the scalar current, J(x), we introduce the form factor where |0 is the vacuum and s, ππ, out| is a two-particle out state with squared center-of-momentum energy, s, projected to zero angular momentum by the current. We then demonstrate in section 2 that the following holds . (1.8) Here the left-hand side is a product of a simple normalization factor, N , with an infinite-volume Euclidean matrix element, modified by the smooth Heaviside-function Θ(z, ∆). The matrix element is built from a single-pion state and a generic operator π −q (0) such that 0| π −q (0) has the quantum numbers of π, −q|. 4 The right-hand side of eq. (1.8) demonstrates the utility of this quantity. Specifically, as we prove in section 2, it is equal to a known linear combination (with one time-independent coefficient Θ(0, ∆) and a second time-dependent function J (0) ), of the real and imaginary parts of the target observable, f (s), up to terms that are suppressed for well-chosen values of t and ∆. The suppression is quantified through an asymptotic series of known geometric functions that can be understood as a generalization of the large t expansion of Maiani and Testa [1]. As with the inverse techniques of refs. [59,60], this result holds for all s and may be competitive with standard methods for s > (4m π ) 2 (m π is the pion mass), where it is challenging to disentangle the multiple open channels from finite-volume information.
A simple cross check of eq. (1.8) is given by setting s = 4m 2 π and carefully taking the ∆ → 0 limit, as detailed in section 2.1. Then one recovers the original result of ref. [1] N π, 0| π 0 (t) In section 2.1 we also discuss the finite t corrections to this and revisit the expansion performed by Maiani and Testa. We express the result in an alternative basis that exhibits faster convergence and give expressions for the next-to-leading order term (depending on the ππ scattering length) and the next-to-next-to-leading order term (depending on the derivative of the scattering amplitude with respect to the virtuality of an off-shell leg).
Second, we consider the same derivation with four-point functions, studying both the standard correlation function and the modification with Θ. In direct analogy to eq. (1.8) above, we deduce that one can extract the 2 → 2 scattering amplitude at all energies from the Θ-correlator. Also of great interest, however, is the standard correlation function at two-particle threshold, which is more straight-forward to implement in the short term. For example in the case of N π scattering, we derive the following result in section 3 Here the left-hand side is a Euclidean matrix element that can be extracted from a four-point function of pion and nucleon fields, all projected to zero spatial momentum. The indices a , a encode isospin and b , b simultaneously indicate isospin and spin such that the coefficients C ab can be chosen to project onto any definite N π quantum numbers (with an implicit sum over indices understood). The left-hand side also includes a simple normalization factor N 2 and the subscript c indicates that the disconnected N → N and π → π contribution is subtracted.
In analogy to the Maiani-Testa correlation function at threshold, this four-point function can be used to directly extract the N π scattering length, a N π , in the channel specified by C ab . The advantage compared to the three-point function is that one recovers a function only of a N π as opposed to a combination of the threshold form-factor [f (4m 2 π )] and the scattering length. Indeed the right-hand side contains two terms, scaling as ta N π and √ ta 2 N π so that a strong constraint on the scattering length may be achieved from the combined fit.
The remainder of this manuscript is organized as follows: In the next section we revisit the derivation of ref. [1] and show how the result is modified for the Θ-correlator. We also present a new form of the large t expansion, in terms of a series of known integrals. We then repeat the analysis for four-point functions in section 3, where we also present the derivation of eq. (1.10). In section 4, we discuss implementation strategies to extract the Θ-correlator, including the role of finite-volume effects. We also include a number of appendices, deriving various technical results used in the main text.
As was shown in ref. [1], for q 1 = q 2 [which, for vanishing total momentum (q 1 + q 2 = 0), translates to q 1 = 0], the quantity c q1q2 (t 2 ) contains growing exponentials as t 2 → ∞, with coefficients depending on scattering amplitudes for which one of the external legs is off the mass shell. To better understand the limitations associated with this, we define the modified quantity where Θ(ω, ∆) is a smooth Heaviside function, e.g.
Θ(ω, ∆) = tanh(2ω/∆) + 1 2 , (2.5) andM = Ĥ2 −P 2 is an operator giving the center-of-momentum-frame energy. In section 4, we describe strategies for accessing this object in a numerical lattice calculation, including the role of the finite-volume boundary conditions. In this section we take c Θ as given and show that it can be used to access the real and imaginary parts of where we have introduced In the same spirit as the inverse methods described in refs. [59,60], this approach is valid for all values of s(q), provided the correlator can be estimated with an appropriate range of ∆ and ω 0 .
Before deriving our main result, we pause here to discuss the intuition motivating eq. (2.4) and its relation to the results of Maiani and Testa. For this discussion we set q 1 = −q 2 = q and ω 0 = ω q .
The key point is that c q,−q (t 2 ) and c Θ q,−q (t 2 |ω q ) can be written as the convolution of the same spectral function, ρ, with two different resolution functions, δ and δ Θ . In particular , (2.9) Figure 1: Sketch of the normalized effective resolution functions, δ(E, t 2 ) and δ Θ (E, t 2 ), defined in eqs. (2.11) and (2.12) respectively. The corresponding correlators can be written as convolution integrals of these resolution functions with the same spectral function, ρ q (E), which contains the time-like information we are after. Thus, the functional forms plotted here give an indication of the energies that predominantly contribute to the correlation functions. The left panel illustrates that c q,−q (t 2 ) is dominated by E ≈ 2m π and that large t 2 sharpens the near-threshold resolution. By contrast, the right panel illustrates the modification involving Θ(E − 2ω q , ∆) (also plotted here in gray) for 2ω q = 5m π . The three curves correspond to the same ∆ value, and the sharpening around E = 5m π is achieved by increasing t 2 , with the dashed curves showing the corresponding factors of e −t2(E−2ωq) . where Here we have included a zero-width Heaviside function (denoted θ(x)) in the definition of δ(E, t 2 ) to further emphasize the similarities. This is allowed as the spectral function, ρ(E), has zero support for E < 2m π . 5 In figure 1 we show the functional forms of δ(E, t 2 ) and δ Θ (E, t 2 ), normalized to unit area. Note that δ(E, t 2 ) is sharply peaked at threshold with a width given by 1/t 2 . For this reason, the large t 2 limit can only access time-like information at threshold, as was famously established in ref. [1]. For δ Θ (E, t 2 ), the peak is shifted and mimics the resolution functions discussed in refs. [59,60], formally allowing one to extract scattering amplitudes at all energies. The key advantage, as compared to the earlier work, is that the effective resolution width can be reduced at fixed ∆ by increasing t 2 . This gives a powerful handle on the target observable and can be expressed as a large t 2 expansion, to which we now turn.
Returning to eq. (2.4), the next step is to insert a complete set of states between π q2 (0) and J(0) to reach 5 We assume throughout that the system has no bound state. where and s(p) = E(p) 2 − P 2 , with the total energy and momentum given by The complete set of states leads to the matrix elements where k is a channel index (referring e.g. to ππ, KK, ππππ), n k is the number of particles in the channel and S k the number of permutations of identical particles. We define k = 1 as the ππ channel, with internal quantum numbers matching π(0)|π, q .
The final key step is to express A k (q 1 , q 2 ; p) in terms of an off-shell amplitude M k (η|µ(q, p)). We use the bold notation here whenever the amplitude has at least one off-shell external leg. This object is defined in more detail in appendix A.1 and for the present argument it suffices to note that M k (0|µ(q, p)) = M k (µ(q, p)) is the usual on-shell two-to-n k scattering amplitude as a function of Lorentz invariants organized in the vector µ(q, p). For example in the case of a two-to-two amplitude The matrix element A k (q 1 , q 2 ; p) admits a simple decomposition in terms of the off-shell amplitude where the distance of the π(0) leg from the mass shell is given by Combining Eqs. (2.13) and (2.20) gives the generalization of Maiani and Testa's result to the case of non-zero total momenta and to the Θ-function modification of the correlator. In appendix A.2 we give explicit expressions in the case of general kinematics and discuss their utility. Here we focus on the role of the Heaviside function and set q 1 = −q 2 = q to reach where we have replaced q → q in quantities that only depend on q 2 = q 2 and have introduced For example, in the case where only the ππ channel contributes, i.e. for E 2 < (4m π ) 2 , this becomes where M s is the scattering amplitude projected to zero orbital angular momentum The on-shell limit is given by Remarkably, the relation between the on-shell G function and Imf (E 2 ) in fact holds for all E 2 , as we prove in appendix A.3.
Returning to the general case, the final step is to perform the change of variables and P indicates the principal-value pole prescription. To reach this expression, one uses the fact that c Θ is real so that one can take the real part of the right-hand side for free. As G is also real this simply replaces: Alternatively, one can demonstrate that the imaginary part of 1/(η − i ) explicitly cancels that of f (s(q 2 )), as we review in appendix A.4, following ref. [1].
Equations (2.23), (2.29) and (2.30) summarize the main new technical results of this section. In the following subsections we explore their implications for extracting scattering information.

Threshold kinematics
As a first step we set ω 0 = ω q − = m π − to reach the threshold case already considered in ref. [1]. For these choices (together with ∆ → 0) the Θ-function has no effect and only the two-pion channel  contributes in G. We reach 31) or, to give an explicit expression in terms of the Euclidean correlator, where e −3mπt1 has been dropped. Here the g n coefficients arise in an expansion of the amplitudes 33) and the functions I (n) (b) are integrals of the kernel κ times the phase space, E 2 /4 − m 2 π /(8πE), both in the limiting case of the threshold amplitude. They can be expressed as follows Expanding in powers of 1/b = 1/(2m π t 2 ) generates the original series presented in ref. [1]. However, as we illustrate in figure 2, the leading terms converge slowly to the full integrals and the latter may prove to be a better basis in describing the correlation function.
Combining g 0 = M s (4m 2 π ) = 32πm π a ππ with Γ(1/2) = √ π gives the known result 6 and thus where a ππ is the two-particle scattering length. Our expressions allow one to easily reach generalizations of this, for example This gives a more explicit expression for the leading off-shellness contaminating this correlator. The η derivative represents a small variation in the virtuality of one external leg away from the mass shell. The quantity is perfectly well defined but depends on the details of the operator π(x) and has no physical meaning for pion scattering.

General kinematics
Next we take ω 0 = ω q > m π and write where Θ(0, ∆) = 1/2 for the specific choice given in eq. (2.5). For the second term we have substituted the expansion and have introduced the basis of integrals Note here that the principal value is only required for n = 0.
Here we differ from eq. (10) of ref. [1] by a factor of 2 in the 1/ √ t 2 -dependent term. We are currently in correspondence with the authors of that work to clarify this difference. Note also that, following ref. [1], we have set our convention for the scattering-length such that aππ < 0 for repulsive interactions at threshold. As with the threshold case, the coefficients g n for n > 1 contain off-shell derivatives of the scattering amplitude.
[See e.g. eq. (2.39).] In figure 3 we plot the functions J (n) (t 2 , ω q , ∆) for fixed values of ω q and ∆, as a function of t 2 . We stress that these expressions require ∆ > 0 in order for the principal-value pole prescription to properly regulate the integral.
We close this subsection by addressing a final subtlety concerning the fact that the sum over J (n) in eq. (2.43) [as well as that over I (n) in eq. (2.31)] is, in fact, a divergent asymptotic series. The detailed properties of this series depend on the coefficients g n , but the divergence is expected as J (n) ∼ n! for large n. We stress however that the task here is to estimate the difference between c Θ and the first two terms on the right-hand side of eq. (2.43). This difference is finite and the divergent asymptotic series provides a representation of it in the usual way, i.e.
Given the clear hierarchy in the J (n) functions, it should be feasible to keep the first few terms, with the g n as fit parameters and remove the contamination in order to extract the form factor.
3 Two-to-two scattering amplitudes We now imitate the results of the previous section, but focusing here on a four-point function built from two nucleon and two pion fields. We begin with the four-point analog of eq. (2.4), defining Here a and a are isospin indices for the pions and b and b are combined spin and isospin indices. The coefficients C ab and C a b are chosen to project to N π states with definite spin and isospin. As above, here a sum over the repeated indices is implied and the tildes on the operators indicate the projection to definite spatial momentum. Note that momenta with a 1 subscript correspond to pion fields and those with a 2 to nucleons, such that 2) where m N is the nucleon mass. Up to the Θ function, this quantity is directly extractable from a standard Euclidean correlator The strategy for implementing Θ(M − M 0 , ∆) is discussed in the next section. 7 Inserting a complete set of states then yields in direct analog to eq. (2.13) above, here with where k = 1 represents the N π channel projected by C ab and only a single delta function arises, since only a single contraction can appear in the disconnected contribution. Note that this decomposition assumes that no two-particle bound state appears in the channel of interest, an assumption we have already noted in footnote 3. If a bound state is present, this can be subtracted from the correlator to reach a new quantity, for which the decomposition here still applies.
Substituting eq. (3.6) into eq. (3.4) yields four terms: one doubly disconnected, two singly disconnected, and the final term proportional to |M| 2 . Here we will assume the doubly disconnected piece is subtracted, using a subscript c to denote the resulting connected correlator. As we prove in appendix A.5, simplifying the remaining terms and setting q 1 = −q 2 = q and q 1 = −q 2 = q then gives

7)
7 Here we denote the inflection point of Θ by M 0 rather than 2ω 0 , which is more convenient since the two hadrons are non-degenerate.
where we have suppressed the k = 1 index, introduced mass labels on the ωs (since the momenta no longer give a distinction) and also introduced E(q) = ω π,q + ω N,q and η(q , q) ≡ E(q) − ω π,q 2 − (ω N,q ) 2 , (3.8) In eq. (3.7) we have also used the fact that the imaginary parts cancel between the various contributions, as we demonstrate in appendix A.5.
In the next paragraph we discuss the final term in eq. (3.7), C(t , t|q , q, M 0 , ∆), a direct analog of the last term in eq. (2.22) in section 2. Before doing so, we note here that the first two terms are individually divergent as q → q . However, these singularities cancel between the two terms such that the final result is finite. Some tedious algebra finds the following result for the q → q limit where the superscript numbers indicate derivatives with respect to the arguments, i.e. η, s, t in the case of M. This is a striking result: up to the contaminations from the final term, to which we turn shortly, the correlator gives an estimator of the on-shell N π → N π scattering amplitude. This is contaminated by an off-shell term with superscript (1,0,0) , but with a different time dependence, such that the physical result can be disentangled. This result holds for q = q and Re M µ(q, q ) simply represents the real part of on-shell amplitude for N (q) + π(−q) → N (q ) + π(−q ).
We turn now to C, the contamination resulting from the connected parts of both A k and A * k . For the simplifying case of q = q, the fully connected quantity is given by Applying the change of variables x = (E − E(q))(t − t) then yields (3.14) Equations (3.10), (3.12) (3.13), and (3.14) summarize the main new technical results of this section. In the following subsections we explore their implications for extracting scattering information. In direct analogy to our treatment in section 2, we divide the discussion into threshold and general kinematics.

Threshold kinematics
For E near threshold and q = q = 0, the function H takes the form where the coordinates in the scattering amplitude are given explicitly by Note here that, although the angular dependence is absent for q = q = 0, the amplitude nonetheless depends on the Mandelstam variable t as an artifact of the off-shell kinematics. Next, we use the fact that the Dirac delta function sets p 2 = p 2 = k N π (E) 2 , with the magnitude of back-to-back momenta for non-degenerate particles given by Introducing the shorthand M N π (E) = M 1 η(0, p)| µ(0, p) and evaluating the remaining phase-space integral, we reach Substituting these simplifications into the definition of C, the four-point function reduces to In this result, the top line depends only on M s , the off-shell scattering amplitude projected to zero orbital angular momentum. This differs from M N π because, for terms that are linear in the scattering amplitude, both the incoming and outgoing spatial momenta are set to zero. As a result the Mandelstam t is set identically to zero and no higher partial waves can contribute. Here we have also used the relation between the threshold scattering amplitude and scattering length for the case of non-degenerate particles In direct analog to our analysis of the three-point function, the final step is to expand the final term of eq. (3.21). We reach (3.23) where in the final term we have introduced b = (m N +m π )(t −t) and δm = (m N −m π )/(m N +m π ). We have also introduced In the case of degenerate, non-identical particles the kernel simplifies to (3.28) and in the case of identical particles one requires an additional factor of (1/2) in these expressions.
Expanding K (0) (b, δm) in powers of t −t gives the direct analog of ref. [1] for the four-point function This was already highlighted in the introduction, see eq. (1.10). The relation appears promising as two separate time dependences, together with a N π and a 2 N π coefficients, may give a strong constraint of this threshold observable. As with the three-point function, one reaches a better description by working in the basis of K integrals. See figure 4.

General kinematics
Finally, the general-energy result for the four-point function, with M 0 = E(q), reads where we have used Θ(0, ∆) = 1/2 and also substituted the derivative of the theta function. The final term is built from the integrals together with the coefficients h n , defined via (3.32) Substituting h 0 = Im M µ(q, q ) , demonstrated in appendix A.5, then gives (3.33) The first two lines represent physical information and the final two represent the contaminations that one must remove via fitting.
This completes our presentation of strategies for extracting scattering information from Euclidean correlators, at both threshold and general kinematics. We now turn to discussing how the modified correlation function, defined with a smooth Θ function inserted with the interpolating operators, Here we consider the functions in the case of two degenerate particles with mass m π and with ω q = 3m π and ∆ = m π . might be extracted in practice, in a numerical lattice calculation.

Implementation strategies
In this section we discuss two possibilities to estimate the modified quantities c Θ and c Θ,N π , defined in eqs. (2.4) and (3.1) respectively. We make one adjustment relative to the previous section by replacing the N π → N π four-point function with the simplified version built only out of pion fields. We denote the latter by c Θ,ππ . With straightforward adjustments of the previous section, in particular setting δm = 0 and including the (1/2) for identical particles, this can be used to extract the ππ → ππ scattering amplitude.
In the previous sections, we have not addressed the fact that numerical lattice calculations are necessarily confined to a finite space-time volume. In this section we consider the effects of a finite cubic spatial volume, with periodicity L in each of the three spatial directions, but continue to neglect the thermal effects of a finite temporal extent. Within this set-up, suppose that one has extracted the matrix element from the relevant three-point function. Here the index i labels the values of non-equivalent back-toback momenta. In practice the combined pion state and pion field should be projected to a definite two-pion isospin along the lines of the N π projections involving the C ab coefficients. For the scalar current, J 0 (t) = d 3 x J(x), the total isospin must be 0 or 2 to give a non-zero result.
Similarly, we consider the two-point function and also matrix elements from four-point functions: It is then useful to assemble these quantities into a matrix of correlation functions with rank dictated by the number of pion momenta available in the calculation. Following refs. [7,67], from the solution of the generalized eigenvalue problem on this matrix, one can extract a series of finite-volume energies, E n (L), and matrix elements
Then the matrix element used in Eq. (2.4) to define the modified quantity c Θ q,−q (t 2 |ω q ) can be constructed as provided that the given value of 2ω qi is smaller than E N (L) to avoid systematic errors due to the truncated sum.
The additional effort in measuring the various correlation functions required to construct M (t) is compensated by the fact that we can study several kinematic points for the 2 → 1 transitions, together with the 2 → 2 amplitudes Note that this procedure in fact relies on finite-L, in order to ensure that a finite (and manageable) number of states appear in the regime where Θ 2ω qi − E n (L), ∆ has support. Once the reconstruction is achieved, the resulting c Θ and c Θ,ππ are expected to have residual volume effects, scaling as O(e −L∆ ). Thus, exactly as with the methods described in refs. [59,60], it is important that L∆ is sufficiently large that these can be removed or included in the systematic uncertainty. Future work, perhaps along the lines of refs. [68,69] is needed to fully analyze the finite-L corrections to these correlators as a function of the Θ width as well as the time and energy coordinates.
Alternatively, for situations where an exclusive study is not possible, e.g. for higher center-of-mass energies where a variational approach would be impractical, approximate solutions to the inverse problem represent a viable alternative. The approach is to define coefficients, w Θ (t, ∆, ω qi |t ), satisfying where the kernels w are found numerically, for each value of t, ∆ and ω qi , using Backus-Gilbert-like methods [61,70], or Chebyshev polynomials [62]. These then satisfy While a large ∆ improves the stability of the numerical solutions, the change of basis from e −Et to Θ(E − 2ω qi , ∆)e −Et may also simplify the complexity of the inverse problem (compared to a target Gaussian or Breit-Wigner peak, for example) thereby potentially reducing systematic errors.
Finally, we note that it may also be useful to construct the coefficients to reproduce 1 − Θ(· · · ) instead, i.e. the reflected Θ-function that decreases with increasing E. The resulting estimator can then be subtracted from the original correlator to reach the desired quantity. This seems like a trivial adjustment as the two sets of coefficients should formally be related as In practice, however, the algorithms used to determine the coefficients depend non-trivially on the target functions and w 1−Θ may give a better overall estimate.

Summary and outlook
In this article we have considered extensions of the work of Maiani and Testa [1]. These divide into two main categories: (i) generalizing the results for threshold correlators and (ii) proposing modified correlators that allow one to extend the reach of scattering extractions above threshold.
In the case of threshold kinematics, the key results are presented in sections 2.1 and 3.1. In the former we review the analysis of the ππJ three-point function (with J a scalar current and π a single-particle interpolator) and consider an alternative expansion in kinematic functions that replaces the original expansion of ref. [1], in powers of inverse time (1/t), and appears to exhibit better convergence at small to moderate t values. We also give explicit expressions for the leading off-shellness contaminating the correlator. In section 3.1, we extend these considerations to fourpoint functions, focusing on N π → N π. This gives a tool to extract the N π scattering length, a N π , from a non-standard fit to the Euclidean time dependence, with terms scaling as a N π t and a 2 N π √ t; see also eq. (1.10) in the introduction.
The second class of extensions, targeting scattering information away from threshold, are detailed in sections 2.2 and 3.2, again with a focus on three-and four-point functions, respectively. The key idea is to consider a modified correlator, in which a smooth step function, denoted by Θ(M −M 0 , ∆), is inserted with the operators. Here ∆ denotes the width of the step,M is the hamiltonian boosted to zero-momentum, and M 0 a free parameter. This modification effectively shifts the threshold to higher energies, so that the Maiani-Testa approach can be used to extract the form factor and scattering amplitude for general kinematics. As we stress in the main text, this method is analogous to those based in reconstructing the spectral function [59,60]. The key distinction here is that the target energy range is isolated by the combined effect of the Θ-function and the Euclidean-time translation operator, e −Ĥt . In particular, the latter naturally provides the damping of contributions above the target energy, see also figure 1. This is formalized by the large time expansions presented in sections 2.2 and 3.2. At leading order these give the real and imaginary parts of the physical observables, separated by coefficients with different time dependencies.
In section 4 we briefly describe the implementation strategies for extracting the Θ-modified correlator from numerical lattice calculations. We outline two methods there, one based on the generalized eigenvalue problem (GEVP) and the other on reconstruction methods such as that due to Backus and Gilbert [70]. In both cases, this approach may have an advantage over the other methods because only the low-energy modes of the correlator need to be modified. In particular, the contamination of higher energies is encoded in the large t expansion of the correlation function and therefore profits from knowledge of the functional form (ultimately arising from Lorentz invariance and the form of the time-translation operator and threshold singularities). This information is not directly used, for example, in refs. [59,60]. In section 4 we also briefly discuss the role of finite volume, and the importance of ∆ in suppressing volume effects.
Looking forward, the next steps are to test these methods in a lattice calculation. The most likely initial applications might include extracting scattering lengths from standard three-and four-point functions, especially on lattice ensembles where this information is available from a finite-volume analysis and can provide a direct comparison. This can be used to test, for example, whether one can achieve a more precise determination due to the fact that the observable appears at leading order in the matrix element, rather than as a shift to the energies. This will set the stage for the more ambitious implementation away from threshold, with ultimate target observables including electroweak decay amplitudes, QED corrections to hadronic matrix elements and long-range matrix elements. This will require some extensions of the formal results presented here, but these are expected to be relatively straightforward, along the lines of ref. [60].

Acknowledgments
We warmly thank John Bulava, Marco Cé, Anthony Francis, Dorota Grabowska, Jeremy Green, Taku Izubuchi, Christoph Lehner, Guido Martinelli, Harvey Meyer, Daniel Robaina, Steve Sharpe and Chris Sachrajda for useful discussions, and for previous and ongoing collaborations that helped to inspire this work. MH was supported by UK Research and Innovation Future Leader Fellowship MR/T019956/1. This work is dedicated to Tommaso Sergio Bruno, who provided invaluable motivation in pushing the manuscript to completion.

A.1 Matrix element decomposition
In this section we prove eq. (2.20) of the main text. Begin with a Minkowski-signature representation of A k (q 1 , q 2 ; p) where we have introduced the four-vectors q 1 = (ω q1 , q 1 ), q 2 = (q 0 2 , q 2 ) and P ≡ E(p), P = ω p1 + · · · + ω pn k , p 1 + · · · + p n k . (A.2) Next we decompose the matrix element on the right-hand side into disconnected and connected components to reach This result serves to define the off-shell amplitude M k as the connected component of the matrix element with a simple pole factor removed. The latter is defined using Note that the self-energy contributions from the external leg are absorbed into the off-shell definition.
From the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula then immediately follows where the right-hand side here denotes the 2 → n k scattering amplitude, a function of 3n k − 4 Lorentz invariants arranged in the vector µ.
We conclude by cancelling the common delta function in eq. (A.3) to reach matching eq. (2.20) of the main text.

A.2 General kinematics in the three-point function
In this section we extend our expressions for the three-point function, presented in section 2, to arbitrary kinematics. We begin by rewriting eq. (2.13) as where the indicates that the quantity has been boosted with velocity β = −(q 1 + q 2 )/(ω q1 + ω q2 ).
To be more concrete we introduce the boost matrix Λ(β) which satisfies We then define which in turn allows us to introduce Then eq. (A.7) holds because dΦ k 2ω q2 , (2ω q2 ) −1 A k and B k are each Lorentz scalars.
It remains only to rewrite the t 2 dependence in the exponential. To do so we substitute where we have introduced the four-vector x µ 2 = (t 2 , 0). The final line holds because q 1 + q 2 = 0 and the Delta function inside dΦ k then sets P = 0.
Thus, we conclude that eqs.

(A.19)
Finally we identify the combination of the integrals and the outer-product as an insertion of the identity on the Hilbert space |k, out, p 1 · · · p n k k, out, p 1 · · · p n k | , thereby completing the proof.

A.5 Decomposition of the four-point function
In this section, we demonstrate eq. (3.7) of the main text. We begin by recalling the basic definitions )t e +(E(p)−ωq 1 −ωq 2 )t × Θ s(p) 1/2 − M 0 , ∆ A k (q 1 , q 2 ; p)A * k (q 1 , q 2 ; p) , (A. 33) where ω q1 = m 2 π + q 2 1 , ω q2 = m 2 N + q 2 2 , and A k (q 1 , q 2 ; p) = δ k,1 (2π) 3 2ω q 1 2ω q 2 δ 3 (p 1 − q 1 ) − 2ω q 2 M k η(q , p)| µ(q , p) * η(q , p) − i . (A.34) The product of A k and A * k generates four terms, schematically represented by The doubly disconnected term is discarded and, for the single delta function term, the integral is removed so that one can send → 0 to recover As the final term is inside the integral, the i must be treated carefully. The first step is to break each pole into real and imaginary parts In the analysis below we will consider q = q and only take the limit q → q on combinations of terms for which we know this is safe. Indeed, for many of the individual terms here the limit does not exist.
(A. 51) In addition, an essentially identical analysis gives the result for [2] and [4] c Θ,N π, [2] q ,q + c Θ,N π, [4] q ,q = −e +(E(q )−E(q))t Θ E(q) − M 0 , ∆ 2ω N,q Re M η(q , q)| µ(q , q) η(q , q) , where the only distinction is the exchange of q ↔ q and t ↔ −t . To conclude the tracking of terms we note that c Θ,N π, [6] q ,q = C(t , t|q , q, ω 0 , ∆) , (A.53) as defined in eq. (3.11) of the main text, and in addition that c Θ,N π, [5] q ,q = 0. This final claim holds because each term is defined away from q = q , with the q → q limit to be performed only on the sum. In this prescription, the double Delta function defining c Θ,N π, [5] q ,q cannot be satisfied. This completes the demonstration of eq. (3.7).