Long-range electroweak amplitudes of single hadrons from Euclidean finite-volume correlation functions

Raúl A. Briceño,1, 2, ∗ Zohreh Davoudi,3, 4, † Maxwell T. Hansen,5, ‡ Matthias R. Schindler,6, § and Alessandro Baroni7, 6, ¶ Thomas Jefferson National Accelerator Facility, Newport News, Virginia 23606, USA 2 Department of Physics, Old Dominion University, Norfolk, Virginia 23529, USA Maryland Center for Fundamental Physics and Department of Physics, University of Maryland, College Park, Maryland 20742, USA RIKEN Center for Accelerator-based Sciences, Wako 351-0198, Japan Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland Department of Physics and Astronomy, University of South Carolina, Columbia, South Carolina 29208, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA (Dated: November 12, 2019)


I. INTRODUCTION
Long-range electroweak matrix elements play a central role in modern hadronic physics, from precision tests of the Standard Model (SM) to investigations into the inner structure of strongly interacting particles. In this work, a subset of such matrix elements are considered, defined with a single incoming and a single outgoing hadron, coupled to two local currents that are displaced in time. Key examples of processes for which these types of matrix elements are needed include Compton scattering, double-β decays, radiative corrections to single-β decays, K-K oscillations, and rare meson decays: (i) Deeply virtual Compton scattering, i.e., the conversion from a virtual photon to a real photon via scattering off a charged hadron (γ Ã h → γh), allows one to extract the generalized parton distributions of the target, as proposed by Ji [1]. Such a process can be studied at leading order in quantum electrodynamics (QED) by performing appropriate Fourier transforms of matrix elements involving timedisplaced electromagnetic currents. (ii) A large experimental effort is dedicated to searches for a lepton-number-violating process in nature, namely the neutrinoless double-β decay of certain nuclear isotopes; e.g., see Refs. [2][3][4][5]. If observed, such a decay would establish that neutrinos are Majorana fermions. In order to understand the implication of a potential observation for extensions of the Standard Model, reliable theoretical constraints on the relevant QCD matrix elements are necessary. Moreover, as a supplement to the direct extraction of nucleon matrix elements, nuclear effective field theories [6][7][8][9][10] indicate that the π þ → π − inversion, as well as short distance nuclear effects, may be important contributions to the nn → ppee conversion process. In fact, calculations of the matrix elements relevant for these processes have already started [11][12][13][14], and a class of these computations require bilocal matrix elements [12][13][14][15][16]. (iii) Interest in constraining radiative corrections to nuclear β decays has grown in recent years in light of a discrepancy among theoretical determinations [17][18][19], leading to ∼2-3σ deviation in the Cabibbo-Kobayashi-Maskawa unitary test from V ud . As an indirect input to calibrate several contributions to the β-decay calculations, constraining the time-ordered product of a weak and an electromagnetic current between a neutron and a proton state, for a number of fixed momentum transfers, would be immensely valuable [20]. (iv) In the strange and charmed mesonic sector, bilocal matrix elements can be used to compute the mass splittings between neutral mesons. Another set of observables based on this class of matrix elements arises in rare decays of such mesons, including K → πl þ l − and K → πνν, relevant to searches of CP violation. These processes are currently being measured, for example, in the NA62 experiment at CERN. 1 At present, the most reliable approach to determine such matrix elements from the underlying theory is lattice QCD, a numerical method for statistically estimating QCD correlation functions using Monte Carlo importance sampling. However, since lattice-QCD calculations are necessarily performed in a finite Euclidean spacetime, 2 the relation between the calculated quantity and the physical observable is not always straightforward.
The formalism for extracting long-range matrix elements from finite-volume Euclidean correlation functions was first derived in Refs. [29][30][31], in order to determine the long-range contributions to the K L -K S mass splitting and the CP-violating parameter ϵ K . This framework has enabled determinations of the mass splitting of neutral kaons via lattice QCD [32][33][34][35]. Long-range matrix elements occurring in rare decays of the kaon, such as in K → πl þ l − and K → πνν processes [36][37][38][39], and hadronic double-β decays [12,13,15,16] have also been studied in recent years with lattice QCD.
In this paper, we generalize the previous work to incorporate matrix elements with multiple two-particle intermediate states propagating between the two currents. The result presented holds for kinematics for which any number of two-hadron states can go on shell. We further accommodate particles with intrinsic spin as well as channels with nonidentical particles. In addition, the two local currents in this work are allowed to have a generic structure, including any number of Lorentz indices and nonzero energy and momentum injection. Finally, the expressions presented account for the unphysical mixing of different angular momentum states due to the reduced symmetry of a cubic volume, as well as the physical mixing of different orbital angular momenta in systems with nonzero spin. As a nontrivial check, given these general considerations, the result obtained is proved to satisfy the unitarity of the physical amplitude to all orders. The general form presented is an essential step toward extensions to multihadron bilocal matrix elements.
In order to explore some of the key ideas of Refs. [29][30][31], as well as the general framework of the present study, it is useful to begin with the spectral decomposition of the Euclidean finite-volume matrix element Z d 3 xe −iq·x hM; LjJ E ðτ; xÞJð0ÞjM; Li ¼ L 3 X n e −τðE n −MÞ hM; LjJð0ÞjP n ; LihP n ; LjJð0ÞjM; Li; where Jð0Þ is a generic local current and J E ðτ; xÞ ≡ e Hτ Jð0; xÞe −Hτ defines its Euclidean time translation. In writing Eq. (1), two important assumptions are made: First, we neglect the effects of the finite temporal direction, focusing only on the spatial volume with periodicity L. This is well motivated as lattice-QCD calculations typically work with T > L such that the finite-T effect is a subleading uncertainty. Second, we have introduced a finite-volume single-hadron state, denoted jM; Li, and formally defined by the L-dependence of its eigenvalue, EðLÞ∶ lim L→∞ EðLÞ ¼ M, where M is the physical mass. 3 For the purpose of the introduction only, the state is assumed to have a vanishing spatial momentum. The introductory example is further simplified by taking the two currents to be the same and to be Hermitian. These restrictions will be removed in arriving at the main result of this work. Finally, it should be stressed that the complete set of states inserted on the right-hand side of Eq. (1) is a 1 Further examples can be found, e.g., in a series of recent USQCD whitepapers [20][21][22][23]. 2 For exploratory studies of real-time dynamics in simple theories with classical computations based on Monte Carlo methods, and with quantum computations based on direct implementation of the Hamiltonian time evolution, see Refs. [24][25][26][27][28]. The generalization of these methods to QCD remains formally and practically challenging. 3 For a fixed value of spatial momentum, as well as fixed internal quantum numbers, this defines a unique state. In theories without a mass gap, e.g., QCD þ QED, this separation fails and the distinction is less straightforward.
proper sum, since the finite-volume boundary conditions lead to a discrete spectrum. The sum runs over an infinite tower of states with the same quantum numbers as Jð0ÞjM; Li. Due to the Fourier transform in Eq. (1), the inserted states are understood to be projected to a definite spatial momentum, P n ¼ q.
The corresponding infinite-volume Minkowski observable that we aim to determine is given by where 2ML 3 accounts for the normalization of jM; Li and an extra factor of L 3 results from performing the spatial Fourier transform. 4 Thus, the task at hand is to relate Eqs. (1) and (2). The basic approach for achieving this can be understood as a two-step procedure. First, one must replace the τ-dependent exponentials (e −E n τ ) with poles [1=ðω − E n Þ]. In the case that ω þ M < E, this is achieved by applying the integral operator R ∞ 0 dτe ωτ . However, as stressed in Refs. [29][30][31], many interesting cases arise where the low-lying finite-volume states are below the target ω value, so that the indicated integral is exponentially divergent.
Various approaches have been considered to treat this problem [29][30][31][36][37][38]. In the present work, we are particularly inspired by the techniques that have been applied for the hadronic vacuum polarization (HVP) contribution to the anomalous magnetic moment of the muon [42][43][44][45][46][47][48], as well as the incredible success in extracting excited finite-volume states to determine hadronic scattering amplitudes [49]. Following these ideas, one can envision solving the generalized eigenvalue problem (GEVP) on a matrix of lattice-QCD correlators, in order to reliably determine as many low-lying energies, E n ðLÞ, and matrix elements, jhP n ; LjJð0ÞjM; Lij 2 , as possible. These can then be subtracted from Eq. (1), rendering the R ∞ 0 dτe ωτ integral convergent. In a second step, the subtracted terms are added back in after integration, but with the exponential time dependence replaced by a pole, as detailed in Sec. II. The result of this construction is an intermediate quantity, denoted by T L ðω; qÞ: These manipulations address the effects of the Euclidean signature, but not yet those arising from the finite volume.
The remaining task is then to determine the correction term ΔT L ðω; qÞ ≡ T ðω; qÞ − T L ðω; qÞ, to be added to the sum over finite-volume poles to reach the final physical result. As a special case, it is instructive to consider ΔT L ðω; qÞ when the value of ω þ M is taken below all single-and multihadron thresholds. In this case, the integral over e ðωþM−E n Þτ converges for all E n . It turns out that in this kinematic region, ΔT L ¼ Oðe −mL Þ, where m is the mass of the lightest degree of freedom, i.e., the pion in QCD.
Finally, the formalism presented can be used to reduce systematic uncertainties of long-range matrix elements of hadrons even below two-particle thresholds by identifying and removing contributions from intermediate states that cannot go on shell. A similar approach has been applied to lattice-QCD calculations of the muonic HVP, where the knowledge of the ππ → ππ and γ Ã → ππ amplitudes allows one to estimate the finite-volume effects, as suggested in Refs. [42,43] and implemented in Refs. [48,81]. This paper is organized as follows: Sec. II presents the main results, including the framework that enables the determination of ΔT L . Section III contains the derivation of ΔT L , relying largely on combining existing ideas in the literature to determine scattering and transition amplitudes in the single-and two-hadron sectors. In Sec. IV, a strong check on the formalism is presented showing that the result obtained is consistent with unitarity constraints on the infinite-volume amplitude. In order to demonstrate an application of the formalism, in Sec. V we present an example with a single-channel intermediate state and evaluate numerically all the building blocks for a specific toy example. We conclude in Sec. VI with a summary and outlook. 4 In fact we will use a different but equivalent form for the infinite-volume amplitude in the following sections. We present this version only because it gives additional intuition on the relation between the finite-and infinite-volume amplitudes. For detailed discussion on the form of T given here, and prospects for directly evaluating Eq. (2) numerically, see Refs. [40,41]. 5 For recent reviews on these formal developments and their numerical implementations, we point the readers to Refs. [49,[78][79][80].

II. THE RELATION BETWEEN FINITE-VOLUME EUCLIDEAN AND INFINITE-VOLUME MINKOWSKI AMPLITUDES
The main result of this work can be compactly written as Here, the dependence of the functions on the four-momenta of the initial and final states is left implicit. The left-hand side is the desired infinite-volume Minkowski amplitude, and the right-hand side is a carefully constructed combination of finite-volume Euclidean and Minkowski quantities, which all can be obtained from finite-volume Euclidean correlation functions. In a nutshell, G L is a finite-volume Euclidean matrix element, G <N L is a reconstruction of its low-lying states (using, e.g., GEVP methods), T <N L is a corresponding sum of the low-lying finite-volume poles and ΔT L is the correction term to remove finite-volume effects. The aim of this section is to provide the exact definitions of these functions and the relations among them. As a final general comment, one must note that for kinematics where intermediate twoparticle states may go on shell, the last two terms, T <N L and ΔT L , contain poles that must exactly cancel. For this reason, the two quantities must be treated consistently. This is indicated by the square brackets and the subscript M, H, referring respectively to the relevant 2 → 2 and 1 → 2ð2 → 1Þ scattering amplitudes. This point will be explained in detail toward the end of this section.
We first introduce the basic kinematic notation used throughout. The three-momentum of the incoming hadron state is denoted by P i and that of the outgoing state by P f . The corresponding four-momentum of the initial state is then given by with an analogous relation for the final state given by i → f. Here, M i and M f are the physical particle masses. The infinite-volume hadronic states are denoted by jP i i and hP f j and satisfy the standard relativistic normalization: Next we introduce two local, Minkowski-signature currents J A ðxÞ and O B ðxÞ, with x μ ≡ ðt; xÞ. Here, A and B are collective indices that specify the quantum numbers of the currents. They can specify scalar, axial, vector, tensor, or other types of currents. The underlying Lorentz structure plays no role in the following discussions. The goal is to relate a finite-volume Euclidean correlation function to the infinite-volume amplitude appearing on the left-hand side of Eq. (4), defined as where again the P i and P f dependence on the left-hand side is implicit. Here T denotes the standard time ordering 6 and the subscript "conn." indicates that only the connected contributions to the matrix element are considered. This distinction is only relevant in the forward limit when In a finite cubic volume, with periodicity L in each of the three spatial directions, the L-dependent shifts to the masses, M i and M f , are exponentially suppressed, scaling as e −mL where m is the mass of the lightest low-energy degree of freedom in the theory. Here we assume mL ≫ 1, such that these corrections can be neglected. Thus the fourvectors P i and P f also label the finite-volume states, denoted by jP i ; Li for the incoming hadron and hP f ; Lj for the outgoing. In the convention of this work, finitevolume states are normalized to unity, keeping in mind that in a periodic cubic volume three-momenta satisfy P ¼ 2πn=L Introducing J E A ðτ; xÞ and O E B ðτ; xÞ as Euclidean counterparts of the local currents in Eq. (7), one can define the lattice-QCD correlator most closely related to T as where T E denotes Euclidean time ordering. In the forward limit, the disconnected contribution, must be subtracted from Eq. (10) in order to obtain the purely connected piece of the matrix element as required by 6 Not to be confused with the temporal length of spacetime that was introduced earlier.
Eq. (7). For notational brevity, the hadron momentum labels P i and P f are also omitted from the arguments of the correlation functions.
To understand the issues in extracting T from G, it is instructive to first perform a spectral decomposition of the latter. Defining one obtains The finite-volume states jn; Li and jn; Li in Eqs. (12) and (13), and therefore the c n and c n coefficients in Eq. (14), differ in general not only due to differing three-momenta but also because the internal quantum numbers of the currents may be distinct. For example, processes such as K → πγ are described by setting one current to the electromagnetic current and the other to the weak Hamiltonian, symbolically For such a process, c n;μ receives contributions from states with zero strangeness (such as ππ states), whereas c n;μ contains intermediate states with strangeness equal to −1 (such as Kπ states). Depending on the detailed choices of states, currents, and kinematics in Eq. (14), finite-volume energies may exist for which E n ðL; where ω is the energy carried away by the current J A . As a consequence, intermediate states can go on shell, generating the long-distance parts of these matrix elements. Such states are responsible for the dominant difference between finite-volume Euclidean and infinite-volume Minkowski correlation functions and are the focus of this work. To separate on-and off-shell states, it is useful to introduce cutoff indices, NðωÞ and NðωÞ, such that for n ≥ NðωÞ and n ≥ NðωÞ, one has E n ðL; P i − qÞ þ ω > E i and E n ðL; P f þ qÞ − ω > E f , respectively; i.e., the finitevolume intermediate states are off shell up to a current energy of ω. Taking NðωÞ and NðωÞ larger than the minimum requirements poses no problem, and as explained in more detail below, will likely be advantageous in practical implementations.
With this discussion, one can define where to keep the notation simple, the superscripts < N and ≥ N are taken to be the representative of both the N and N dependence of the functions, and the ω dependence of N and N is left implicit. As mentioned in the Introduction, the c n (c n ) coefficients for the intermediate states 0 to N − 1 (N − 1) can be separately evaluated in a dedicated lattice-QCD calculation of three-point functions formed with optimized operators; see Eqs. (12) and (13). In Eq. (16), the subtracted integral is convergent by construction, as NðωÞ and NðωÞ are chosen such that e ωτ ½G L ðτ; qÞ − G <N L ðτ; qÞ decays with increasing jτj. The result of the integration carries no memory of the Euclidean signature and thus brings us closer to the stated goal of recovering T . However, the approach is clearly incomplete, since the intermediate states labeled from 0 to N − 1 (N − 1) are not accounted. Additionally, the finite-L effects have yet to be addressed.
One can now define which can naively be constructed using exactly the coefficients and energies defining G <N L ðτ; qÞ. As explained below, in practice one must take a slightly different approach to properly treat finite-volume effects. Note that T ≥N L ðω; qÞ, defined in Eq. (15), satisfies the same decomposition as that given here for T <N L ðω; qÞ, but with the sums running from N, N to ∞. Thus, the combination gives the sum over all finite-volume poles from 0 to ∞ and differs from the target quantity, T ðω; qÞ, only by the finitevolume effects encoded in ΔT L ðω; qÞ.
As stressed in the beginning of this section, the twoparticle poles within T <N L must be exactly canceled by those in ΔT L . This is essential because (i) these poles are artifacts of the particular volume, L, and cannot be part of the physical quantity T , and (ii) the formalism holds for values of ω arbitrarily close to, indeed exactly coinciding with, the poles in a given lattice volume. Note that the same is not true for the poles within T ≥N L ðω; qÞ, which are safely above the range of allowed ω values. The requirement of exact cancellation of divergences in T <N L and ΔT L means that these two quantities must be treated in a consistent way. To explain this properly, we first need to give some details on the construction of ΔT L ðω; qÞ.
Both the generalized Lüscher and Lellouch-Lüscher formalisms are by now well understood, and here only the key equations are provided for completeness. First, the Lüscher formalism provides a relationship between the finite-volume spectrum in the two-particle regime and the infinite-volume amplitude. This can be written as [58] det ½F ½x ðL; where the superscript ½x is set equal to either J jP i i or OjP i i and serves as an indicator of the internal quantum numbers for the two-particle states. Here, M ½x is the infinite-volume 2 → 2 scattering amplitude and F ½x is a known kinematic function, given explicitly in Eq. (43) of Sec. III. As explained in detail in that section, the matrix and determinant space used here is a Kroneckerproduct space of ðorbital angular momentumÞ ⊗ ðspinÞ ⊗ ðflavor channelsÞ. Depending on the time ordering of the currents, the quantization condition is to be evaluated at qÞ is the four-momentum of the current J A . One must then identify the set of solutions, in P 0 for which the left-hand side vanishes. This gives the finite-volume spectra, denoted by E n ðL; P i − qÞ (for ½x ¼ ½J jP i i) and E n ðL; P f þ qÞ (for ½x ¼ ½OjP i i).
From the same building blocks that define the quantization condition, one can construct a matrix, F , that plays an important role throughout this work: where ½x and P ≡ ðE; PÞ are taken as generic representatives of the two-particle quantum numbers and momenta. This matrix has poles, E n , whenever Eq. (19) These factors allow one to relate finite-and infinitevolume matrix elements. For the present setup the relevant relations are where various infinite-volume 1 þ J → 2 transition amplitudes are introduced, e.g., The amplitudes defined with an "out" state are understood as column vectors and those with an "in" state as row vectors, acting on the space of angular momentum, spin, and flavor of the two-hadron state. Now given a set of finite-volume energies, E n ðL; P i − qÞ and E n ðL; P f þ qÞ, and the various finite-volume matrix elements noted above, one can constrain all quantities needed to remove the finite-volume effects from T L ðω; qÞ, to extract T ðω; qÞ. The correction is given by where each term on the right-hand side has the structure ½row vector × ½matrix × ½column vector, defined in the space of two-particle degrees of freedom. In practice, the determination of the scattering amplitudes, M, as well as the transition amplitudes, H, will generally require fits of all finite-volume data to a set of parametrizations. Given a particular fit, one can then recover the finite-volume matrix elements, and thus the coefficients c n and c n , as well as the corresponding finitevolume energies. It follows that the sum over low-lying poles, T <N L ðω; qÞ, can be evaluated in two ways, either directly from the finite-volume data or from the energies and matrix elements predicted through a global fit. Given the statistical nature of the lattice-QCD data, for this analysis it is crucial that the second approach is used. Only in this way does one assure an exact cancellation of the poles in ΔT L and T <N L is reached. This requirement of a common parametrization for these two pieces is emphasized in Eq. (4) by the M; H subscript. The cancellation of poles will also be illustrated in Sec. V for a particular example.
This completes the discussion of the procedure involved in extracting the infinite-volume amplitude T using the inputs that are defined in a finite Euclidean spacetime, as summarized in Eq. (4). The quantities that must be evaluated in implementing this procedure, as well the relationships among them, are depicted in Fig. 1 for further clarity.

III. FINITE-VOLUME CORRECTIONS TO THE FOUR-POINT CORRELATION FUNCTION
In this section, the expression for the additive finitevolume function, ΔT L ðω; qÞ, given in Eq. (30), is derived. The derivation requires a diagrammatic representation of finite-volume correlation functions, first laid out for systems with identical scalar particles by Lüscher in Ref. [50], and Lellouch and Lüscher in Ref. [72], and presented in a purely quantum-field-theoretic context by Kim, Sachrajda, and Sharpe (KSS) in Ref. [52]. The approach of KSS is also the starting point of Ref. [31] in the analysis of long-range effects in K-K mixing. These techniques have since been generalized to accommodate any number of open twoparticle channels with arbitrary masses and spin [43,[53][54][55][56][57][58][73][74][75]. Here we adopt this formalism to extend the work of Ref. [31].
Consider the infinite-volume Minkowski-signature correlation function where one of the four fields is left at the origin of position space to avoid the overall momentum-conserving delta function. The Lehmann-Symanzik-Zimmermann (LSZ) reduction formula implies that this correlation function Eq. (15) Eq. (16) Eqs. (17), (30) Eq. (4) has poles associated with the incoming and outgoing single-particle states, with quantum numbers set by Ψ and Ψ 0 . By amputating the external legs and placing the momenta P i and P f on shell, one arrives at an alternative expression for iT , first defined in Eq. (7), where the fields are normalized such that hP i jΨð0Þj0i ¼ h0jΨ 0 ð0ÞjP f i ¼ 1, and as before the P i and P f dependence of the amplitude here and in Eq. (34) is left implicit. The finite-volume counterpart of the correlation function in Eq. (31) can be written as with the momentum and position coordinates still carrying Minkowski signature. The integrals over time coordinates are as above, but those over space coordinates are restricted to the finite cubic volume, as indicated by the L subscript. A similar LSZ reduction can be implemented here to reach This provides an alternative definition for T L ðω; qÞ, the finite-volume counterpart of the Compton amplitude introduced in Eq. (18). At this point, the goal is to quantify the difference between Eqs. (32) and (35). These quantities correspond directly to the correlation functions studied by KSS [52] and generalized by two of us in Ref. [75]. For completeness of the presentation, some details of the derivation will be repeated here. The first step is to derive a diagrammatic representation of the infinite-volume correlation function that explicitly displays all two-particle intermediate states.
The result, illustrated in Fig. 2, can be represented algebraically in the following compact form: where we have chosen to work directly with T instead of C and have thus amputated all single-particle operatordependent terms on both sides of the equations.  (29), and then discarding all diagrams that fall into disconnected pieces when any two lines, carrying the total four-momentum of the system, are cut. The quantity M, referred to as the Bethe-Salpeter kernel, is defined in exactly the same way, but with the 2 → 2 hadronic amplitude in place of H ½X 2→1ð1→2Þ . 7 The quantum numbers associated with M and I (as well as those of S introduced below) are deduced from those of H ½X 2→1ð1→2Þ present in each term.
In Eq. (36), the various two-particle irreducible (2PI) quantities are combined via integrals over two-particle loops, denoted by ⊗ I ⊗. This symbol is thus built from two fully dressed propagators, the requisite symmetry factor and the loop integral. The ⊗ symbol is meant to stress that the quantities are connected in a complicated way, in particular by integration over the loop momentum. This results in H 2→1 , M, and H 1→2 being evaluated at offshell values of the momenta as well. As an explicit example, one has where Δ a1 and Δ a2 are the fully dressed propagators for particles 1 and 2 in channel a, in accord with the quantum 7 This quantity was named B in Ref. [82]. numbers of state OjP i i. ξ a ¼ 1=2 if the particles in channel a are identical and 1 otherwise. Momentum dependences on the left and right of the semicolon in the arguments refer to the momenta of incoming and outgoing pairs of hadrons, respectively. From the definitions of ⊗ I ⊗ and the various 2PI quantities, three important identities directly emerge: with analogous relations for J ↔ O in Eqs. (38) and (39). These are used in the derivation to express the finitevolume function, T L , in terms of the physical amplitudes appearing on the left-hand sides. Finally, T in Eq. (36) is defined by all diagrams not included in the other terms. This includes all connected 2PI contributions to T as well as diagrams that contain a single-particle intermediate state, as shown in Fig. 2.
The purpose of the skeleton expansion for T is twofold. First, it allows one to identify all sources of singularities, as well as imaginary contributions to T . This will play a key role in Sec. IV where the unitarity of the amplitudes is established. Second, it enables identifying all powerlike L-dependences in its finite-volume analog, T L , defined in Eq. (34). In particular, as discussed in more detail in Refs. [52,[55][56][57][58]73,75], one can show that Here all quantities are as above except for ⊗ S ⊗. This is defined exactly as ⊗ I ⊗ but with the replacement R d 3 k=ð2πÞ 3 → ð1=L 3 Þ P k ; see the example in Eq. (37). The sum runs over all three-momenta allowed by the finitevolume boundary conditions, in particular k ¼ 2πn=L where n is a vector of integers.
Having set up the diagrammatic expansions for T L and T , one can now substitute ⊗ S ⊗¼⊗ I ⊗ þF, where F is a finite-volume cut, defined by this relation and given more explicitly in Eq. (43) below. In words, the substitution rule reads: replace each sum with an integral plus a sum-integral difference, encoded by the symbol F. Two key simplifications then arise. First, the off-shell contributions from H 1→2 , M, and H 2→1 lead to only exponentially suppressed volume effects that can be neglected. As a result, all functions adjacent to F are projected to their on-shell momenta. Second, the propagators Δ a1 and Δ a2 are expanded about their on-shell point meaning that only the physical mass enters F. In particular, for spinning particles, the spin structure of the propagator becomes simple helicity projectors acting on the neighboring functions.
Non-singular below three-particle thresholds FIG. 2. A diagrammatic expansion of the Compton amplitude defined in Eqs. (31) and (32). Here, the aim is to display explicitly all two-hadron intermediate states in the s-channel, as represented in Eq. (36). In addition, single-particle intermediate states are separated, though these play a less important role in the finite-volume analysis of this section. The two-current vertex, as well as the 1 → 2 kernel, contain no two-particle singularities below three-particle production thresholds. Representative contributions to the former and the latter are shown for a given theory in the bottom-right and bottom-left panels, respectively. The black circle represents the 2 → 2 scattering amplitude, M.
Substituting ⊗ S ⊗¼⊗ I ⊗ þF in Eq. (41) and summing all ⊗ I ⊗ dependences, and upon restoring the superscripts on all functions, one deduces This is the main result of the derivation, analogous to Eqs. (44) and (84) of Refs. [52,75], respectively. Summing the geometric series in n then directly gives Eq. (30) of the previous section. Figure 3 provides a diagrammatic summary of the procedure that isolates the finite-volume effects in the four-point function to all orders. At this point it remains only to define F ½XjP i i , the finite-volume cut matrix defined on the Kroneckerproduct space of ðorbital angular momentumÞ ⊗ ðspinÞ ⊗ ðflavor channelsÞ. 8 In order to provide an explicit form, one can change to the basis of total angular momentum and use the compact notation introduced in Ref. [75]. Let Jm J denote the total angular momentum and its azimuthal component and S and l the total spin and orbital angular momentum. Also, label the channel space using indices a and a 0 . Then the F ½XjP i i function can be written as [58,75] F ½XjP i i fJg;fJ 0 g ðP; LÞ ≡ ξ a δ aa 0 δ SS 0 where ξ a is the symmetry factor defined after Eq. (37). Note that the superscript ½XjP i i on F is related to the quantum numbers of allowed two-hadron channels in this relation. The particles in channel a are labeled with the numbers 1 and 2, and their corresponding masses are m a1 and m a2 . Then the magnitude of back-to-back three-momentum in the center-of-mass frame is given by q Ã a , which is the solution to In Eq. (43), we have also introduced ω a1 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , with the latter defined as the spatial part of the four-vector reached by boosting ðω a1 ; kÞ with boost velocity β ¼ −P=E. E is the total energy of two hadrons in the lab frame, i.e., E ¼ γE Ã . 9 Equation (43) provides the final missing piece for the main result of this work, Eq. (4).
Second, to obtain the generalized Lellouch-Lüscher relation, one must match the residues at the poles defining T L , determined in two different ways. To this end, consider a special case of Eq. (35) for which J ¼ O, the initial and final states are the same, and P i ¼ P f . In this case, one can first aim for an expression for T L in terms of finite-volume matrix elements. This follows directly from Eqs. (12), (13), (17), and (18). Choosing the kinematics so that E i − ω is arbitrarily close to a finite-volume energy, E n , one obtains To reproduce a second expression for T L near the finite-volume pole, one notes that the poles are encoded in Eq. (30) and thus Equating these two expressions reproduces Eqs. (22) Before closing this section, we elaborate further on the scope of the results derived. In particular, we note that the general angular momentum, spin and flavor structure of the finite-volume function, F, defined in Eq. (43), is a new feature that substantially complicates the extraction of Compton-like amplitudes using lattice QCD. However, the main challenges already arise in the determination of coupled hadronic scattering amplitudes (via the generalized Lüscher formalism) as well as one-to-two transition amplitudes (using the generalization of the Lellouch-Lüscher relation). As a result of significant numerical and algorithmic progress in recent years, implementation of these ideas is already well underway.
The central complication, whenever multiple flavor or angular-momentum channels need to be considered, is that one no longer has a direct one-to-one mapping between finite-and infinite-volume observables. In such cases, the only known approach is to determine a large set of finitevolume energies and matrix elements, for multiple choices of box size and total momentum, and to subsequently perform a global amplitude analysis using the relevant finite-volume formalism. This was first implemented for hadronic amplitudes with partial waves unphysically mixing in the finite volume [83], and was later extended to cases where multiple flavor channels are kinematically allowed and must be disentangled [59,60,65,66]. Most recently, the methods have been applied to particles with intrinsic spin, inducing also an infinite-volume coupling between different orbital angular momenta [67,68] and to transition amplitudes in which states are coupled via the spin and momentum of the external current [84,85]. This progress provides strong empirical evidence that there is no significant obstacle for implementing the proposed formalism for Compton-like amplitudes in the kinematic region considered in this work.

IV. UNITARITY CHECK
This section provides an additional check on Eq. (4) by proving that the result is consistent with the S-matrix unitarity. As was recently shown in the context of threeparticle scattering in Ref. [86], all-orders perturbation theory can be used to directly extract all imaginary contributions to the scattering amplitude in a given kinematic region. Since the diagrammatic description emerges from a unitary theory, the result of summing over all contributions must automatically respect any consequence following from the unitarity relation, i.e., S † S ¼ 1, in both the finite and the infinite volumes. Nevertheless, it is instructive to see how the final result directly satisfies the expected constraint.
Let us first work out the consequences of unitary for T . To this end, consider an all-orders expression for the Compton amplitude in the infinite volume, as depicted in Fig. 4. For the purpose of this check, the two local FIG. 3. Shown is the full four-point correlation function in a finite volume. In the second line only the s-channel diagrams are shown explicitly, and the u-channel counterparts follow by switching the current vertices. The notation is similar to that used in Fig. 2. The vertical dashed lines denote a contribution from the finite-volume F function, defined in Eq. (43). The open circles denote the overlap of a single-hadron state with the vacuum.
currents are assumed to be the same and to be Hermitian , and the initial and final states are set to coincide with P i ¼ P f . 10 When necessary, the dependence on the kinematic variable ffiffi ffi s p ¼ E Ã is made explicit throughout this section, while the full dependence of the functions on kinematic variables, including on the momentum transfer of the currents, will be suppressed for brevity.
In the kinematic region in which only two-particle states can go on shell, one can analyze all diagrams by systematically isolating their imaginary contributions. These are present only in the two-particle s-channel loops and are proportional to the phase-space factor ρ, defined as The phase-space factor, ρ a ðsÞ, arises from the imaginary part of the iϵ pole prescription in the two-particle loops. To capture all contributions to the amplitude, it is convenient to split the full iϵ integral into its real and imaginary parts. Performing this for a purely hadronic 2 → 2 scattering amplitude yields the standard relation where K is the K matrix, defined to have the same diagrammatic expansion as the scattering amplitude, but evaluated with the principal-value prescription for the s-channel loops. Similarly, the principal-valued version of H 1→2 can be denoted by H 1→2 , which is related to H 1→2 via A similar relation for H 2→1 can be realized in terms of H 2→1 , As with the various finite-volume quantities in the previous section, Eqs. meaning that the amplitudes have the same complex phase in each partial wave. This is Watson's theorem, and by returning to the general matrix space, one recovers a simple generalization for the case of multiple two-particle channels [73].
In direct analogy to these more standard results, one can further determine an expression for the Compton amplitude that is consistent with unitarity to all orders in perturbation theory. This, in the kinematic region below three-hadron production thresholds, divides into a purely real principalvalued version of T , denoted T, as well as a series of ρ cuts, analogous to those appearing in Eqs. (49)-(51), above. In addition, one must include a single-particle piece represented by the second and third terms on the right-hand side of Fig. 4. To express this, let us denote by fðq 2 Þ the 1 → 1 matrix element with a single current insertion and label the corresponding particle mass by m. Combining all terms gives FIG. 4. Shown is the representation of the Compton-like amplitude to all orders in perturbation theory, in which the contributions from imaginary parts (proportional to the phase-space factor ρ) are isolated. The black dots are either the full 1 → 1 vertex with a single current insertion, if, or the 1 → 2 or 2 → 1 scattering amplitudes, iH. The gray circles define iH and iK, respectively, which are the counterparts of transition and scattering amplitudes when the principal-valued prescription has been used in the s-channel loops. Further detail on these quantities is given in the text.
where s ¼ ðP − qÞ 2 and u ¼ ðP þ qÞ 2 are the usual Mandelstam variables, and the last term indicates the inclusion of terms in which the incoming and outgoing currents are interchanged, as shown in Fig. 4. At this point, it is apparent that iT shares certain features with the previously considered amplitudes, M and H 1→2ð2→1Þ . First, besides the single-particle pole, the only singularity in the amplitude arises from the phase-space factor ρ a ðsÞ ∼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s − ðm a1 þ m a2 Þ 2 p . This threshold singularity is common to all three amplitudes. Second, in the case that a single two-particle channel is open, the third term in Eq. (53) breaks into a sum over the angular momentum modes and the phase of each term matches that of M ðlÞ and H ðlÞ .
Below the lowest-lying two-particle threshold, iρ and thus also M and H 1→2ð2→1Þ become purely real. Similarly, the only possible imaginary contribution to the Compton amplitude below threshold is due to the iϵ prescription within the single-particle pole. Above the lowest lying twoparticle threshold, iρ, M, H 1→2ð2→1Þ , and T are all complex valued. In particular, from Eqs. (49)-(51) and (53) one can readily show that where ρ a ðsÞ ≡ ρ a ðsÞΘðs − ðm a1 þ m a2 Þ 2 Þ, with Θ being the Heaviside step function.
Having established the unitarity constraint on the various infinite-volume amplitudes entering the formalism, it can now be shown that Eq. (4) is consistent with these expressions. To this end, it is sufficient to demonstrate ImΔT L ¼ ImT , which, together with the fact that G ≥N and T <N are real, demonstrates the desired consistency. One first notes that where all s dependences are dropped for brevity and we have used the definition of F given in Eq. (20). Here, we use the definition ImX ≡ ðX − X † Þ=ð2iÞ for a generic matrix, X, together with the fact that H 2→1 M −1 and M −1 H 1→2 are real valued and equal, entry by entry. The next step is to make use of the fact that M −1 þ F, and thus also its inverse, are Hermitian matrices. This leads to the substitution where, in the second line, we have used  60) and (61) into Eq. (58) leads to as desired.
In summary, we have provided an expression for the Compton amplitude, in terms of purely real building blocks, that is exactly consistent with unitarity, and we have used this to express ImT in terms of transition amplitudes and phase-space factors. This expression is shown to be consistent with Eq. (4), the main result of this work. The decomposition of T is built from the singleparticle form factor and propagator, the phase-space factor ρ, the amplitudes H 1→2 , H 2→1 , and M, as well as a realvalued short distance piece, denoted by T. In fact, since all but the last function can be obtained from two-and threepoint correlation functions in a finite-volume study, it is natural to think of T as the target observable, the determination of which requires the method summarized by Eq. (4). In the next section we discuss this point in further detail and highlight one subtlety that arises when considering amplitudes for which the corresponding K matrix has real-valued poles.

V. NUMERICAL IMPLEMENTATION
In this section, we provide a numerical example to illustrate the approach outlined in Sec. II above and summarized by Eq. (4). To reduce technical complications, the example involves kinematics for which only a single channel is open. The channel is composed of two identical LONG-RANGE ELECTROWEAK AMPLITUDES OF SINGLE … PHYS. REV. D 101, 014509 (2020) scalars, each with mass m. It is further assumed that only the lowest partial wave (l ¼ 0) contributes. The initial and final single-hadron states are taken to be the same and are put to rest. Finally, the external current is taken to be a scalar, and as a result, the form factor fðQ 2 Þ and transition amplitude Hðs; Q 2 Þ are scalar functions. The subscripts 1 → 2 and 2 → 1 on H will be dropped in this section. The dependence on kinematic variables is displayed using a Lorentz-invariant notation, where s ¼ ðP − qÞ 2 is the invariant mass of the system and q 2 ¼ −Q 2 , where q is the momentum transferred by the first current. As stressed in the previous section, the formalism can only be applied to amplitudes that exactly satisfy unitarity. In the case of the hadronic scattering amplitude, MðsÞ, this means one must use Eq. (49) with a real-valued KðsÞ. With the standard relation to the scattering phase shift δðq Ã Þ, and the use of a Breit-Wigner parametrization for the latter, K is fully constrained. Here, m R is related to the mass of the resonance and g characterizes the coupling to two-particle states. In the numerical results, the following values are assumed: m R ¼ 2.5m and g ¼ 3.0. As can be seen in Fig. 5, these values lead to a standard peaklike structure near s ¼ m 2 R and a visible cusp at threshold, s ¼ 4m 2 .
The truncation of MðsÞ to a single partial wave reduces the two-hadron energy quantization condition, Eq. (19), to a simple algebraic relation, that can readily be solved numerically. The corresponding finite-volume spectrum is shown in the upper panels of Fig. 6 for two choices of the total three-momentum. From the unitarity constraint on the 1 → 2ð2 → 1Þ transition amplitude, KðsÞ must be combined with the real-valued quantity, Hðs; Q 2 Þ, as in Eq. (51), to obtain H. An all-orders perturbation-theory argument requires that for any pole singularity of KðsÞ, there must be an associated pole singularity in Hðs; Q 2 Þ [87]. This motivates rewriting Hðs; Q 2 Þ as [84,85] Hðs; where Bðs; Q 2 Þ is a smooth real function. Here we take Bðs; Q 2 Þ ¼ fðQ 2 Þ=3, where fðQ 2 Þ is the single-particle form factor, which can be parametrized by a simple monopole form, The resulting functional form of Hðs; Q 2 Þ for a range of kinematics is shown in the second panel of Fig. 5.
For the remainder of the discussion, we consider kinematics for which only the s-channel contributions to T are physical. As a result, only these can lead to powerlaw finite-volume effects, and the u-channel contributions can be ignored. Given the specified forms of MðsÞ and Hðs; Q 2 Þ, one can now determine ΔT L using Eq. (30). The result is plotted in the middle panels of Fig. 6 for a volume with mL ¼ 8, and two values of the total spatial momentum of the current. The initial and final states are fixed to be at rest. Note that the real part of ΔT L exhibits a series of poles that coincide with the energy levels shown in Fig. 6. The imaginary piece has a cusp at threshold and has the same qualitative peak structure as MðsÞ and Hðs; Q 2 Þ. This is consistent with the fact that ImΔT L satisfies Eq. (62). Furthermore, given Eqs. (12) and (22), T <N L can be constrained within the parametrization presented. The lower panels of Fig. 6 plot T <N L þ ΔT L for the two spatial momenta considered, exhibiting the exact cancellation of the poles that was emphasized at the end of Sec. II.
Next, we numerically demonstrate that the imaginary part of ΔT L is equal to that of T . To proceed, T in Eq. (53) needs to be constructed in terms of smooth functions, and for this purpose, the issue of K-matrix poles, mentioned briefly at the end of Sec. IV, must be addressed. Note that when Eq. (66) is applied, the third term in brackets in Eq. (53) develops real-valued poles corresponding to those of the K matrix, In the full Compton amplitude, this must be exactly canceled by a pole in T, motivating one final reparametrization where Sðs; Q 2 Þ is a smooth function in the region of interest. We conclude that, in this kinematic region, the analytic structure of the Compton amplitude is given in terms of known functions, plus an additional smooth contribution, Sðs; Q 2 Þ, which will be set to vanish for this toy example. This completes our construction of a parametrization for the Compton amplitude, with the result plotted in Fig. 7. This numerically confirms the result derived in Sec. IV and summarized in Eq. (62) that ImΔT L ¼ ImT .
To complete the illustration of the proposed formalism, in Fig. 8, we plot the subtracted correlation function, e ωτ ½G L ðτÞ − G <N L ðτÞ, for a set of kinematics and for various subtractions, as a function of Euclidean time. This defines the integrand in Eq. (16) leading to T ≥N L . While this quantity can be obtained directly from Eq. (4) given the constrained values of T and T <N L þ ΔT L (by performing an inverse integral transform), one can also approximate the function by summing up contributions from a sufficiently large number of lowest intermediate finite-volume states; see Eq. (14). Nine such states are used to obtain the results of Fig. 8. The function exhibits the expected behavior, with an increased number of subtractions accelerating the falloff at large Euclidean separations. It must be stressed that in the case where two-particle states can go on shell, a number of subtractions are required to render the integral convergent. However, in practice it may be profitable to also subtract states in the off-shell region to accelerate the convergence and improve the overall uncertainty on the resulting amplitude. This is demonstrated in FIG. 7. Real (red line) and imaginary (blue line) parts of the Compton amplitude. The initial and final states are fixed to be at rest. The darker lines correspond to currents that have zero momentum, and the faded lines correspond to q ¼ 2π½001=L, where mL ¼ 8 to match the functions illustrated in Fig. 6. the first panel of Fig. 8, where the chosen value of ω does not provide sufficient energy to produce any on-shell intermediate states in this example, but the states close to threshold cause only slowly damping functions of τ if not subtracted.
As a final remark, we emphasize that this numerical example only serves as a quantitative demonstration of the various features of the building blocks of the master equation (4) in a toy model, which, however, was built consistent with physically motivated scenarios. While the approach in this section was to construct the final-volume quantities from the knowledge of infinite-volume amplitudes, in the realistic use of the formalism presented, the known quantities will be the finite-volume two-, three-and four-point Euclidean correlation functions, and the infinitevolume Minkowski observable will subsequently be deduced. In particular, one must note that in the parametrization of this section, the smooth function Sðs; Q 2 Þ is the only unknown infinite-volume quantity whose value cannot be fixed by the knowledge of two-and three-point functions in the theory, and only a determination of the four-point correlation function as defined in Eq. (10) can constrain its value, following the procedure of this work. 11

VI. CONCLUSION AND OUTLOOK
The formalism presented in this work offers a path from Euclidean finite-volume correlation functions of timedisplaced local electroweak currents to long-range contributions to hadronic amplitudes in an infinite Minkowski spacetime. Given the complicated nature of the desired amplitudes, it should not be too surprising that the relation derived in this paper requires a detailed understanding of various building blocks on both ends of the mapping. Figure 9 summarizes the conclusions of this work and provides guidance on how lattice-QCD quantities may be used to access the physical infinitevolume amplitudes of interest.
The main result of this work is summarized in Eq. (4). In short, an extract of the original correlation function is identified that is independent of the time signature in the theory, denoted in Fig. 9 as the subtracted one-body matrix element. This is achieved by removing contributions to the finite-volume four-point correlation function that arise from the lowest-lying intermediate states, including all the states that can go on shell. A closed form for the necessary additive piece, ½T <N þ ΔT M;H in Eq. (4), is then provided. Not only does this form remove all power-law finitevolume effects below three-hadron productions thresholds but also it restores the correct analytic structure of the infinite-volume Minkowski amplitude. This additive piece can be evaluated separately from dedicated lattice-QCD studies of two-and three-point functions.
The approach of this work generalizes the formalism of Refs. [29][30][31]. In particular, the results presented hold not only for single-hadron long-range electroweak transitions but also for transitions involving the vacuum in the initial/ final states, such as for the matrix element incurred in studying the QCD structure of the photon [88]. Arbitrary quantum numbers, such as spin, flavor, partial waves, and total momentum, are incorporated in the formalism, and the possibility of multiple coupled partial-wave or flavor FIG. 8. Shown is the integrand of Eq. (15) for the example of this section for a range of values of N and ω, where initial and final states are fixed to be at rest. The color coding in the upper two panels corresponds to that of the lower panel. 11 A similar consideration is presented in Refs. [12,13] in the context of matching lattice-QCD results to an effective field theory description of a nuclear double-β decay. channels in intermediate two-hadron states is accounted for. An explicit map of the workflow for future numerical implementation of the formalism is shown in Fig. 1, with a reference to quantities that are defined in various equations throughout this paper.
The general framework of this work can further serve as guidance on how to address more complex scenarios, such as considering two-hadron initial and final states (relevant for neutrinoful and neutrinoless double-β decay), or extending the kinematic reach of the problem so that more than two hadrons can be produced on shell in the intermediate states. Extension to kinematic regions beyond the three-hadron productions may be possible given the recently developed technologies for the determination of three-hadron observables from a finite-volume study [82,87,[89][90][91][92][93][94][95][96][97][98][99][100][101]. More immediately, one can imagine extending these ideas for processes such as γ Ã γ Ã → ππ, which would be relevant for dispersive analyses of the hadronic light-by-light contribution to muon g − 2 [102].