Form factors of two-hadron states from a covariant finite-volume formalism
Alessandro Baroni Raúl A. Briceño Maxwell T. Hansen Felipe G. Ortega-Gama


JLAB-THY-18-2878 
CERN-TH-2018-263 






[]abaro008@odu.edu


Department of Physics and Astronomy University of South Carolina, 
 
712 Main Street, Columbia, South Carolina  29208, USA






[]rbriceno@jlab.org

Thomas Jefferson National Accelerator Facility, 12000 Jefferson Avenue, Newport News, Virginia 23606, USA
 Department of Physics, Old Dominion University, Norfolk, Virginia 23529, USA




[]maxwell.hansen@cern.ch
Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland





[]felortga@jlab.org

Thomas Jefferson National Accelerator Facility, 12000 Jefferson Avenue, Newport News, Virginia 23606, USA
Department of Physics, College of William and Mary, Williamsburg, Virginia 23187, USA











In this work we develop a Lorentz-covariant version of the previously derived formalism for relating finite-volume matrix elements to 
2 + 𝒥→2 transition amplitudes. We also give various details relevant for the implementation of this formalism in a realistic numerical lattice QCD calculation. Particular focus is given to the role of single-particle form factors

in disentangling finite-volume effects from the triangle diagram that arise when 𝒥 couples to one of the two hadrons. This also leads to a new finite-volume function, denoted G, the numerical evaluation of which is described in detail. As an example we discuss the determination of the ππ + 𝒥→ππ amplitude in the ρ channel, for which the single-pion form factor, F_π(Q^2), as well as the scattering phase, δ_ππ, are required to remove all power-law finite-volume effects. The formalism presented here holds for local currents with arbitrary Lorentz structure, and we give specific examples of insertions with up to two Lorentz indices.

 




Form factors of two-hadron states from a covariant finite-volume formalism

    Felipe G. Ortega-Gama
    December 30, 2023
============================================================================
























§ INTRODUCTION



In recent years, interest in hadron spectroscopy has increased significantly, primarily due to various experimental discoveries of unconventional excitations.[For recent reviews of the experimental and theoretical understanding of these states we point the reader to Refs. <cit.>.] This has led to an abundance of theoretical proposals as to the underlying nature of the unexpected states. Possible explanations range from multi-hadron molecules to compact multi-quark configurations, to kinematic singularities arising from specific Feynman-diagram topologies <cit.>. In many cases, experimental data alone is not sufficient to distinguish between available explanations, and thus many questions remain unresolved. 

In some cases, theoretical calculations can provide access to experimentally-unavailable quantities that may shed light onto the structure of the QCD spectrum. With this goal in mind, in this work we present a detailed framework that will allow for the rigorous lattice-QCD calculation of transition amplitudes, mediated by electroweak or other external currents, involving two hadrons each in the initial and final states. We abbreviate our process of interest by 2 + 𝒥→2, where each bold-faced 2 counts the QCD-stable hadrons in the state and 𝒥 represents a generic, local external current.

The approach discussed here is based on prior formalism developed by two of us in Ref. <cit.>.[This, in turn, was inspired and guided by the work of Refs. <cit.>.] In Sec. <ref> we present a slightly modified version of this formalism in which all infinite-volume quantities are Lorentz covariant and the single-particle matrix elements that enter, abbreviated 1 + 𝒥→1, are expressed in terms of standard Lorentz-invariant form factors. After extracting the 2 + 𝒥→2 transition amplitudes, one can proceed to determine form factors as well as distribution functions[Distribution functions are accessed in lattice calculations via spatially non-local operator insertions <cit.>. These may suffer from further finite-volume effects associated with the size of the operators as discussed in Ref. <cit.>. This class of effects is not addressed by the present formalism and must be treated separately] of bound states or resonances that couple to the asymptotic states. From the form factors and distribution functions, in turn, one can obtain various structural parameters, e.g. the charge or even the gluonic <cit.> radius of a given state.

The primary reason why a non-trivial formalism is required to extract multi-hadron observables from lattice QCD calculations is that the latter are performed in a finite spatial volume, usually a cube defined with periodic boundary conditions on the quark and gluon fields. This complicates the determination of scattering and transition amplitudes, because there is no simple relation between the finite-volume QCD eigenstates and the asymptotic multi-particle states that arise in the infinite-volume limit of the theory. However, in certain cases, it is possible to derive relations between finite- and infinite-volume observables. These have been implemented with great success to access a wide variety of scattering quantities directly from numerical LQCD calculations. See Ref. <cit.> for a recent review.

The most well-established such relation is that derived by Lüscher in Refs. <cit.> over three decades ago. In these seminal papers he showed how elastic two-particle scattering amplitudes can be extracted from the finite-volume energy spectrum below the lowest lying three- or four-particle threshold. Since then, the idea has been generalized to all possible two-body systems, in particular to multiple two-particle channels built form any number of particle species,  including particles with any intrinsic spin <cit.>. These formal ideas, together with significant algorithmic developments, have resulted in a proliferation of scattering amplitudes determined directly from lattice QCD <cit.>. 
A key limitation to the methods currently being implemented is the restriction to two-particle states, but the formal extension to three-particle systems has received significant attention recently and is progressing <cit.>.

Similar developments have resulted in both perturbative <cit.> non-perturbative <cit.> relations between finite-volume matrix elements and electroweak amplitudes. These have already been implemented in a variety of LQCD calculations <cit.>. In particular, Refs. <cit.> extracted the ργ^⋆→π electromagnetic form factor by determining the energy dependence of the corresponding amplitude, ππγ^⋆→π.

 

 

In this work we turn our attention to the prospect of determining 2 + 𝒥→2 transition amplitudes from finite-volume matrix elements. This was previously considered in Ref. <cit.>. In contrast to that work, here we restrict attention to kinematics such that only one two-particle channel is open, and take the two particles in the channel to be scalars. In addition, we only consider flavor-conserving external currents, so that the initial and final two-particle states contain the same particles.  Just as in Ref. <cit.>, the two-particle states are composed of QCD-stable (pseudo-)scalars. Relaxing these restrictions to provide a Lorentz-covariant formalism for any number of two-particle channels, including those with intrinsic spin, is expected to be straightforward.

As in Ref. <cit.> in this work we derive a mapping between finite-volume matrix elements of two-particle states and the infinite-volume  2 + 𝒥→2 amplitude.
The result is summarized by the flow-chart shown in Fig. <ref>. We find that, given the following quantities:


  * the two-particle finite-volume spectrum,

  * the 1 + 𝒥→1 form factors, 

  * the finite-volume two-particle matrix elements of 𝒥,

one can systematically constrain the 2 + 𝒥→2 amplitude in the kinematic window in which only the accommodated channel contributes. Our relation requires the generalized Lellouch-Lüscher factors <cit.>, that enter multiplicatively in the conversion, as well as
a new finite-volume function, denoted G, that appears in an additive correction, together with the single-particle form factor as well as the two-to-two scattering amplitude.

A simple limiting case of our result is the one in which the single-particle form factors vanish. In this limit the finite-volume artifacts associated with the G function also vanish and one recovers a Lellouch-Lüscher-like relation in which the correction factor appears twice, once each for the initial and final two-particle states.
However, when the single-particle form factors are nonzero, the term containing G is expected to give the dominant finite-volume effects. In particular, the analysis of Ref. <cit.> showed that, in the case of weak interactions, the finite-volume effects on the ground-state 2 + 𝒥→2 matrix element are given by an expansion in powers of 1/L.[In fact, the authors of Ref. <cit.> consider n + 𝒥→n matrix elements.]  In these expressions, the diagrams that appear as our G give 1/L^2 corrections while all other terms contribute with additional powers of 1/L.



 The purpose of this work is to improve certain technical details of the formalism and to provide more concrete information on the implementation procedure. We begin by providing a covariant version of the formalism in Sec. <ref>, where we also discuss three key examples involving ππ states. Then, in Sec. <ref>, we explain in detail our approach for evaluating G and outline why this is more challenging than the more-standard finite-volume functions relevant for two-to-two scattering. 

In addition to the standard threshold singularities, in this section we discuss and illustrate the presence of triangle singularities in the G function.  In Sec. <ref> we conclude and provide an outlook for future studies. Finally we include two appendices: Appendix <ref> gives various details relevant for the derivation of the improved formalism used here. Appendix <ref> includes various technical aspects regarding the evaluation of the finite-volume functions discussed in the main text.























§ COVARIANT REPRESENTATION OF THE FORMALISM


In this section we revisit the formalism derived in Ref. <cit.> and present a modified form in which all infinite-volume quantities are Lorentz covariant. 
We focus here only on the final result, and in Appendix <ref> we explain the (minimal) modifications to the original derivation that lead to this new form. 

This section is divided into three subsections: In Sec. <ref> we review the required notation and give the quantization condition, as well as the generalized Lellouch-Lüscher matrix, for two-particle states in a finite volume. Then, in Sec. <ref>, we provide a full description of our covariant 2 + 𝒥→2 formalism. Finally, in Sec. <ref>, we consider a handful of specific examples to show how our general expressions reduce for a particular system with specified quantum numbers.









 §.§ Kinematics and the quantization condition


 
We denote the 4-momentum of the incoming state in the finite-volume frame by P^μ_i ≡ (E_i , P_i) and that of the outgoing state by P^μ_f ≡ (E_f , P_f). The center-of-momentum (c.m.) energies corresponding to these are then given by

    E_i^⋆≡√(P_i^2) = √(E_i^2 - P_i^2) ,    and   
    E_f^⋆≡√(P_f^2) = √(E_f^2 - P_f^2) .

This also defines the metric used for the Minkowski-signature 4-vector dot products throughout. Following the notation of Ref. <cit.> we use ⋆ to indicate quantities defined in either the incoming or the outgoing c.m. frame. As explained below, we often use an i or an f index in addition to the ⋆, in order to completely specify the frame.

In this work we accommodate all values of 3-momenta allowed by the periodic boundary conditions, i.e. P_i = 2 πd_i/L and P_f = 2 πd_f/L where d_i and d_f are 3-vectors of integers. The energies and 3-momenta can differ between the initial and final states due to the momentum carried by the external current, P_f^μ - P_i^μ [see Fig. <ref>(a)]. The physical quantities discussed below depend on Lorentz scalars. For most systems we will primarily be sensitive to spacelike values of the momentum transfer, motivating us to introduce

    Q^2 ≡  -(P_f - P_i)^2  ,

which is positive for spacelike P_f^μ - P_i^μ.


As mentioned above, we restrict our attention here to values of E_i^⋆ and E_f^⋆ such that only a single two-particle channel can propagate. Within the single channel considered, we accommodate both identical and non-identical scalars and allow these to be non-degenerate in the latter case, with physical masses m_1 and m_2. We assume, however, that the current, 𝒥, is flavor conserving so that the same two particles appear in the initial and final state.[Given the results presented below, implementing the covariant modification to the multi-channel expressions of Ref. <cit.> with flavor-changing currents should be straightforward, albeit tedious and likely leading to index-heavy notation.]



We now turn to the kinematic variables describing individual particles within the two-particle channel. For the remainder of this subsection, take P^μ≡ (E, P) to simultaneously represent the initial and final state 4-momenta. Denoting the 3-momentum of particle 2 (with mass m_2) by k, the corresponding on-shell 4-vector is k^μ = (ω_k2, k), where

    ω_ k2≡√( k^2 + m_2^2) ,

with k=|k|. 

Next note that, in order to satisfy the specified total 4-momentum (P^μ), particle 1 must carry P^μ - k^μ≡ (E - ω_k2, P - k). Thus, for general k, one cannot simultaneously require that the particle momenta sum to P^μ and that particle 1 is on shell. The latter holds only when the temporal component, E - ω_k2, is equal to

    ω_P k1≡√((P -k)^2 + m_1^2) .


To better understand when the on-shell condition (E - ω_k2 = ω_P k1) is satisfied, it is useful to introduce [Λ_β]^μ_ν≡Λ^μ_ν(β) as a boost matrix with boost velocity β≡P/E. We then define

    k^⋆μ≡ (ω_k2^⋆, k^⋆) ≡[Λ_β]^μ_ν k^ν ,

and observe

    ω_k2^⋆ = √(k^⋆ 2 + m_2^2) ,

where k^⋆≡|k^⋆|.


[It is worth emphasizing that the definitions of k^⋆ and k^⋆ depend on (E, P), k and m_2 but not on m_1. This asymmetry in the definition is removed when both particles are on shell.] 


By contrast, the 4-momentum of particle 1 boosts to

    P^⋆μ - k^⋆μ = (E^⋆ - ω_k2^⋆, - k^⋆) ≡[Λ_β]^μ_ν (P^ν - k^ν)  ,

where we have used [Λ_β]^μ_ν P^ν = (E^⋆, 0). We deduce that the c.m. frame on-shell condition is E^⋆ = ω_k1^⋆ +  ω_k2^⋆ where

    ω_k1^⋆≡√(k^⋆ 2 + m_1^2) .


The advantage of working in this frame is that the on-shell condition reduces to a constraint on the value of k^⋆. In particular the particle pair is on shell if and only if k^⋆ = q^⋆ with the latter defined by

    E^⋆≡√(q^⋆ 2 + m_1^2) + √(q^⋆ 2 + m_2^2) .

Finally, the initial or final state-indices must be applied to all of these quantities once the total 4-momentum is associated with a particular state. In particular if we take P^μ→ P_i^μ, then the corresponding quantities above become ω_P_i k2, [Λ_β_i]^μ_ν, k^⋆μ_i, k^⋆_i, k^⋆_i, ω_k1i^⋆, ω_k2i^⋆, E^⋆_i, q^⋆_i. The only quantities that do not inherit a frame index are the finite-volume frame momentum, k, and the corresponding energy, ω_k2.



With this notation in hand we now give the quantization condition for two scalar particles in a finite volume. This is written as a determinant condition involving the on-shell two-particle scattering amplitude, ℳ(P^2), represented as a diagonal matrix in angular-momentum space. For any fixed values of P and L, the finite-volume energy spectrum is given by solutions to  <cit.>

    [ℳ(P^2)^-1 + F(P,L) ]= 0  .

This holds in the region 0 < E^⋆ < E^⋆_th, where E^⋆_th is the energy of the lowest-lying multi-particle threshold that we ignore, which could be a two-, three-, or four-particle threshold. The relation holds up to neglected corrections of the form e^- m L, where m is the physical mass of the lightest degree of freedom in the theory.

The precise definition of ℳ (P^2) is given by

    ℳ_ℓ m ;   ℓ' m'(P^2) = δ_ℓℓ'δ_m m'ℳ^(ℓ)(P^2)  ,    with   ℳ^(ℓ)(s) = 8 π E^⋆/ξ1/q^⋆δ^(ℓ)(q^⋆) - i q^⋆|_ E^⋆ 2 = s ,

where δ^(ℓ)(q^⋆) is the scattering phase shift, and ξ = 1/2 if the particles are identical and ξ=1 otherwise. 

The remaining ingredient is F(P,L), a matrix of finite-volume geometric functions defined by

    F_ℓ m ;   ℓ' m'(P,L)  ≡ξ [1/L^3∑_𝐤 - ∫d^3 k/(2π)^3 ]
     𝒴_ℓ m(k^⋆ )   𝒴^*_ℓ' m'(k^⋆)   
      /2 ω_Pk1 2 ω_k2(E-ω_k2 -ω_Pk1 + i ϵ ) ,

 where 
 
    𝒴_ℓ m(k^⋆)    ≡√(4 π)  Y_ℓ m(k̂^⋆) (k^⋆/q^⋆)^ℓ .

In Appendix <ref> we review an efficient method to evaluate this, based on analytic expressions for the integrals defined using the cutoff functions introduced in Ref. <cit.>. 

To close this subsection, we introduce one additional finite-volume matrix

    ℛ(E_n, P) ≡lim_E → E_n[   (E - E_n)/F^-1(P,L) + ℳ(P^2)]  .

This object, introduced in this form in Ref. <cit.>, is the generalization of the Lellouch-Lüscher factor <cit.> that relates finite-volume matrix elements of two-particle states to the corresponding infinite-volume decay and transition amplitudes. 













 §.§ Relating finite-volume matrix elements with 2 + 𝒥→2 transition amplitudes


We are now ready to present our improved finite-volume formalism. The approach that we advocate here differs from that of Ref. <cit.> in two key ways. 

First, the separation of singularities, required to disentangle finite-volume effects in the 2 + 𝒥→2 amplitude, is done here using Lorentz invariant poles of the form 1/(k^2 - m^2). In the previous work we instead used 1/[2 ω_k(k^0 - ω_k)]. As long as we consistently modify the pole form everywhere, it turns out that either choice is valid. The advantage of the present approach is that it ensures all infinite-volume quantities are Lorentz covariant and also simplifies the form of the new finite-volume function, G, arising due to the triangle diagram shown in Fig. <ref>(b) and defined explicitly in Fig. <ref>(d) and in Eq. (<ref>) below.

Second, we treat the single-particle matrix elements in a simpler way here than we did in Ref. <cit.>. Our approach requires isolating the 1 + 𝒥→1 matrix element in order to express the finite-volume effects of the triangle diagram [Fig. <ref>(b)]. In our previous publication, a complicated scheme was presented to project the matrix element on shell. Though correct, we have come to realize that this procedure is unnecessary. The reason, as we explain in more detail below, is that one can decompose the matrix elements into kinematically determined tensor structures and form factors. Projecting the kinematic factors on shell is unnecessary, and removing this step gives a simpler approach that leads to the same infinite-volume observables.



We begin by introducing notation for the physical 2 + 𝒥→2 matrix element that we are after [see also Fig. <ref>(c)]

    _μ_1 ⋯μ_n(P_f, k'; P_i, k) ≡⟨   P_f, k'; out|𝒥_μ_1 ⋯μ_n(0)  |  P_i, k; in⟩_ conn. .

Here the initial and final states are standard two-particle asymptotic states with the usual relativistic normalization convention and 𝒥_μ_1 ⋯μ_n(0) is a generic local current insertion. As is shown in Fig. <ref>(a), the initial state is built from particle 1 [with mass m_1 and on-shell momentum (P_i - k)^μ satisfying (P_i-k)^2 = m_1^2] and particle 2 [mass m_2, momentum k^2 = m_2^2]. The final state is built from the same pair, now carrying momenta P_f-k' and k'. Following the discussion of the previous sub-section, the c.m. frame 3-momenta are denoted by k_i^⋆ and k_f^⋆ and have magnitudes equal to q_i^⋆ and q_f^⋆ respectively. The label “conn.” emphasizes that only fully connected diagrams, those shown in Fig. <ref>(c), contribute to the definition of the 2 + 𝒥→2 matrix element.

A consequence of the on-shell constraints is that, once total energy and momenta are fixed, the two-particle states only have directional degrees of freedom, k̂^⋆_i and k̂^⋆_f. However, in contrast to the scattering amplitude, ℳ, for  a decomposition in spherical harmonics is not useful. The directional degrees of freedom sweep across pole singularities due to the diagrams in the second line of Fig. <ref>(c), implying that the decomposition is only defined in the sense of distributions. More importantly, these long-distance singularities guarantee that higher partial waves will not be suppressed.   

The issue is easily resolved by removing the singular terms before decomposing in harmonics. This was already discussed in detail in Ref. <cit.> where the quantity _ was first introduced. In this work we define an alternative, Lorentz-covariant version of the divergence-free amplitude with the same symbol [see again the second line of Fig. <ref>(c)][This subtraction assumes that the current couples only to particle 1. In the case that the current couple to both particles, two additional terms must be subtracted. In particular, if the particles are identical one must always subtract four terms in which the propagators carry the four possible values of external momenta. See Eq. (<ref>) as well as Ref. <cit.> for explicit expressions.]

    _ df;μ_1 ⋯μ_n≡_μ_1 ⋯μ_n - 
    i  (P_f,k',k) i/(P_f-k)^2 - m_1^2_μ_1 ⋯μ_n
    -
     _μ_1 ⋯μ_ni/(P_i - k')^2 - m_1^2
    i '(P_i, k',k)
       ,

where _μ_1 ⋯μ_n is the single-particle matrix element of the external current, defined in detail below, and

    (P_f,k',k)    ≡∑_ℓℳ^(ℓ)(P_f^2)   ( 2 ℓ + 1  ) P_ℓ (k̂'^⋆_f ·k̂_f^⋆)   ( k^⋆_f/q^⋆_f)^ℓ , 
    '(P_i, k',k)    ≡∑_ℓℳ^(ℓ)(P_i^2)   ( k'^⋆_i/q^⋆_i)^ℓ  ( 2 ℓ + 1  ) P_ℓ (k̂'^⋆_i ·k̂_i^⋆)   .

Here P_ℓ (cosθ) are the standard Legendre polynomials, satisfying ∑_m=-ℓ^ℓ 4 π Y^*_ℓ m(x̂) Y_ℓ m(ŷ) = (2 ℓ + 1) P_ℓ(x̂·ŷ). 
Unlike , _ df admits a uniformly convergent decomposition in spherical harmonics

    _ df;μ_1 ⋯μ_n(P_f, k'; P_i, k) ≡ 4π Y^*_ℓ' m'( k̂'^⋆_f   )    _ df;μ_1 ⋯μ_n; ℓ' m'; ℓ m( s_f, s_i, Q^2)         Y_ℓ m( k̂^⋆_i   )      ,

 where s_i ≡ P_i^2 = E_i^⋆ 2, s_f ≡ P_f^2 =  E_f^⋆ 2, and the repeated harmonic indices on the right-hand side are summed.

The subscript “df”, short for divergence-free, refers only to the lack of kinematic singularities arising from a long lived state between the 2→2 and 1 + 𝒥→1 transitions. In particular, as we show in Sec. <ref> below, _ df does have so-called triangle singularities as a function of E_i^⋆ and E_f^⋆ associated with diagram of the kind depicted in Fig. <ref>(b). These pose no problem to the harmonic decomposition but must, of course, be understood in order to successfully extract and interpret both _ df and .

The scattering amplitude, ℳ, is defined with powers of (k^⋆/q^⋆), referred to as barrier factors. These must be included due to the manner in which the factors of w_μ_1 ⋯μ_n and ℳ arise in the triangle diagram of Fig. <ref>(b). In particular, the loop is summed over all finite-volume momenta in such a way that the current insertion, w_μ_1 ⋯μ_n, as well as the external factors of ℳ, are sampled at off-shell values of momenta (i.e. values for which k^2 ≠  m_i^2). Nonetheless, the power-law finite-volume effects are governed by the on-shell values of w_μ_1 ⋯μ_n and ℳ only, and the off-shell contributions can be absorbed into other infinite-volume quantities. The catch here is that the on-shell projection, effected via k^⋆→ q^⋆, amounts to replacing k^⋆ℓY_ℓ m(k̂^̂⋆̂) with q^⋆ℓY_ℓ m(k̂^̂⋆̂). The latter has spurious singularities near k^⋆ = 0 and thus more care is needed. The inclusion of barrier factors resolves the issue.



The single-particle matrix element, _μ_1 ⋯μ_n, appearing in Eq. (<ref>), is a function of (P_f - k, P_i - k) in the first term and (P_f - k', P_i - k') in the second. Using the first term as a reference, the explicit definition is given in a three step processes. First,
the fully on-shell version is defined via a single-particle matrix element

    ^on_μ_1 ⋯μ_n (P_f-k,P_i-k) ≡⟨ P_f-k, m_1 |𝒥_μ_1 ⋯μ_n(0) | P_i - k, m_1 ⟩ .

Second, this is formally continued to an off-shell quantity in the context of some generic effective field-theory. The latter object is then decomposed into kinematic tensors and form factors

    ^off_μ_1 ⋯μ_n (P_f-k,P_i-k)
       =∑_j K^(j)_μ_1 ⋯μ_n(k , P_f,P_i) f^(j)[Q^2; (P_f - k)^2; (P_i-k)^2]  ,

where the sum runs over all possible tensor structures for the given current.
Third, and finally, a partial on-shell projection is performed to define the version of _μ_1 ⋯μ_n that appears in Eq. (<ref>). In this step we set the virtualities within the form factors to be on shell [(P_f - k)^2, (P_i-k)^2 ⟶ m_1^2] and also set k^0 = ω_k2 everywhere. We reach

    _μ_1 ⋯μ_n (P_f-k,P_i-k)    ≡∑_j K^(j)_μ_1 ⋯μ_n(m,k,P_f,P_i)   |_k^0=ω_k2 f^(j)(Q^2) .

This definition is not completely on shell because, within 
K^(j)_μ_1 ⋯μ_n only, the virtualities (P_i-k)^2 and (P_f-k)^2 may differ from m_1^2.

In what follows we will consider sums and integrals over the spectator momentum, k. With this in mind, it is convenient to rewrite K^(j)_μ_1 ⋯μ_n as a sum of terms that isolate the dependence on this quantity

    K^(j)_μ_1 ⋯μ_n(m,k,P_f,P_i) ≡∑_n'=0^n   K^ω_μ_1⋯μ_n'(m,k)    C^(j)_μ_n'+1…μ_n(P_f,P_i)  ,

where C^(j)_μ_n'+1…μ_n(P_f,P_i) has no indices for n'=n and 

    K^ω_μ_1⋯μ_n(m,k)  ≡ k_μ_1⋯ k_μ_n  |_k_0=ω_k2 .

Here the superscript ω indicates that the 4-momenta in K^ω are on shell. This will be important in Sec. <ref> below, where we introduce various off-shell versions of K in our formulas for evaluating G.

 Having defined all of the infinite-volume quantities that enter our formalism, we now turn our attention to the finite volume. As mentioned in the introduction, we restrict attention to a finite cubic spatial volume, with periodicity L applied to the fields in each of the three spatial directions. 
In this set-up, we consider a matrix element in which the local current 𝒥 is sandwiched between two finite-volume states, each of which has the quantum numbers of the two-particle system. 
As we demonstrate in Appendix <ref>, following the derivation of Ref. <cit.>, this LQCD observable is related to the infinite-volume 2 + 𝒥→2 transition amplitude via

    |  ⟨ P_f, L |𝒥^μ_1 ⋯μ_n(0)  |  P_i, L ⟩ |^2 
    =1/L^6Tr[
     ℛ( P_i) 
    ^μ_1 ⋯μ_n(P_i,P_f,L)
    ℛ( P_f)  
    ^μ_1 ⋯μ_n(P_f,P_i,L)
    ]  ,

 where  directly determines  via

    ^μ_1 ⋯μ_n(P_f,P_i,L)  - ^μ_1 ⋯μ_n(s_f,s_i,Q^2) ≡∑_j
      ∑_n'=0^n 
      C^(j), μ_n'+1…μ_n(P_f,P_i) f^(j)(Q^2)
    
    ×ℳ(s_f)  G^μ_1 ⋯μ_n'(P_f,P_i,L)  ℳ(s_i)   .
 
 In the case of distinct particles that both couple to the current one must subtract two terms. In the second the particle labels 1 and 2 are swapped everywhere in the definition of G. In addition, the f^(j)(Q^2) will take on different values if the particle species differs. These should be given a species label in the case that the current couples to both. See Eq. (<ref>) as well as Ref. <cit.> for explicit expressions.
  
 In these expressions we have suppressed angular-momentum indices on ℛ, ^μ_1 ⋯μ_n, ^μ_1 ⋯μ_n, ℳ and G^μ_1 ⋯μ_n'. Each object carries the set ℓ' m'; ℓ m and these are contracted between adjacent factors in the usual matrix multiplication. The trace in Eq. (<ref>) is also over this angular-momentum space. Note, by contrast, that the index set μ_1 ⋯μ_n is not summed but rather fixed to common values for all objects appearing in these equations.
 
Finally, the matrix G^_μ_1 ⋯μ_n'(P_f,P_i,L) is defined diagrammatically in 
 Fig. <ref>(d) and has the explicit form

    G_μ_1 ⋯μ_n; ℓ_f m_f ;     ℓ_i m_i (P_f,P_i,L)
       ≡ [1/L^3∑_𝐤 - ∫d^3 k/(2π)^3 ]
    𝒴_ℓ_f m_f(k^⋆_f)  
     D(m,k) K^ω_μ_1⋯μ_n(m,k)   𝒴^*_ℓ_i m_i(k^⋆_i)   
    
       ,

where

    D(m,k)   ≡1/2ω_k21/  (P_f-k)^2 - m_1^2 + i ϵ1/  (P_i-k)^2 - m_1^2 + i ϵ|_k_0=ω_k2 .

This differs from the form presented in Ref. <cit.> due to the aforementioned modifications: The poles are Lorentz invariant and the 1 + 𝒥→1 matrix element is expressed in terms of tensor structures leading to K^ω. Note that the modifications to G are directly connected to those in the definition of , Eq. (<ref>). We have altered these two intermediate objects in such a way that  is unchanged. 









 §.§ Examples


In this final subsection we show how the construction outlined above may be applied to specific, phenomenologically well-motivated examples.



  §.§.§ (π^+ π^0)_J=1 + j_μ→ (π^+ π^0)_J=1 



We begin with the electro-magnetic form factors of a charged ρ meson. The ρ decays predominantly to the vector-isovector ππ state. Indeed for heavier than physical light-quark masses (such that 4 M_π > M_ρ) and in the iso-spin symmetric theory, this is the only possible QCD decay channel. If the light-quark mass is further increased, the two-pion threshold eventually exceeds the ρ mass (2 M_π > M_ρ) and the latter becomes a stable particle.
In this case, one can extract the form factors of the ρ directly from finite-volume matrix elements. 
See, for example, Ref. <cit.>. 

To determine the analogous observable at quark masses for which the ρ→ππ decay occurs, it is necessary to first extract the (π^+ π^0)_J=1 + j_μ→ (π^+ π^0)_J=1 amplitude for a wide range of kinematic points. As depicted in Fig. <ref>, by fitting the dependence of the initial- and final-state energies to a functional ansatz, one may analytically continue these amplitudes to the complex-valued pole to obtain the ρ form factors. Detailing the steps of this continuation will be the focus of future work. Here we focus on the extraction of the (π^+ π^0)_J=1 + j_μ→ (π^+ π^0)_J=1 amplitude for real ππ energies.

By interpolating isospin-one initial and final states (I=1,m_I = ± 1), we project the system to a sector where all even angular momenta vanish. Thus, regardless of the values of P_f and P_i, we will always have a finite-volume irrep that contains J^P = 1^-, with the next contamination coming from J ≥ 3. Taking the latter to be negligible, we can approximate total angular momentum as a good quantum number. However, even in this simple case, the azimuthal component, m_J, is not a good quantum number in the finite volume. In general the positive and negative helicity states mix, but one can readily construct linear combinations of these that are invariant under transformations of the cubic group or its little groups <cit.>. 

Considering first the incoming state, we restrict attention to a specified set of finite-volume quantum numbers: Λ_i and μ_i, labeling the irrep and row respectively, of the little group defined by P_i. In addition we assume that, within this irrep, the interactions are dominated by ℓ_i = 1 and a particular m_i value so that the matrices can be truncated to single entries.  
Doing the same for the final-state, Eq. (<ref>) becomes

[Note, the procedure for subducing the matrix elements onto the appropriate symmetry group is discussed in detail in Ref. <cit.>. Although some of the details of the formalism has changed, this aspect remains the same.]

    𝒲_L, df;Λ_fμ_f ; Λ_iμ_i^μ(P_f,P_i,L) =  e^ i δ^I=1_ππ(q^⋆_f)⟨ E_f,n', P_f, L ;Λ_fμ_f|    j^μ(0)  |  E_i,n, P_i, L;Λ_iμ_i ⟩/|ℛ^I=1_Λ_fμ_f(P_f, L)  ℛ^I=1_Λ_iμ_i(P_i, L) |^1/2  e^   i δ^I=1_ππ(q^⋆_i) ,

where δ^I=1_ππ is the elastic scattering phase shift for ππ→ππ in the ρ channel and

    |ℛ^I=1_Λμ(P, L) |^1/2≡1/L^3√(q^⋆/16 π E^⋆) [ E^⋆/E∂/∂ E^⋆ ( ϕ_Λμ^d(E^⋆, L) + δ^I=1_ππ(q^⋆)  )  ]^-1/2 ,

with ϕ_Λμ^d(E^⋆, L) defined for the irreps that couple to ℓ=1 ππ in, for example, Ref. <cit.>. We stress that all instances of E_f and E_i in these expressions are to be evaluated at any pair of finite-volume energies, E_f,n'(L) and E_i,n(L) respectively, satisfying Eq. (<ref>). 

In Eq. (<ref>), 𝒲_L, df;Λ_fμ_f ; Λ_iμ_i^μ is the subduced  version of 𝒲_L, df^μ. As discussed in Ref. <cit.>, this can be obtained from 𝒲_L, df;ℓ_fm_f ; ℓ_im_i^μ by rotating these into the helicity basis and then using the subduction matrices <cit.>. 
This subduction procedure requires no approximation. However each irrep couples to an infinite tower of partial waves, and only by neglecting these above a certain maximum value does one reach useful expressions.

 

Compared to Eq. (<ref>), in Eq. (<ref>) we have dropped the trace, since we are ignoring all but one partial wave, and have solved for ^μ. For the latter step there is a potential sign ambiguity that one must address. Note that ℛ^I=1_Λμ(P, L) = |ℛ^I=1_Λμ(P, L) | e^- 2 i δ^I=1_ππ(q^⋆) as is shown, for example, in Eqs. (132)-(134) of Ref. <cit.>. The phases in ℛ^I=1_Λμ(P_f, L) and ℛ^I=1_Λμ(P_i, L) precisely generate the Watson phases within ^μ(P_f,P_i,L) as they must, since the finite-volume matrix element is real. The remaining ± ambiguity is constrained by the known value in the Q^2 → 0 limit, but in certain cases a remaining ambiguity may survive.



The final step is to convert this to the infinite-volume quantity ^μ via

    𝒲^μ_ df;Λ_fμ_f ; Λ_iμ_i(s_f,s_i,Q^2) 
       =  
    𝒲^μ_L, df;Λ_fμ_f ; Λ_iμ_i(P_f,P_i,L)  
       
    -    f_π^+(Q^2)  ℳ^I=1(s_f)   [ 
      (P_f + P_i)^μ  G_Λ_fμ_f ; Λ_iμ_i(P_f,P_i,L)  - 2   G^μ_Λ_fμ_f ; Λ_iμ_i(P_f,P_i,L)     ]  ℳ^I=1(s_i)  .

The general form of Eq. (<ref>) is overly complicated for this application but still applies with f^(1)(Q^2) = f^(2)(Q^2) = f_π^+(Q^2) corresponding to the usual (spacelike) pion form factor. Here we have also used the standard relation

    ⟨ P_f - k; M_π| j_μ(0) | P_i - k; M_π⟩ = (P_f + P_i - 2 k)_μ f_π^+(Q^2)  .

In these expressions we are neglecting the electro-magnetic form factor of the neutral pion which is expected to be small but non-zero for Q^2 ≠ 0. Finally we remark that a factor of i may appear in this relation depending on the conventions for Euclidean or Minkowski gamma matrices. As the same gamma factors appear in all terms of our formalism changing conventions just amounts to multiplying both sides of μ-indexed equations by a common factor.





  §.§.§ (π^+ π^0)_I_i + j_μ→ (π^+ π^0)_I_f 


The electromagnetic current has both I=0 and I=1 components, but G parity guarantees that all matrix elements of the iso-scalar part between two ππ states must vanish. If we take the angular momentum to be unconstrained then the incoming ππ state may, in general carry isospin I_i=0,1,2. The iso-vector current then couples the fixed incoming isospin as follows:
0 → 1, 1 → 0,1,2 and 2 → 1,2. If we further restrict attention to | m_I | = 1 states (i.e. π^+ π^0 states) then this reduces to 1 → 1,2 and 2 → 1,2. An alternative way to distinguish these possibilities is by fixing orbital angular momentum: J=0 ⇔ I=2 and J=1 ⇔ I=1. 

In this way we identify four possible transitions involving π^+ π^0 states, the p-wave to p-wave matrix element considered above as well as s to s, s to p and p to s. It turns out that all four of these transitions are described by Eqs. (<ref>) to (<ref>) provided that we can neglect the effects of finite-volume mixing with J ≥ 2 states. The only modifications are that one must generally project to different irreps Λ to access the J=0 components and that the scattering amplitudes must correspond to the isospin of the state















    𝒲_L, df;Λ_fμ_f,Λ_iμ_i^μ(P_f,P_i,L) =  e^ i δ^I_f_ππ(q^⋆_f)⟨ E_f,n', P_f, L ;Λ_fμ_f|    j^μ(0)  |  E_i,n, P_i, L;Λ_iμ_i ⟩/|ℛ^I_f_Λ_fμ_f(P_f, L)  ℛ^I_i_Λ_iμ_i(P_i, L) |^1/2  e^   i δ^I_i_ππ(q^⋆_i) ,

with

    𝒲^μ_ df;Λ_fμ_f,Λ_iμ_i(s_f,s_i,Q^2) 
    =  
    𝒲^μ_L, df;Λ_fμ_f,Λ_iμ_i(P_f,P_i,L)  
    
    -    f_π^+(Q^2)  ℳ^I_f(s_f)   [ 
      (P_f + P_i)^μ  G_Λ_fμ_f,Λ_iμ_i(P_f,P_i,L)  - 2   G^μ_Λ_fμ_f,Λ_iμ_i(P_f,P_i,L)     ]  ℳ^I_i(s_i)  .




  §.§.§ Gluonic structure


Thirty years ago, Jaffe and Manohar identified a structure function that provides a measure of the gluon distribution within hadrons <cit.>. This has since lead to lattice calculations of the leading moments of these distributions, for example within the ϕ-meson. Thus far, the calculations are restricted to heavy quark masses where the ϕ is stable within QCD <cit.>. Similarly, calculations of gluonic moments for light nuclei are already underway, again for values of the light quark masses that lead to the nuclei being deeply bound <cit.>.  The formalism presented here will allow for future calculations closer to the physical point by accommodating the finite-volume effects of loosely-bound as well as resonant states. 

In order to extract the leading moment of the gluon structure function, one must evaluate the traceless part of the product of two gluon energy-momentum tensors [𝒪_μν∼𝒢_μβ 𝒢_βν]. As one might expect, this is more complicated than the case considered above in part because it is a rank-two tensor. 

A starting point in extracting gluonic moments of resonances from LQCD would likely be to consider the ρ, discussed in Sec. <ref> above.

In this case, the relation between the finite-volume matrix elements and the transition amplitude 
is very similar to Eq. (<ref>). The only distinction arises in relating 𝒲^μν_ df and 𝒲^μν_L, df. To do so one must determine the scalar (G), vector (G^μ), and tensor (G^μν) contributions to the finite-volume G-function and combine these with the relevant gluonic form-factors of the single-pion state. 
 


This concludes our discussion of the covariant formalism. 
The aim of the section was to provide a procedure by which three inputs: (i) single-particle form factors, (ii) 2→2 scattering amplitudes and (iii) finite-volume kinematic functions, can be combined with finite-volume two-particle matrix elements to extract the infinite-volume 2 + 𝒥→2 transition amplitudes. In this recipe the ingredient that remains most obscure is the new finite-volume function G^μ_1 ⋯μ_n(P_f,P_i,L), defined in Eq. (<ref>). Thus, in the next section, we give a detailed description of how this can be efficiently evaluated numerically.





















§ EVALUATING G(P_F,P_I,L)


Our aim is to evaluate 

    G_σ (P_f,P_i,L)
       ≡ [1/L^3∑_𝐤 - ∫d^3 k/(2π)^3 ]
    𝒴_ℓ_f m_f(k^⋆_f)  
     D(m,k) K^ω_μ_1⋯μ_n(m,k)   𝒴^*_ℓ_i m_i(k^⋆_i)   
    
       ,

where we have introduced the collective index σ≡ [μ_1 ⋯μ_n;  ℓ_f m_f ;   ℓ_i m_i]. We also take the convention that if σ is written as a low (high) index, then all of the Lorentz indices it contains are also understood to be low (high). The sum is straightforward to calculate numerically once a cutoff function has been included. 
We comment here that the ultraviolet divergences match between the sum and the integral, meaning that the difference has an unambiguous limit as the cutoff is removed and thus that G^σ is a universal quantity with no scheme dependence. 

Evaluating the integral part of G^σ turns out to be significantly more challenging. 
The integrand contains singularities associated with on-shell intermediate states and, although these are perfectly integrable given the i ϵ pole prescription,
numerical evaluations converge very slowly for standard numerical techniques. Thus it is highly advantageous to find analytical representations to the extent possible. 

For the case of P_i = P_f, it turns out that one can provide exact analytical expressions for the integral, as discussed in the following subsection. For the generic case, with P_i ≠ P_f, we have not managed to obtain fully analytic results. Instead, we write the integral as the sum of two terms. The first includes all singularities and can be evaluated analytically to the level of Feynman parameters. The second is defined with a smooth integrand such that standard numerical integration is effective. Our approach for evaluating G^σ in the case of P_i ≠ P_f case is detailed in Sec. <ref> below.









 §.§ P_i = P_f


The function G^σ simplifies considerably when the initial and final momenta coincide, P_i = P_f ≡ P, i.e. when the momentum transfer vanishes, P_f^μ - P_i^μ = 0. A particularly helpful feature of these kinematics is that a natural preferred frame emerges: the simultaneous rest frame of the incoming and outgoing particles, in which the spatial part of P vanishes. 

Another consequence of P_i = P_f is that the product of poles within D(m,k) becomes a double pole of the form

    D(m,k) = 
    1/2ω_k2(1/  (P-k)^2 - m_1^2 + i ϵ)^2_k^0 = ω_k2 .

Focusing on the factor within parenthesis, note that this can be rewritten as

    1/  (P-k)^2 - m_1^2 + i ϵ  =   ω_q2^⋆/  E^⋆ 1/(q^⋆2-k^⋆2+ i ϵ)   
       - 1/  4 ω_q2^⋆  E^⋆+ 𝒪(q^⋆ 2-k^⋆ 2)   ,

where k^⋆ and q^⋆ are defined with respect to the P rest frame.

If D(m, k) were to only contain a single pole, then, after acting with the sum-integral difference, only the leading-order, singular term would be relevant. This is because the sum-integral difference of smooth functions 
leads only to exponentially suppressed volume dependence that we neglect. However, in this case the first two terms in the expansion are important as they generate double and single poles upon squaring. This leads to

    G_σ (P,P,L)
     = 
       [1/L^3∑_𝐤 - ∫d^3 k/(2π)^3]  1/2ω_k2𝒴_ℓ_f m_f(k^⋆) 
      K^ω_μ_1⋯μ_n(m,k)     𝒴^*_ℓ_i m_i(k^⋆)      D(m, k^⋆, q^⋆)  ,

where

    D(m, k^⋆, q^⋆) ≡ω_q2^⋆ 2/  E^⋆ 21/(q^⋆2-k^⋆2+ i ϵ)^2   
       -
         1/ 2 E^⋆21/(q^⋆2-k^⋆2+ i ϵ) .


At this stage it is useful to decompose the angular dependence within the tensors into a single set of spherical harmonics

    𝒴_ℓ_f m_f(k^⋆) 
      K^ω_μ_1⋯μ_n(m,k)     𝒴^*_ℓ_i m_i(k^⋆) 
      ≡√(4 π)/(q^⋆ )^ℓ_i + ℓ_f∑_J M𝒞_σ, JM(β , k^⋆)     k^⋆ JY_J M(k̂^⋆)  .

As we explain in Appendix <ref>, 𝒞_σ, J M (β, k^⋆) can be efficiently calculated by writing the factors within K^ω as boosted c.m. frame vectors, k^μ = [Λ_-β]^μ_ν k^⋆ν. Such factors can then be written as spherical harmonics and, using Clebsch-Gordon coefficients, these can be combined with the external factors of 𝒴_ℓ_f m_f(k^⋆) and  𝒴^*_ℓ_i m_i(k^⋆) to identify the a final harmonic basis. As a final note, we stress that it is possible to unambiguously separate the dependence on β = P/E and k^⋆ within 𝒞_σ, i.e. one can vary k^⋆ while holding β constant. This will be important for the manipulations performed below.

The construction of  𝒞_σ, JM (β, k^⋆) is discussed in detail in Appendix <ref>. As a specific simple example, here we consider the case σ = [μ; 10; 10]. The numerator within G_σ is then

    (q^⋆)^2𝒴_10(k^⋆) 
    k_μ𝒴^*_10(k^⋆) = 3 (k^z ⋆ )^2  [Λ_-β]_μ^ν k_ν^⋆ .

The current insertion, k_μ, can be written as a combination of ℓ=0 and ℓ=1 spherical harmonics. Combining this with the two ℓ=1 harmonics from the external states, one finds that 𝒞_σ, JM is zero for J>3. 

The J M = 0 0 component only has a non-zero contribution in the k_0^⋆ component. Isolating this contribution and substituting the definition of the boost matrix, we reach

    𝒞^σ_00(β, k^⋆) =
    k^⋆2ω_k2^⋆P^μ/E^⋆ .

The remaining nonzero coefficients, arising for J ≤ 3, take a more complicated form in general, but simplify for P∝ẑ. If we additionally focus on the μ=z component, then only three coefficients survive

    𝒞^σ_10(β, k^⋆)  =3√(3)/5E/E^⋆  k^⋆ 2 ,       𝒞^σ_20(β, k^⋆)  =2√(5)/5P^z/E^⋆ω _k2^⋆ ,       𝒞^σ_30(β, k^⋆)  =6 √(7)/35E/E^⋆ .

We revisit this case below to show how it enters our final construction for G^μ_  10 ;10.  

Returning to the main line of argument, we now substitute Eq. (<ref>) into the expression for G_σ (P,P,L), Eq. (<ref>), to reach

    G_σ (P,P,L) =  √(4 π)/(q^⋆ )^ℓ_i + ℓ_f∑_J M [1/L^3∑_𝐤ω_k2^⋆/ω_k2 - ∫d^3 k^⋆/(2π)^3]  𝒞_σ, JM(β , k^⋆)/2ω_k2^⋆    k^⋆ JY_J M(k̂^⋆)   D(m, k^⋆, q^⋆)       .


The next step is to expand the k^⋆ dependence in the first factor within the summand about the pole location 

    𝒞_σ, JM(β , k^⋆)/2ω_k2^⋆ = 𝒞_σ, JM(β , q^⋆)/2ω_q2^⋆   - (q^⋆ 2 - k^⋆ 2) ∂/∂ q^⋆ 2𝒞_σ, JM(β , q^⋆)/2ω_q2^⋆  + 𝒪 [(q^⋆ 2 - k^⋆ 2)^2  ]  .

In the second term we have rewritten the derivative with respect to k^⋆ 2 (to be evaluated afterwards at k^⋆ 2 = q^⋆ 2) directly as a derivative with respect to q^⋆ 2. This is possible because 𝒞_σ, JM(β , k^⋆) has no implicit q^⋆ 2 dependence and it is formally possible to vary q^⋆, and thus E^⋆, while holding β constant. 

Combining this with the definition of D(m, k^⋆, q^⋆) one can show

    𝒞_σ, JM(β , k^⋆)/2ω_k2^⋆D(m, k^⋆, q^⋆) =  ω_q2^⋆/ 2 E^⋆ 2𝒞_σ, JM(β , q^⋆)/(q^⋆2-k^⋆2+ i ϵ)^2   
    -  [  ( 1/ 2 E^⋆2  + ω_q2^⋆ 2/  E^⋆ 2∂/∂ q^⋆ 2 )  𝒞_σ, JM(β , q^⋆)/2ω_q2^⋆ ]  1/q^⋆2-k^⋆2+ i ϵ
    
        + 𝒪 [(q^⋆ 2 - k^⋆ 2)^0  ] .

Remarkably, the operator in parenthesis vanishes when acting on 1/ω_q2^⋆ so that this reduces to

    𝒞_σ, JM(β , k^⋆)/2ω_k2^⋆D(m, k^⋆, q^⋆)     =  ω_q2^⋆/ 2 E^⋆ 2𝒞_σ, JM(β , q^⋆)/(q^⋆2-k^⋆2+ i ϵ)^2   
    -   ω_q2^⋆/ 2 E^⋆ 2∂_q^⋆ 2𝒞_σ, JM(β , q^⋆)  /q^⋆2-k^⋆2+ i ϵ  
        + 𝒪 [(q^⋆ 2 - k^⋆ 2)^0  ]  , 
        =   ω_q2^⋆/ 2 E^⋆ 2∑_n=1^2     (- ∂_q^⋆ 2)^2-n 𝒞 _σ, JM(β , q^⋆)/(q^⋆2-k^⋆2+ i ϵ)^n     + 𝒪 [(q^⋆ 2 - k^⋆ 2)^0  ]  .

It follows that Eq. (<ref>) can be rewritten as

    G_σ (P,P,L) =  1/(q^⋆ )^ℓ_i + ℓ_fω_q2^⋆/ 2 E^⋆ 2∑_n=1^2    ∑_J M [(- ∂_q^⋆ 2)^2-n 𝒞 _σ, JM(β , q^⋆) ]   [1/L^3∑_𝐤ω_k2^⋆/ω_k2 - ∫d^3 k^⋆/(2π)^3]    √(4 π)  k^⋆ JY_J M(k̂^⋆)  /(q^⋆2-k^⋆2+ i ϵ)^n    .









To reach our final form we make two additional modifications. First we introduce a cutoff function on the sum-integral difference to enable effective numerical evaluation. Second, we re-express our results in terms of dimensionless quantities 𝐫^⋆=𝐤^⋆ L/(2π) and x=q^⋆ L/(2π). Then, shuffling around terms and introducing a new geometric function, we conclude

    G_σ (P,P,L) =  1/(q^⋆ )^ℓ_i + ℓ_fω_q2^⋆/ 2 E^⋆ 2∑_n=1^2 L^2n-3/(2 π)^2n∑_J M(2 π)^J/L^Jlim_α→ 0𝒵^(n)_J M(P,L,α)   (- ∂_q^⋆ 2)^2-n 𝒞 _σ, JM(β , q^⋆)       ,

where

    𝒵^(n)_J M(P,L,α) =     [ ∑_𝐫ω_k2^⋆/ω_k2 - ∫   d^3 r^⋆]    
      √(4 π)  r^⋆ JY_J M(r̂^⋆)     / (x^2-r^⋆2+ i ϵ)^n   e^- α (r^⋆2 - x^2)^n .

These two equations give the main result of this sub-section. In Appendix <ref> we give some details about the evaluation of 𝒵^(n)_J M(P,L,α). 
 In Fig. <ref> we plot 𝒵^(n)_J M(P,L,α) for P  = (2 π / L) [001], for various values of J and two different volumes. 

We give two final comments concerning the new kinematic function, 𝒵^(n). First we note an advantage of the decomposition over a single spherical harmonic performed in Eq. (<ref>). It is now  straightforward to use the symmetries of the finite-volume system to identify for which values of J M, 𝒵^(n)_JM will be nonzero. This is discussed in detail in Appendix <ref>.

Second, we remark that the cutoff function used here is designed with the property that the 𝒪(α) correction cancels the pole, and thus generates a smooth quantity within the sum-integral difference. If, for example, one were to instead use e^- α (r^⋆2 - x^2) for all n values, this would still be formally valid, but would lead to corrections of the form α/L^k, making it more difficult to estimate the α→ 0 limit. In fact these can be systematically subtracted and, as we shown in Appendix <ref>,
 this proves to be an efficient alternative approach for evaluating lim_α→ 0𝒵^(n)_J M.






 
We close this sub-section by returning to the specific case discussed above, σ = [μ=z; 10; 10] and P = (2 π/L) [00d_z]. Suppressing the arguments of 𝒵^(n), the final result for G^σ can be written as

    G^σ(P,P,L) =   - 1/ q^⋆ 2ω_q2^⋆/ 2 E^⋆ 21/ 4 π^2 L [       𝒵^(1)_00P^z/E^⋆ (   ω_q2^⋆  + q^⋆2/2 ω_q2^⋆ )+     2 π/L𝒵_10^(1)3√(3)/5E/E^⋆  +    (2 π)^2/L^2𝒵^(1)_202√(5)/5P^z/E^⋆1/2 ω _q2^⋆ ]  
    
    +
     1/ q^⋆ 2ω_q2^⋆/ 2 E^⋆ 2L/(2 π)^4[     𝒵^(2)_0 0P^z/E^⋆  q^⋆2ω_q2^⋆     +    2 π/L𝒵^(2)_1 03√(3)/5E/E^⋆  q^⋆ 2 +    (2 π)^2/L^2𝒵^(2)_2 02√(5)/5P^z/E^⋆ω _q2^⋆  +     (2 π)^3/L^3𝒵^(2)_3 06 √(7)/35E/E^⋆ ]     .

Note that, in the case of mass-degenerate particles, 𝒵^(1)_JM=0 for all odd J. If the particles are at rest in the finite-volume frame then all odd-J functions vanish as does J=2. This holds for both n=1 and n=2 for both degenerate and non-degenerate particles. [See again Appendix <ref>.]




In Fig. <ref> we plot the real and imaginary parts of G^μ_10; 10 (P,P,L) for P  = (2 π / L) [001]. 

 
 
 




 


 §.§ P_i ≠ P_f


We now turn our attention to the more challenging general case of P_i ≠ P_f. Note that this is realized if any of the four components of the 4-vectors differ, in particular also for P_i = P_f but E_i ≠ E_f.
As with P_i = P_f, here the evaluation of the sum is straightforward, while the integral is significantly more challenging. One of the major complications is that the two poles do not coincide in general as one varies the integration variable, k, but may overlap on a two-dimensional subspace for certain choices of external momenta. The contribution of this double-pole submanifold to the integral must be treated with care.


Though we have not found a fully analytic determination of the integral entering G^σ, we do have a recipe that gives the desired quantity accurately and with high efficiency. The approach is to rewrite the three-dimensional integral in terms of a d^4 k integral plus a second d^3 k integral with a smooth, singularity-free integrand. The four-dimensional integral, which carries all of the singularity structure, can then be reduced 
to a one-dimensional integral over a Feynman parameter. The second, smooth term can directly be evaluated using standard numerical integrators. 

To give the relevant expressions, we first introduce an extension of the cutoff function entering the definition of 𝒵_JM^(n) in the previous subsection

    H( α̅, k) ≡ e^- α̅(k^⋆ 2_i - q^⋆ 2_i) (k^⋆ 2_f - q^⋆ 2_f) = e^-   α   (r^⋆ 2_i - x^ 2_i) (r^⋆ 2_f - x^ 2_f)  ,

where α̅= L^4/(2 π)^4 α and r^⋆ 2_i(2π/L)^2=k^⋆ 2_i, x^ 2_i(2π/L)^2=q^⋆ 2_i, etc. 
We then write

    G_σ (P_f,P_i,L) = lim_α̅→ 0 [  𝒮_σ(α̅,P_f,P_i,L)
    -
     ℐ_σ(α̅,P_f,P_i)   ] ,

where

    𝒮_σ(α̅,P_f,P_i,L)
        ≡1/L^3∑_𝐤 
    H(α̅,k)  𝒴_ℓ_f m_f(k^⋆_f)  
     D(m,k) K^ω_μ_1⋯μ_n(m,k)   𝒴^*_ℓ_i m_i(k^⋆_i)   
      , 
    ℐ_σ(α̅,P_f,P_i)   
       ≡∫d^3 k/(2π)^3    H(α̅,k) 
     𝒴_ℓ_f m_f(k^⋆_f)  
     D(m,k) K^ω_μ_1⋯μ_n(m,k)   𝒴^*_ℓ_i m_i(k^⋆_i)   
      .

The sum can be evaluated directly as written, 
and thus we make no further modifications to 𝒮_σ. The remainder of this section is dedicated to ℐ_σ. 




  §.§.§ Separation into ℐ_𝒜; σ(  P_f, P_i) and ℐ_𝒩; σ( α̅, P_f, P_i)




As summarized above, our approach is to split the integral into a singular part that can be evaluated semi-analytically, denoted ℐ_𝒜; σ(  P_f, P_i), and a smooth remainder that is well-suited to numerical evaluation, denoted ℐ_𝒩; σ(α̅,  P_f, P_i). To proceed we define

    D_c(m,  k) ≡i/k^2-m_2^2+iϵ1/  (P_f-k)^2 - m_1^2 + i ϵ 1/  (P_i-k)^2 - m_1^2 + i ϵ ,

and then introduce D_r(m, k) via

    ∫d k^0/2 π D_c(m,k) ≡ D(m, k) + D_r(m, k)  .

In Appendix <ref> we give an explicit expression for D_r. This term mops up the contributions from the two P-dependent poles in D_c. Here the subscripts c and r stand for covariant and remainder, respectively. The idea is to perform the integral of D_c semi-analytically and that of D_r numerically, and then to take the difference.

In these relations we have neglected the possible factors of k_μ and the spherical harmonics multiplying D(m, k). To include these, we first introduce a tensor, M, defined such that

    M^ν_1 ⋯ν_N_ [μ_1 ⋯μ_n;  ℓ_f m_f ;     ℓ_i m_i ]  K^ω_ν_1⋯ν_N(m,k) ≡𝒴_ℓ_f m_f(k^⋆_f)   K^ω_μ_1⋯μ_n(m,k)   𝒴^*_ℓ_i m_i(k^⋆_i)     ,

where N = n+ℓ_f+ℓ_i. A more explicit definition, together with various examples, is given in Appendix <ref>. This simply amounts to recasting the factors of k^⋆ within the harmonics as boosted factors of k^μ.

To incorporate K^ω_ν_1 ⋯ν_N, note that D_r receives two contributions, one each from the poles at k^0 = E_f + ω_ P_f  k1 - i ϵ and k^0 = E_i + ω_ P_i  k1- i ϵ. (The third pole, at k^0 = ω_k2 - i ϵ, generates the term we are after, D(m,k).) We thus define

    D_rf(m, k)    ≡_E_f + ω_ P_f  k1d k^0/2 π D_c(m,k)  , 
    
    K^f_ν_1⋯ν_N(m, k)    ≡  k_ν_1⋯ k_ν_N|_k^0 = E_f + ω_ P_f  k1 ,

where the integral here is a closed clockwise contour encircling the pole indicated. The definitions with an i subscript are given by making the replacement f → i everywhere, and explicit expressions for D_ri(m, k) and D_rf(m, k) are given in App. <ref>. With these quantities in hand, Eq. (<ref>) generalizes to

    d k^0/2 π D_c(m,k) K_ν_1 ⋯ν_N(k) = D(m, k) K^ω_ν_1 ⋯ν_N(m, k) +𝒦_r; ν_1 ⋯ν_N(m, k)  .

Here  the integral on the left-hand side is a closed clockwise contour encircling the three poles below the real axis and we have also introduced

    𝒦_r; ν_1 ⋯ν_N(m, k) ≡  D_rf(m, k)K^f_ν_1 ⋯ν_N(m, k) + D_ri(m, k) K^i_ν_1 ⋯ν_N(m, k) .


The next step is to address the issue of ultraviolet convergence for these integrals. Equation (<ref>) is manifestly convergent, due to the inclusion of the cutoff function H(α̅ , k). 
But to reach an integral that can be evaluated analytically it is convenient to introduce a second form of ultraviolet regularization. We explain the approach first for the special case of two indices, N=2. Here the integral has a logarithmic divergence that can be removed by taking Eq. (<ref>) and subtracting from this the same equation defined with m_1,2→Λ

    ∫d k^0/2 π [ D_c(m,k) - D_c(Λ,k)  ] K_ν_1  ν_2(k) =  [ D(m, k) K^ω_ν_1  ν_2(m, k) -   D(Λ, k) K^ω_ν_1  ν_2(Λ, k)  ]  
    
    +   [ 𝒦_r; ν_1  ν_2(m, k)  -  𝒦_r; ν_1  ν_2(Λ, k)    ]  .

The regularization scale, Λ, is chosen so that the integrands that depend on it are smooth functions of k with no need for an i ϵ prescription.  This holds for any Λ satisfying 2 Λ >   max[E^⋆_f,E^⋆_i] though in practice it is useful to take the cutoff well above this minimum value.  



On the left-hand side of Eq. (<ref>) we have used that, for N=2, the k^0 integral can be extended to the entire real axis, with a vanishing arc at negative complex infinity. We additionally note that, as a result of the subtraction, the left-hand side and also both square bracketed terms on the right-hand side vanish as 1/|k|^5 in the limit |k|→∞. These thus give convergent integrals with respect to d^3 k.

This approach can be extended to any number of k_μ factors, simply by forming more complicated linear combinations to cancel all divergent powers

    ∫d k^0/2 π∑_j=0^n_j c_j D_c(Λ_j,k)   K_ν_1 ⋯ν_N(k) =   ∑_j=0^n_j c_j D(Λ_j, k) K^ω_ν_1 ⋯ν_N(Λ_j, k) +   ∑_j=0^n_j c_j 𝒦_r; ν_1 ⋯ν_N(Λ_j, k)     ,

where we have introduced c_0=1 and Λ_0= {m_1, m_2}. As above, for j>0 we take Λ_j such that the corresponding integrands are smooth functions of k (Λ_j>0 >  max[E^⋆_f,E^⋆_i]). In all cases, the linear combinations are constructed such that (i) the k^0 integral extends to the entire real axis with a vanishing contribution from the arc at negative complex infinity, and (ii) the left-hand side and each of the two sums on the right-hand side give convergent integrals with respect to d^3 k. In the following we sometimes refer to this as a Pauli-Villars-like regulator. We give a general algorithm for forming these linear combinations in Appendix <ref>.

The final step is to multiply both sides of Eq. (<ref>) by the cutoff function, and solve for the desired integral, defined in Eq. (<ref>). We deduce

    ℐ_σ(α̅,P_f,P_i)   = ℐ_𝒜; σ(  P_f, P_i)  + ℐ_𝒩;σ(α̅, P_f, P_i)   ,

where

    ℐ_𝒜; σ(  P_f, P_i)     ≡  M^ν_1 ⋯ν_N_σ∫d^4 k/(2 π)^4 K_ν_1 ⋯ν_N(k)
     ∑_j=0^n_j c_j D_c(Λ_j,k)   
     , 
    
    ℐ_𝒩; σ(α̅, P_f, P_i)    ≡  M^ν_1 ⋯ν_N_σ [   ∫d^3 k/(2π)^3  
     [H(α̅, k) - 1]
      ∑_j=0^n_j c_j  [ D(Λ_j, k) K^ω_ν_1 ⋯ν_N(Λ_j, k) +𝒦_r; ν_1 ⋯ν_N
    (Λ_j, k)  ]  
     
    
       
     -  ∫d^3 k/(2π)^3   H(α̅, k)  ∑_j=1^n_j  c_j D(Λ_j, k) K^ω_ν_1 ⋯ν_N(Λ_j, k)
     
    
       
     -  ∫d^3 k/(2π)^3   H(α̅, k)   ∑_j=0^n_j c_j   𝒦_r; ν_1 ⋯ν_N(Λ_j, k)    
      ]  .

Equations (<ref>)-(<ref>) are the main results of this subsection.

We emphasize here that the integrands in the definition of ℐ_𝒩 are smooth for all real k, and the integrals are ultraviolet convergent. For the second and third terms this follows from the fact that H(α̅, k) decays exponentially, together with the observation that 𝒦_r is a smooth function, as is D(Λ_j, k) for j>0. For the first term, smoothness is guaranteed because [H(α̅, k) - 1] vanishes at the pole and ultraviolet convergence follows from the careful construction of the linear combination. As the integrands are smooth and the integrals are convergent, ℐ_𝒩; σ(α̅, P_f, P_i) is well suited to numerical evaluation. 

It is instructive to consider a few specific examples of this construction, beginning with n=0, ℓ_i=ℓ_f=0, in which no factors of k_μ appear in the numerator of ℐ_σ. In this case the integrals are directly convergent, without any need for additional subtraction terms (i.e. the sum over j in Eq. (<ref>) reduces to the j=0 term). Equations (<ref>) and (<ref>) then reduce to 

    ℐ_𝒜(  P_f, P_i)     ≡∫d^4 k/(2 π)^4   D_c(m,k)    , 
    ℐ_𝒩(α̅, P_f, P_i)    ≡∫d^3 k/(2π)^3 [ D(m, k)   +  D_rf(m, k)   + D_ri(m, k)  ]   [H(α̅, k) - 1]  
        -  ∫d^3 k/(2π)^3  H(α̅, k)    [ D_rf(m, k)   + D_ri(m, k)  ]    .

This no-index version of ℐ_𝒩 is plotted in the left panel of Fig. <ref>. 

 We emphasize here that ℐ_𝒩, and thus also the original integral, diverges in the α̅→ 0 limit. This is because in the original integral the covariant propagators are evaluated at on-shell k (k^2 = m^2), so the propagators scale as 1/|k| and the integrand as d^3 k / |k|^3. In other words the convergence of ℐ_𝒜 is always better than that of the original integral by two powers of k, resulting from the off-shell integration of k^μ. For fixing the subtraction scheme in Eqs. (<ref>) and (<ref>), it is only necessary that ℐ_𝒜 be rendered finite. [See also Appendix <ref>.]









We close with one final example: σ = [μ =0; 10; 10], corresponding to a factor of 𝒴_10 k_0 𝒴^*_10 in the numerator. This leads to an ℐ_𝒜;σ integral with an integrand scaling as d^4k k^2/k^6, i.e. diverging as logΛ. Performing a single subtraction of the same integral with m →Λ is therefore sufficient to render the result finite. In fact, to improve the numerical evaluation of ℐ_𝒩;σ, and to test our general method, here we choose to make two subtractions. As explained in Appendix <ref>, one possible choice is to add an integral evaluated at Λ = 3m (with coefficient -35/27) and a second at Λ = 6m (with coefficient 8/27). Implementing this in Eqs. (<ref>) and (<ref>) leads to convergent forms of ℐ_𝒜;σ and ℐ_𝒩;σ respectively, with integrands scaling as d k/k^5. ℐ_𝒩;σ in this scheme is plotted in the right panel of Fig. <ref>. 

 

As we include additional factors of k_μ the expressions complicate, first because we need additional terms in the sum over j (to maintain convergent integrals) and second because the numerical integrals depend on multiple vector components. However, we find that no conceptual issues arise and the task amounts to coding Eq. (<ref>) with an efficient numerical integrator. We give some details on our approach in Appendix <ref>, but consider ℐ_𝒩;σ(α̅, P_f, P_i) as a numerically known function for the remainder of the main text. 

Thus it remains only to evaluate ℐ_𝒜; ν_1 ⋯ν_N(  P_f, P_i), to which we now turn.



  §.§.§ Evaluating ℐ_𝒜; σ(  P_f, P_i)


We first use the tensor M^ν_1 ⋯ν_N_σ, introduced in Eq. (<ref>) above, to define a version of ℐ_𝒜 with no spherical-harmonic indices

    ℐ_𝒜; σ(  P_f, P_i) ≡ M^ν_1 ⋯ν_N_σℐ_𝒜; ν_1 ⋯ν_N (P_f,P_i)  .

The integrals on the right-hand side can then be written as

    ℐ_𝒜; ν_1 ⋯ν_N(  P_f, P_i)   ≡lim_δ→ 0∑_j=0^n_j c_j ℐ_𝒜; ν_1 ⋯ν_N(  P_f, P_i,  Λ_j, δ)   ,

where

    ℐ_𝒜; ν_1 ⋯ν_N(  P_f, P_i, m, δ)   ≡∫d^4 - δ k/(2 π)^4 - δi/k^2-m_2^2+iϵ1/ (P_f-k)^2 - m_1^2 + i ϵ 1/  (P_i-k)^2 - m_1^2 + i ϵ  k_ν_1⋯ k_ν_N .

Here we have used the fact that the sum over j gives a convergent integral, and is thus equal to the δ→ 0 limit of the integral in 4 - δ dimensions. Then, at fixed delta, one can exchange the orders of summation and integration, leading to Eq. (<ref>).

To evaluate the integral in Eq. (<ref>), we first perform the standard Wick rotation on the k^0 integration contour (counter-clockwise to the imaginary axis), and similarly rotate the time components P_i,0 and P_f,0. We then make the variable redefinitions




    k_0 ≡ ik_E,0, P_f,0≡ i P_E,f,0  , P_i,0 = i P_E,i,0 ,

to reach

    ℐ_𝒜; ν_1 ⋯ν_N( P_f, P_i, m, δ) 
       =
    ξ_ν_1⋯ν_N∫d^4 - δ k_E/(2π)^4 - δ1/k^2_E+m_2^21/  (P_E,f-k_E)^2 + m_1^2    1
     /  (P_E,i-k_E)^2 + m_1^2    k_E,ν_1 ⋯  k_E,ν_N , 
    
       =
    ξ_ν_1⋯ν_N∂/i∂χ_E^ν_1⋯∂/i∂χ_E^ν_Nℐ^χ(P_f, P_i, m,δ)
    |_χ=0 ,

where ξ_ν_1⋯ν_N≡(i)^δ_ν_1,0+⋯+ δ_ν_N,0. Note here that the indices are not contracted between ξ_ν_1⋯ν_N and the momenta but rather are held fixed on both sides of the equation.

In the second step we have introduced the generating functional

    ℐ^χ(P_f, P_i, m, δ)
       ≡∫d^4 - δ k_E/(2π)^4 - δe^iχ_E · k_E/k^2_E+m_2^21/  (P_E,f-k_E)^2 + m_1^2   1/  (P_E,i-k_E)^2 + m_1^2    .

As we show in Appendix <ref>, this reduces to


    ℐ^χ(P_f, P_i, m, δ)
     = 
     ∫_0^1dx∫_0^1-xdy
      e^- iχ· (x P_f+y P_i)
       ∑_n=0^∞ (χ^2 )^n/ 4^n n!Γ(1-n + δ/2)/(4π)^2 - δ/2  M(m,x,y)^2n - 2 - δ ,

 where
  
    M(m,x,y)^2     =  (1-x-y)m_2^2+ (x+y)m_1^2-x(1-x) s_f-y(1-y) s_i + x y  (Q^2 + s_f + s_i)-iϵ .

In these results we have analytically continued back to real P_i,0 and P_f,0 and expressed all quantities in terms of the 4-vectors P_f amd P_i as well as the Lorentz invariants Q^2= - (P_i-P_f)^2, s_i=P_i^2, and s_f=P_f^2.
The corresponding analytic continuation of Eq. (<ref>) is given by

    ℐ_𝒜^ν_1 ⋯ν_N( P_f, P_i, m, δ)      =
    ∂/-i∂χ_ν_1⋯∂/-i∂χ_ν_Nℐ^χ(P_f, P_i, m,δ)
    |_χ=0 .

Taken together, Eqs. (<ref>)-(<ref>) 
give the main result of this subsection.  


As above, it is instructive to consider ℐ_𝒜 for N=0, i.e. with no factors of k_μ in the numerator. In this case Eq. (<ref>) is evaluated with no χ-derivatives and gives
  
    ℐ_𝒜( P_f, P_i, m, δ)  
    
    
         =   Γ(1 + δ/2)/(4π)^2 - δ/2∫_0^1dx∫_0^1-xdy  1/M(m,x,y)^ 2 + δ .

 
 As noted at the end of the previous subsection, this integral is convergent in the δ→ 0 limit, and simplifies to 
  
    ℐ_𝒜( P_f, P_i, m, 0) 
    
       =
     1/(4π)^21/s_i∫_0^1dx∫_0^1-xdy
     
    1/[y-y_+(x)][y-y_-(x)] ,

where we have substituted 
 
    M(m,x,y)^2 = s_i  [y-y_+(x)][y-y_-(x)]  ,

 with

    y_±(x)    ≡1/2( A ±√(A^2 + B + iϵ))  , 
    
    A    ≡ 1+m_2^2-m_1^2+x(q^2- s_f- s_i)/s_i , 
    
     B    ≡ - 4m^2_2 -x (m^2_2-m^2_1) -  x(1-x) s_f/s_i .



The final analytic step is to evaluate the integral with respect to y.

We do so via the identity

    ∫_0^1-xdy/y-(f(x) ± i ϵ) =
     log|1 - x -f(x)/f(x)| + i arctan( 1 - x-Re f(x) /Im f(x)  ±ϵ )  +  i arctan( Re f(x) /Imf(x)  ±ϵ )   ≡ℒ_±ϵ[f(x)]
     
     ,

where in the second equality we have introduced a short-hand for the result of the integral. Note that, as long as Imf(x) is nonzero, we can safely set ϵ = 0 in these expressions.  In the case that the imaginary part does vanish, then we use the relation lim_ϵ→ 0arctan(a/ϵ) = sign[a] π/2.



Applying this to Eq. (<ref>) we deduce
  
    ℐ_𝒜( P_f, P_i, m, 0)  =
    ∫_0^1dx   F^(1)(x)  ,

where F^(1)(x) [also given in Eq. (<ref>)] is
 
    F^(1)(x) ≡1/(4π)^2 s_iℒ_+ ϵ[y_+(x)] - ℒ_- ϵ[y_-(x)]/y_+(x)-y_-(x) .
 
In this case Im y_+(x) = - Im y_-(x) = Im√(A^2 + B)/2. Thus the substitution lim_ϵ→ 0arctan(a/ϵ) = sign[a] π/2 is only needed when the argument of the square root is positive. As we explain in detail in Appendix <ref>, evaluating the remaining integral over x reveals that ℐ_𝒜(P_f,P_i) has triangle singularities that arise whenever P_f and P_i take on values for which all three particles in the triangle of Fig. <ref>(b) can go on shell.





More precisely, we show in Appendix <ref> that the singularity locations are governed by the discriminant of the polynomial A^2 + B ≡ a x^2 + b x + c, given by

    X (s_f,s_i,Q^2) ≡ b^2 - 4 a c = m_1^2 ((s_f-s_i)^2+Q^2(2 m_2^2 -m_1^2 + s_f+s_i))-Q^2 (m_2^4-m_2^2 (Q^2+s_f+s_i)+s_f s_i)   .

Critical kinematics are realized whenever X (s_f,s_i,Q^2) = 0, so that A^2 + B = (x-x_c)^2, and in addition, x_c and y_c ≡ y_+(x_c) fall in the integrated range. It can be shown that these conditions are equivalent to the ones found using Landau's
singularity classification <cit.>.

At values of P_f and P_i satisfying these conditions, the real part of ℐ_𝒜 has a step-function discontinuity of height

    Disc(ℐ_𝒜)
       = 1/16√((P_i· P_f)^2 -   s_f s_i) ,

and the imaginary part shows a logarithmic divergence.
In Fig. <ref> we illustrate how these singularities form in the ϵ→ 0 limit of the i ϵ pole prescription. 
In Fig. <ref> we show the singularity structure as a function of E_f^⋆ for various fixed values of E_i^⋆. In particular one sees that, for sub-thresold E_i^⋆, ℐ_𝒜 is a smooth function away from the two-particle production threshold. As E_i^⋆ approaches 2m a step forms in Reℐ_𝒜 and Imℐ_𝒜 develops a log divergence. Then, as E_i^⋆ is further increased, the location of the singularity in E_f^⋆ moves towards and eventually collides with the two-particle threshold. 

To complete this subsection we would like to comment on the behavior of this singularity for some special set of kinematics. First, in the case of identical initial and final 3-momenta, i.e. P_i=P_f, ℐ_𝒜 does not have any other singularities apart from those arising at threshold, and therefore the G-function will not exhibit a triangle singularity. This is consistent with our analysis of the P_f=P_i case, and with the numerical example shown in Fig. <ref>. In other words, given that all the external momenta in the triangle diagram are time-like, the condition of all three internal propagators to be on shell cannot be realized. A second example is the special case of m_1 = m_2 and s_i = s_f = s. Solving X(s,s,Q^2)=0 in this simple case leads to a singular manifold given by

    Q^2 = -4s (1-s/4 m^2)  .

In this case it is easier to visualize that this condition is equivalent to that of all three intermediate particles in the triangle diagram, Fig <ref>(b), going on-shell. (See also Refs. <cit.>.)










 
This concludes our discussion of ℐ_𝒜 within the main text.  In Appendix <ref> we extend the results here by explicitly evaluating ℐ_𝒜; ν_1, ℐ_𝒜; ν_1  ν_2 and ℐ_𝒜; ν_1  ν_3  ν_3. For these integrals we find that the above-threshold discontinuities persist, but are milder when factors of k_μ appear in the numerator.
 
 
 

  §.§.§ Examples of G^σ(P_f,P_i,L)











  
 
 Having discussed the integral entering G^σ(P_f,P_i,L) in great detail we are now ready to put everything together and evaluate the complete function. We do so for two different examples of external kinematics. First, in Fig. <ref>, we consider the case of P_i = P_f = (2 π/L)[001] with σ = [μ=0; 10; 10]. Then, in Fig. <ref> we take P_i = [000] and P_f=(2 π/L)[001] with σ = [00;00], i.e. with no numerator factors. In both cases we fix mL=6 and plot Re G^σ for all values of E_i^⋆ and E_f^⋆ in the region of interest. 
 
 As explained in the figure captions, each example illustrates important issues and features that arise. In Fig. <ref> we consider various diagonal slices of the E_i^⋆, E_f^⋆ plane. We find that the result for E_i^⋆ = E_f^⋆ is in perfect agreement with the P_i = P_f result determined by combining the various 𝒵_JM functions. This provides a strong check on the two different methods.  
 Figure <ref> also illustrates that double poles arise along the E_i^⋆ = E_f^⋆ line, but these split to single poles as the slice is rotated away from this singular choice.
 
 In addition, Fig. <ref> illustrates the results of using the Pauli-Villars-like regulator to separate the integral into ℐ_𝒜 and ℐ_𝒩. As discussed towards the end of Sec. <ref>, the original integral contains a factor scaling as 𝒴_10 k^μ𝒴_10∼|k|^3 in the numerator leading to a logΛ divergence in ℐ_𝒜. Following Appendix <ref> we handle this by evaluating Eqs. (<ref>) and (<ref>) with {Λ_1, Λ_2} = {3m, 6m} and {c_1, c_2} = {-35/27,8/27}. This removes not only the divergence but also a dk/k^3 term in the integrand to further optimize the numerical convergence of ℐ_𝒩. We already gave the result for ℐ_𝒩 in this prescription in the right panel of Fig. <ref>. Here the result is combined with the sum and ℐ_𝒜 to reach G^σ. With all building blocks summed together, the dependence on the Pauli-Villars parameters, Λ_i, cancels (between ℐ_𝒜 and ℐ_𝒩) as does the dependence on the smooth cutoff parameter, α̅, (between ℐ_𝒩 and the sum). 
 
Turning now to Fig. <ref>, this result displays two additional features of G^σ. First we see that, when P_i and P_f differ, the non-interacting two-particle poles appear in different locations for E_i^⋆ and E_f^⋆. Interactions shift the finite-volume energies away from these singularities so that G^σ, like the Lellouch-Lüscher factors, will generally be evaluated away from the divergent locations. However, as with all finite-volume kinematic functions, this implicit knowledge of the non-interacting spectrum is a key ingredient in the all-orders correction of the scattering-state volume effects. 

 



Second, we see the appearance of triangle singularities inherited through ℐ_𝒜. Such features are simply part of the correct definition of G^σ. Indeed, because the singularity structure is directly induced by the infinite-volume diagram of Fig. <ref>(b), it also appears within the infinite-volume 2 + 𝒥→2 transition amplitude itself. The steps, cusps and log divergences of G^σ are present in both 𝒲 and 𝒲_df for exactly the same kinematics. Thus, understanding the features is crucial to extracting and interpreting the infinite-volume observables that we are after. We discuss the role of discontinuities within the transition amplitude in more detail in future work. 















 


§ CONCLUSION 


In this work we have presented a modified version of the finite-volume formalism for studying 2 + 𝒥→2 transition amplitudes. This is closely related to the approach of Ref. <cit.>, but differs in that all infinite-volume quantities are Lorentz covariant and the 1 + 𝒥→1 matrix elements have been reformulated in terms of standard form factors. As explained in Sec. <ref> and Appendix <ref>, the new result is reached by making minor adjustments to the derivation presented in Ref. <cit.>. For example, in that work finite-volume effects are expressed as sums over poles of the form 1/[2 ω_k (k^0 - ω_k)], and here the same effects are expressed via invariant poles, 1/(k^2 - m^2). 

These changes lead to a modified form of the finite-volume function, denoted G, with an added benefit that the new form is easier to evaluate numerically. As described in Sec. <ref>, the Lorentz covariant structure allow us to write the integral appearing in G in terms of Feynman parameters. This reduction is also crucial to revealing the analytic structure of G, including the triangle singularities described in great detail Sec. <ref> and Appendix <ref>, and illustrated in Figs. <ref>, <ref> and <ref>. We also recall that the form of G presented in Ref. <cit.> carries four sets of spherical harmonic indices, resulting from a cumbersome description of the 1 + 𝒥→1 matrix elements. By contrast, our improved expression carries only the angular momentum indices of the external states, together with Lorentz indices to describe the current insertion.




To avoid proliferation of flavor and channel indices, in this work we restricted attention to kinematics for which a single channel of two scalar particles is open.  
Accommodating multiple channels is straightforward given the results of Ref. <cit.>. 



Incorporating particles with spin has yet to be considered for these types of observables and is the subject of future work. 

This final generalization is of great importance given the phenomenological interest in two-nucleon matrix elements <cit.>. From

 our previous experience with spinning particles <cit.>, we expect that the extension will be relatively straightforward. 

Looking to less trivial extensions, it would be of great interest to extended these ideas in order to develop an approach for extracting non-local matrix elements of two-particle systems. This would make it possible to extract distribution functions of resonant states, following the methods of Refs. <cit.>, and would open the possibility for lattice QCD calculations of two-body contributions to double-beta decays <cit.>. Matrix elements of non-local operators suffer from different types of finite-volume artifacts. These depend crucially on whether the operators are displaced in Euclidean time, as in Ref. <cit.>, or in a spatial direction, as considered for example in Ref. <cit.>. 
Finally, extensions of this work to energies for which three or more particles can go on-shell should be feasible in the future, especially given the recent progress in understanding the finite-volume spectrum of three-particle states <cit.>. 

Returning to the present formalism, several open questions persist that we plan to address in future work. For example,  is defined using the partial-wave basis, while the form factors of resonances or bound states are more naturally described using Lorentz decomposition. It is always possible to relate the partial-wave and Lorentz bases. However, due to the reduction of rotational symmetry in the finite-volume, we do not expect a one-to-one correspondence between the finite-volume matrix elements and the different Lorentz components of . For example, in the ρ channel with nonzero spatial momenta in the finite-volume frame, the different helicity components mix to different finite-volume irreps. This means that the components are sampled by different finite-volume quantization conditions and thus at different energies. As a result, just as is done in the analysis of coupled-channel scattering <cit.>, it will be necessary to perform global fits of the matrix elements using Eq. (<ref>). This requires a detailed understanding of the analytic structure of these amplitudes in which triangle singularities play a crucial role.
















§ ACKNOWLEDGEMENTS


RAB acknowledges support from U.S. Department of Energy contract DE-AC05-06OR23177, under which Jefferson Science Associates, LLC, manages and operates Jefferson Lab, and the U.S. Department of Energy Early Career award contract DE-SC0019229. AB acknowledges support from U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Award No. DE-SC0010300.
The authors would like to thank A. Jackura, as well as J. Dudek, R. Edwards, D. Wilson, and the rest of the Hadron Spectrum Collaboration for useful discussions.


























§ DERIVATION OF THE COVARIANT FORMALISM








In this appendix we describe how the derivation of Ref. <cit.> must be modified to give the covariant version of the formalism, presented in Sec. <ref> above. As in our previous derivation, we restrict attention to a finite cubic spatial volume, implemented by requiring all fields to have periodicity L in the three spatial directions. 

Within this set-up we introduce a finite-volume three-point function, C_L^2 → 2(P_f,P_i), defined as the sum over all possible diagrams connecting the initial and final states to the inserted current. See also Fig. <ref>(a). For E_i^⋆, E_f^⋆ below the next multi-hadron threshold (such that only a single two-particle channel can propagate), all volume effects scaling as a power of 1/L are captured by the skeleton expansion shown in the figure. Here the label V within the loops stands for volume and denotes that the diagrams are defined with the spatial momenta summed over the discrete set allowed by the periodic boundary conditions, k = 2 πn/L with n a 3-vector of integers. The corresponding diagrams in an infinite volume are given by replacing these sums with integrals and are represented by a loop with no label.


The power-like volume effects of C_L^2 → 2(P_f,P_i) are encoded in the skeleton expansion of Fig. <ref>(a), built from fully dressed hadron propagators (indicated by the simple black lines) and two-particle Bethe-Salpeter kernels (indicated by open circles). The vertices with a current insertion are given by the same diagrammatic set defining propagators and kernels, but with the current attached at all possible locations. In the kinematic window of interest, the difference between the finite- and infinite-volume definitions of propagators and kernels as e^- m L. 




 
The set-up here is identical to that of Ref. <cit.>. Indeed the only modifications required are in the evaluation of the two-particle loop, in which the current couples to one of the two particles. The relevant diagram is shown in Fig. <ref>(b). To explain how the analysis is altered, we begin by recalling the finite-volume reside of this diagram, given in Eq. (27) of Ref. <cit.>

    𝒢_L,μ_1 ⋯μ_n   ≡ [1/L^3∑_𝐤 - ∫d^3 k/(2π)^3 ]
    1/2ω_k2   iℒ(P_f, k)   Δ(P_f - k)  ^off_μ_1 ⋯μ_n (P_f-k,P_i-k)  Δ(P_i - k)   iℛ^†(P_i, k) |_k^0=ω_k2 .




Here ℒ(P_f, k) and ℛ^†(P_f, k) are generic endcap functions to be replaced with Bethe-Salpeter kernels or the overlap to the interpolators in the final derivation. 
In contrast to Ref. <cit.>, we define 𝒢_L,μ_1 ⋯μ_n with these endcaps accompanied by factors of i. The difference arises because we formulate the derivation here with Minkowski momenta (in contrast to the Euclidean conventions of the previous publication). In this set-up the i factors give a more natural extension to the Bethe-Salpeter kernels, multiplied by this factor due to the weight, e^iS, in the quantum path integral. Our second notational modification is to represent the 1 + 𝒥→1 insertion with a set of Lorentz indices, and to make explicit that the quantity is off-shell, i.e. that (P_f - k)^2, (P_i - k)^2 ≠ m_1^2. [See also Eq. (<ref>) above.] Our third and final alteration is to restrict attention to a single channel, thus removing the a and b indices from Eq. (27) or Ref. <cit.>.

We now express ℒ(P_f, k), ℛ^†(P_f, k) and ^off_μ_1 ⋯μ_n (P_f-k,P_i-k) in terms of their on-shell counterparts, plus corrections. For the endcap functions this is done exactly as in Ref. <cit.>, by first defining

    ℒ(P_f,  q^⋆_f  k̂^⋆_f)=∑_ℓ_f m_f √(4 π) Y_ℓ_f m_f (k̂^⋆_f)ℒ_ℓ_f m_f (P_f),
         ℛ^†(P_i,  q^⋆_i  k̂^⋆_i)
    =∑_ℓ_i m_i √(4 π) Y^*_ℓ_i m_i (k̂^⋆_i)ℛ^†_ℓ_i m_i(P_i)  ,

and then recombining the components with the 𝒴_ℓ m harmonics defined in Eq. (<ref>)

    ℒ_on(P_f, k^⋆_f ) ≡∑_ℓ_f m_f 𝒴_ℓ_f m_f ( k^⋆_f)ℒ_ℓ_f m_f (P_f)  ,
         ℛ_on^†(P_i, k^⋆_i )
    ≡∑_ℓ_i m_i 𝒴^*_ℓ_i m_i ( k^⋆_i)ℛ^†_ℓ m(P_i)  .

With these in hand we introduce the δ operator as follows

    ℒ(P_f,k)|_k^0=ω_k2    = ℒ_on(P_f, k^⋆_f ) + [ℒδ](P_f, k)  , 
    ℛ^†(P_i,k) |_k^0=ω_k2    = ℛ^†_on(P_i, k^⋆_i ) + [δℛ^†](P_i,k)  .

These results match Eqs. (41) and (42) of Ref. <cit.>, up to the minor notational differences discussed above. A key property that we will use below is that [ℒδ](P_f, k) vanishes like (P_f - k)^2 - m^2_1 in the on-shell limit, and [δℛ^†](P_i,k) vanishes like (P_i - k)^2 - m^2_1.
  
We now imitate this separation with the 1 + 𝒥→1 insertion, and in doing so introduce the first major difference as compared to our earlier work. Beginning with Eq. (<ref>) of the main text, above, we introduce the shorthand k_f ≡ P_f-k and k_i ≡ P_i-k to write

    ^off_μ_1 ⋯μ_n (k_f,k_i)     = ∑_j K^(j)_μ_1 ⋯μ_n(k , P_f,P_i)  f^(j)(Q^2, k_f^2, k_i^2)  .

We now follow the approach introduced in the main text by only projecting the scalar form factors, f^(j), to their on-shell values. To understand the idea note that

    f^(j)(Q^2, k_f^2, k_i^2)=  f^(j) [ Q^2;   P_f^2 + k^2 - 2 E^⋆_f ω^⋆ f_k2     ;   P_i^2 + k^2 - 2 E^⋆_i ω^⋆ i_k2   ]   ,

where ω^⋆ f_k2 is the temporal component of k^μ⋆ f, given by boosting k^μ = (ω_k2, k) to the P_f = 0 frame, and ω^⋆ i_k2 is the f → i analog. 

The key point is that, when expressed in this way, f^(j) has no dependence on k̂_f^⋆ or k̂_i^⋆. Thus there is no need to decompose in harmonics, nor to include barrier factors. The on-shell projection is simply 

    f^(j)(Q^2 )    ≡  f^(j)(Q^2, m_1^2, m_1^2)    = f^(j) [ ⋯   ] |_k^⋆_f = q^⋆_f, k^⋆_i = q^⋆_i , 
    
      f^(j)_off,on ( Q^2 , k_f^2 )    ≡ f^(j)(Q^2,k_f^2, m_1^2) =  f^(j) [ ⋯   ] |_  k^⋆_i = q^⋆_i , 
    
       f^(j)_on,off ( Q^2 , k_i^2 )    ≡ f^(j)(Q^2,m_1^2,k_i^2) =  f^(j) [ ⋯   ] |_k^⋆_f = q^⋆_f .

To avoid clutter we have suppressed the arguments on the right-hand side, identical to those of Eq. (<ref>). We stress here that the on-shell projections are subtle in that k is used to define two separate variables k^⋆_f and k^⋆_i. The separation is unambiguously given by whether the original k appears in k_f = P_f-k or k_i = P_i-k. With the induced k_f^⋆ and k_i^⋆ dependence, it is possible to separately project the initial and final states on shell via k^⋆_i → q^⋆_i and k^⋆_f → q^⋆_f respectively.

We next form linear combinations of the on- and partially-off-shell form factors in direct analogy to Eq. (37) of Ref. <cit.> to reach

    f^(j)(Q^2,k_f^2,k_i^2)
    =f^(j)(Q^2)
    +
    δ f^(j)(Q^2)
    +
     f^(j)(Q^2)δ
    +
    δ f^(j)(Q^2)δ .

The individual terms on the right-hand side are defined by analogy to Eqs. (38)-(40) of Ref. <cit.> and have the explicit form

    δ f^(j)(Q^2)
       ≡  f^(j)_off,on ( Q^2 , k_f^2 ) - f^(j)(Q^2)  ,  
    
     f^(j)(Q^2)δ   ≡     f^(j)_on,off( Q^2 , k_i^2 )  - f^(j)(Q^2)  ,  
    δ f^(j)(Q^2)δ   ≡ 
    f^(j)(Q^2, k_f^2,k_i^2) + f^(j)(Q^2)
    -  f^(j)_off,on ( Q^2 , k_f^2 )  
    -  f^(j)_on,off( Q^2 , k_i^2 )  - f^(j)(Q^2)  .

Note that functions with a δ on the left side (right side) also depend on k_f^2 (k_i^2), but we keep this dependence implicit to avoid clutter. As in Ref. <cit.>, the key point here is that δ appearing on either side of the function indicates that the latter vanishes in the on-shell limit, scaling as k_f^2 - m_1^2 (for δ on the left) and as k_i^2 - m_1^2 (for delta on the right). 
 
 It now remains only to separate the propagators by defining
 
    Δ(k^2) ≡𝒟(k^2) + 𝒮(k^2)  ,

 where
 
    𝒟(k^2) ≡i/k^2 - m_1^2 + i ϵ ,

 is the non-interacting covariant propagator. Our conventions are such that Δ(k^2) has unit residue at the single particle pole, implying 𝒮(k^2) is smooth and finite near k^2 = m_1^2. This form of separation, in which 𝒟(k^2) remains Lorentz invariant, is the second key distinction relative to our earlier formalism.
 
At this stage we have separated the endcaps [Eqs. (<ref>), (<ref>)] the single-particle form factor [Eq. (<ref>)] and the covariant propagator [Eq. (<ref>)] into on-shell terms plus corrections. Substituting these four identities into Eq. (<ref>) then gives the analog of Eq. (43) of Ref. <cit.>. The final step is to rearrange terms according to singularity structure. Of the 64 terms (reach by multiplying 4 binomials as well as the form-term separation of f^(j)) all but 17 are smooth at both poles, and thus give only exponentially suppressed contributions to the sum-integral difference. Those that are singular break into three classes, 8 with only the 1/(k_f^2 + m^2 - i ϵ) singularity, 8 with 1/(k_i^2 + m^2 - i ϵ), and a single maximally singular term that leads to the appearance of G^σ(P_f,P_i,L).

To give an explicit form for the single-pole terms we need to introduce one final piece of notation. We define the operator δ_df via

    [  [ℒ(P_f, k) Δ(k_f^2) _μ_1 ⋯μ_n(k_f, k_i)  ] δ_df ] ≡∑_jK^(j)_μ_1 ⋯μ_n(m,k,P_f,P_i) 
    × [ℒ(P_f,k)     Δ(k_f^2)  f^(j)( Q^2 , k_f^2,k_i^2 ) 
      -   ℒ_on(P_f, k^⋆_f ) 𝒟 (k_f^2)  f^(j)_on,off( Q^2 ,k_i^2 ) ]  ,

 and similarly for δ_df acting on the left.
 
 We can make use of this to rewrite Eq. (<ref>) as 
 
    𝒢_L,μ_1 ⋯μ_n   = 
     -∑_j
     f^(j)(Q^2) 
     i ℒ_ℓ_f m_f  (   [1/L^3∑_𝐤 - ∫d^3 k/(2π)^3]      
    𝒴_ℓ_f m_f(k^⋆_f) 
    D(m,k)    
     K^(j)_μ_1 ⋯μ_n(m,k,P_f,P_i) 
    𝒴^*_ℓ_i m_i(k^⋆_i)     )
      iℛ^†_ℓ_i m_i 
    
     +      [  [i ℒ(P_f, k) Δ(k_f^2) _μ_1 ⋯μ_n(k_f, k_i)  ξ^-1 ] δ_df ]    
       iF(P_i,L)    i ℛ^†
     +
    i ℒ
       iF(P_f,L)   
      [ δ_df [  ξ^-1_μ_1 ⋯μ_n(k_f, k_i) Δ(k_i^2) i ℛ^†(P_i, k)   ]   ]   ,

where D(m,k) is defined in Eq. (<ref>).
On the second line the square bracketed quantities have been decomposed in spherical harmonics and carry implicit indices that are contracted with F(P,L), defined in Eq. (<ref>). We have included the inverted symmetry factor, ξ^-1 = 2 for identical particles, to compensate the factor in the definition of F. Note that this arises naturally from the fact that the current couples to each of the two particles when these are identical. This implies that 𝒲_df is defined with four subtraction terms, given by coupling the current to each of the four external propagators. When 𝒲_df is then projected to definite angular momenta, these terms become pairwise redundant leading to the factors of 2. In the case of non-identical particles that both couple one must sum over the two choices of species mass within Δ(k_f^2) as well as the alternative mass assignments within D(m, k), i.e. {m_1, m_2}→{m_2, m_1}.

Using Eq. (<ref>), we can write the first term on the right-hand side in terms of G_μ_1 ⋯μ_n; ℓ_f m_f ;     ℓ_i m_i(P_f,P_i,L), defined in Eq. (<ref>). Following the steps outlined in Ref. <cit.>, one finally arrives at Eq. (<ref>) – the relationship between the finite-volume matrix element and i_ df;μ_1 ⋯μ_n, defined in Eq. (<ref>). The latter emerges through the identity δ_df𝒲δ_df≡𝒲_df.








 


We close by giving explicit expressions for the case that the current has a non-negligible coupling to both of the two particles. When two distinct species, 1 ≠ 2, each admit a 1 + 𝒥→1 transition, then the definition of _ df, Eq. (<ref>), is replaced by an expression with four subtractions

    _ df;μ_1 ⋯μ_n≡_μ_1 ⋯μ_n - 
    i  (P_f,k',k) i/(P_f-k)^2 - m_1^2_1;μ_1 ⋯μ_n
    -
     _1;μ_1 ⋯μ_ni/(P_i - k')^2 - m_1^2
    i '(P_i, k',k)  
    
    - 
    i  (P_f,k',P_f - P_i + k) i/(P_f-P_i+k)^2 - m_2^2_2; μ_1 ⋯μ_n
    -
     _2; μ_1 ⋯μ_ni/  (P_i - P_f + k')^2 - m_2^2
    i '(P_i, P_i - P_f + k',k)   .

Here the _2; μ_1 ⋯μ_n in the bottom line have arguments (k', P_i - P_f + k') and (k', P_i - P_f + k') respectively. We have no freedom to choose these once the external momenta are fixed.  

The extra subtractions also lead to an additional G-dependent term in the relation between    and . But in this expression, as a result of the sum, we do have the freedom to re-label the coordinates. With this in mind it is most convenient to rewrite _2; μ_1 ⋯μ_n as a function of (P_f - k, P_i -  k) and decompose exactly as in the main text

    _2;μ_1 ⋯μ_n (P_f - k, P_i -  k)    ≡∑_j K^(j)_2;μ_1 ⋯μ_n(m,k,P_f,P_i)   |_k^0=ω_k1 f_2^(j)(Q^2)  , 
    K^(j)_2;μ_1 ⋯μ_n(m,k,P_f,P_i)    ≡∑_n'=0^n   K^ω_2; μ_1⋯μ_n'(m,k)    C^(j)_2;μ_n'+1…μ_n(P_f,P_i)  .

This then allows us to write 

    ^μ_1 ⋯μ_n(P_f,P_i,L)  - ^μ_1 ⋯μ_n(s_i,s_f,Q^2) ≡
    ∑_j
      ∑_n'=0^n 
      C_1^(j), μ_n'+1…μ_n(P_f,P_i) f_1^(j)(Q^2)
    
    
     [ ℳ(s_f)  G_12^μ_1 ⋯μ_n'(P_f,P_i,L)  ℳ(s_i)  ] 
    
     +
     ∑_j
      ∑_n'=0^n 
      C_2^(j), μ_n'+1…μ_n(P_f,P_i) f_2^(j)(Q^2)
     [ ℳ(s_f)  G_21^μ_1 ⋯μ_n'(P_f,P_i,L)  ℳ(s_i) ]
        ,
 
 where G_12 is exactly equal to the quantity defined in Eqs. (<ref>) and (<ref>) of the main text and G_21 is the same but with m_1 ↔ m_2 everywhere.
 
Turning to the case of identical particles here the relation between between _ df  and  again has four subtractions, exactly as in Eq. (<ref>). The four terms continue to be distinct due to the four different momentum assignments. However, the relation between    and  is exactly as Eq. (<ref>) of the main text, a single G-dependent term with no symmetry factors. This follows from the fact that identical particles lead to a unique diagram of the form shown in Fig. <ref>(b), and that this has no symmetry factor, even in the case that the three hadrons in the loop are identical.

 







 



§ DETAILS FOR EVALUATING THE FINITE-VOLUME FUNCTIONS  


In the following subsections we collect various details relevant for the evaluation of the two finite-volume functions that enter our formalism F_ℓ m ;   ℓ' m'(P,L) and G^σ(P_f,P_i,L).



 §.§ Index gymnastics 


We first discuss various identities for rearranging spherical-harmonic and Lorentz indices in the evaluation of G_σ (P_f,P_i,L). Begin with the case of P_i=P_f, in particular with Eq. (<ref>) of the main text. Multiplying both sides by (q^⋆ )^ℓ_i + ℓ_f and substituting the definition of 𝒴_ℓ m(k^⋆), this reduces to 

    4 π (k^⋆)^ℓ_i + ℓ_f  Y_ℓ_f m_f(k̂^⋆) 
      k_μ_1⋯  k_μ_n   Y^*_ℓ_i m_i(k̂^⋆)    |_k^0 = ω_k2 
      =  √(4 π)∑_J M𝒞_σ, JM(β , k^⋆)     k^⋆ JY_J M(k̂^⋆)  .

The first aim of this subsection is to use this result to derive a useful expression for evaluating 𝒞_σ, JM(β , k^⋆).
  
Substituting the relation k_μ = [Λ_-β]_μ^ν  k^⋆_ν, 
we reach

    [Λ_-β]_μ_1^ν_1⋯[Λ_-β]_μ_n^ν_n [ 4 π (k^⋆)^ℓ_i + ℓ_f   Y_ℓ_f m_f(k̂^⋆) 
      k^⋆_ν_1⋯  k_ν_n^⋆   Y^*_ℓ_i m_i(k̂^⋆)   ]   |_k^0 = ω_k2 
      =  √(4 π)∑_J M𝒞_σ, JM(β , k^⋆)     k^⋆ JY_J M(k̂^⋆)  .

At this stage, the factor in square brackets is a simple polynomial in the coordinates k_x^⋆, k_y^⋆ and k_z^⋆ with additional dependence on the magnitude entering through ω_k2^⋆. Thus, for a given set of indices, one can readily determine an explicit expression, and then use the orthogonality of the Y_ℓ ms, to deduce

    𝒞_σ, JM(β , k^⋆) = 
    [Λ_-β]_μ_1^ν_1⋯[Λ_-β]_μ_n^ν_n √(4 π) (k^⋆)^ℓ_i + ℓ_f - J∫ d Ω^⋆Y^*_J M(k̂^⋆)   [ Y_ℓ_f m_f(k̂^⋆) 
      k^⋆_ν_1⋯  k_ν_n^⋆  Y^*_ℓ_i m_i(k̂^⋆)       ] |_k^0 = ω_k2 .


As a specific example, we return to the case of σ = [μ=z; 10; 10] and P = [00d_z], already discussed in the main text.  For these indices, Eq. (<ref>) reduces to

    𝒞^σ_JM(β , k^⋆) =    3/√(4 π) (k^⋆)^2 - J 2 πδ_M,0∫_-1^1 d cosθ^⋆ Y^*_J 0(θ^⋆)   [ cos^2θ^⋆ (
     β^z γ ω_k2^⋆  + γ   k^⋆cosθ^⋆ )     ]   .

Here we have used the fact that, with μ fixed to z, we only need to sum over one row of the boost matrix, [Λ_-β]^z_ν
= (β^z γ, 0, 0, γ).  The integrals are now trivial to evaluate. For example, the JM=00 component reduces to

    𝒞^σ_00(β , k^⋆) = β^z γω_k2^⋆3/  4 π (k^⋆)^2  2 π∫_-1^1 d cosθ^⋆ cos^2 θ^⋆ 
     = k^⋆2ω_k2^⋆P^z/E^* ,

where we have used β^z = P^z/E and γ = E/E^⋆. This matches Eq. (<ref>) of the main text.

 An alternative approach for determining 𝒞_σ, JM(β , k^⋆) is to note that the product of two spherical harmonics can be written in terms of Clebsch-Gordon coefficients
 
    Y_ℓ_1 m_1(k^⋆) Y_ℓ_2 m_2(k^⋆)
       =
    ∑_ℓ_3 m_3√((2ℓ_1+1)(2ℓ_2+1)/(2ℓ_3+1) 4π)⟨ℓ_10ℓ_20|ℓ_30⟩⟨ℓ_1m_1ℓ_2m_2|ℓ_3m_3⟩Y_ℓ_3 m_3(k^⋆)  ,

and that each factor of k^⋆_ν can be rewritten using the identity
    
    k^⋆ν    = √(4π)(ω^⋆_k2Y_00(k̂^⋆),
    k^⋆Y_1-1(k̂^⋆)-Y_11(k̂^⋆)/√(6),
    k^⋆ Y_1-1(k̂^⋆)+  Y_11(k̂^⋆)/-i√(6),
    k^⋆Y_10(k̂^⋆)/√(3))
    
    ≡∑_ℓ m<2T^ν_ℓ m(k^⋆)   Y_ℓ m(k̂^⋆)  ,

where the last equality defines T^ν_ℓ m(k^⋆).
Substituting this into the left-hand side of Eq. (<ref>) leads to sums over products of spherical harmonics on that side of the equation. These can be pairwise combined using Clebsch-Gordon coefficients, until the left-hand side is reduced to a sum over a single harmonic. Then one can use the orthogonality of spherical harmonics to match this, term by term, to the right-hand side and thereby determine the values of 𝒞_σ, JM(β, k^⋆).



We next consider the case of P_i ≠ P_f. As explained in the main text, here we find it more useful to convert indices in the other direction, i.e. to trade all dependence on spherical harmonics for additional momentum 4-vectors. The key distinction between this case and that discussed above is that we no longer have a natural c.m. frame. The rest frames of P_i and P_f differ, and expressing the integrand in either frame leads to ugly expressions. 
This, together with the need to reach covariant expressions that we can evaluate semi-analytically, led us to introduce M_σ^ν_1 ⋯ν_N, in Eq. (<ref>) above. The definition can be re-expressed as

    M^ν_1 ⋯ν_N_ [μ_1 ⋯μ_n;  ℓ_f m_f ;     ℓ_i m_i ]  k_ν_1⋯ k_ν_N    = 4 π/(q^⋆_f)^ℓ_f(q^⋆_i)^ℓ_i [ ( k^⋆_f  )^ℓ_f    Y_ℓ_fm_ f(k̂^⋆_f) ]
          k_μ_1⋯ k_μ_n [  ( k^⋆_i  )^ℓ_i   Y^*_ℓ_i m_ i(k̂_i) ]  .

The second, and final aim of this subsection, is to derive a useful expression for M_σ^ν_1 ⋯ν_N.

First we introduce a new set of tensors, denoted 𝒯, that allow us to express the right-hand side in terms of 4-vectors

    M^ν_1 ⋯ν_N_ [μ_1 ⋯μ_n;  ℓ_f m_f ;     ℓ_i m_i ]  k_ν_1⋯ k_ν_N    =  𝒯^α_1 ⋯α_ℓ_f_[ℓ_f m_f  ]𝒯^*γ_1 ⋯γ_ℓ_i_[ℓ_i m_i  ]  k^⋆ f_α_1⋯  k^⋆ f_α_ℓ_f        k_μ_1⋯ k_μ_n    k^⋆ i_γ_1⋯  k^⋆ i_α_ℓ_i .

The exact definition of the  𝒯^α_1 ⋯α_ℓ_f_[ℓ m   ] can be inferred by comparing Eqs. (<ref>) and (<ref>).
 Note that they contain the √(4 π)/(q^⋆)^ℓ prefactors and also encode the combinations of k^⋆ components needed to form the various spherical harmonics. For example √(4 π) k^⋆ Y^*_10(k̂^⋆)/q^⋆ = -√(3) k_μ=3^⋆/q^⋆ implies 𝒯^μ_[10] = √(3) (0,0,0,-1)/q^⋆, since k_μ^⋆ =(ω^⋆_k2, -k^⋆). The single-index tensor, 𝒯^α_[ℓ  m   ], is closely related to T^ν_ℓ m, introduced above.  
 
 The final step 
  is to boost all 4-vectors to the finite-volume frame. We deduce  

    M^ν_1 ⋯ν_N_ [μ_1 ⋯μ_n;  ℓ_f m_f ;     ℓ_i m_i ]   =     𝒯^α_1 ⋯α_ℓ_f_[ℓ_f m_f  ]𝒯^*γ_1 ⋯γ_ℓ_i_[ℓ_i m_i  ]
    ×[Λ_β_f]_α_1^ν_1⋯[Λ_β_f]_α_ℓ_f^ν_ℓ_f×δ_μ_1^ν_ℓ_f+1⋯δ_μ_n^ν_ℓ_f+n×[Λ_β_i]_γ_1^ν_ℓ_f+n+1⋯[Λ_β_i]_γ_ℓ_i^ν_N .






 §.§ Evaluating F_ℓ m ;   ℓ' m'(P,L) 



We next turn to the finite-volume matrix, F(P,L). For convenience we repeat the definition given in Eq. (<ref>) above 

    F_ℓ m ;   ℓ' m'(P,L)  
       ≡ξ [1/L^3∑_𝐤 - ∫d^3 k/(2π)^3 ]
     𝒴_ℓ m(k^⋆ )    𝒴^*_ℓ' m'(k^⋆)   
      /2 ω_Pk1 2 ω_k2(E - ω_Pk1-ω_k2  + i ϵ ) .

Our aim here is to rewrite this in terms of 𝒵^(1). With this in mind, we first observe
   
    F_ℓ m ;   ℓ' m'(P,L)      =
    ξ/2  E^⋆∑_JM B_JM^[ℓ m; ℓ' m']1/(q^⋆)^J [1/L^3∑_𝐤 - ∫d^3 k/(2π)^3 ]  ( k^⋆/q^⋆ )^ℓ + ℓ' - J2ω_k2^⋆/2ω_k2√(4 π)  (k^⋆)^J Y_JM(k̂^⋆ )   / q^⋆ 2 - k^⋆ 2 +iϵ ,

 where we have used the fact that 2 ω_Pk1 (E - ω_Pk1-ω_k2  + i ϵ ) can be replaced with 2 E^⋆ (q^⋆ 2 - k^⋆ 2 +iϵ)/2 ω_k2^⋆ up to smooth terms in the integrand that lead to exponentially suppressed volume dependence.   The definition of B_JM^[ℓ m; ℓ' m'] can be inferred by comparing Eqs. (<ref>) and (<ref>) and is given more explicitly by
 
    B_JM^[ℓ m; ℓ' m']   ≡√(4 π)∫ d Ω^⋆   Y^*_JM(k̂^⋆)    [ Y_ℓ m(k̂^⋆)          Y^*_ℓ' m'(k̂^⋆)  ]   ,
       
     =
    (-1)^m'√((2ℓ+1)(2ℓ'+1)/(2J+1))⟨ℓ 0ℓ'0|J 0⟩⟨ℓ mℓ -m'|JM⟩,

where we have used Eq. (<ref>) to reach the last equality. 
  
The final step here is to note that the additional barrier factor, (k^⋆/q^⋆)^ℓ + ℓ' - J, appearing in Eq. (<ref>), can in fact be set to 1. This is justified because (k^⋆)^J Y_JM(k̂^⋆ ) is an analytic function and because ℓ + ℓ' - J ≥ 0 for all nonzero B_JM^[ℓ m; ℓ' m']. The latter point directly follows from the explicit expression of the B_JM^[ℓ m; ℓ' m'] in terms of the Clebsch-Gordon coefficients. As a result, the difference (k^⋆/q^⋆)^ℓ + ℓ' - J - 1 cancels the pole, leading to a smooth summand and a suppressed sum-integral difference.  

 Removing this extra factor, and re-expressing the sum and integral with dimensionless coordinates, we conclude
 
    F_ℓ m ;   ℓ' m'(P,L)     = 
    ξ/  8 π^2 L E^⋆∑_JM B_JM^[ℓ m; ℓ' m'](2 π)^J/ (q^⋆ L)^Jlim_α→ 0𝒵^(1)_J M(P,L,α)  .





 

 §.§ α-dependence of  𝒵^(n)_J M

Here we explain our choice of cut-off function used in the definition of 𝒵^(n)_J M, Eq. (<ref>), recalled here for convenience

    𝒵^(n)_J M(P,L,α) =     [ ∑_𝐫ω_k2^⋆/ω_k2 - ∫   d^3 r^⋆]    
      √(4 π)  r^⋆ JY_J M(r̂^⋆)     / (x^2-r^⋆2+ i ϵ)^n   e^- α (r^⋆2 - x^2)^n .

For n ≥ 2 the integral and sum are individually convergent in the limit α→ 0. Nonetheless, evaluating the sums for various non-zero α and extrapolating α→ 0 turns out to be more efficient than saturating the α = 0 expression directly.  

By including the power of n in the cutoff function, e^- α (r^⋆2 - x^2)^n, we ensure that differentiating with respect to α gives a smooth summand. This, in turn, implies that ∂_α𝒵^(n) vanishes up to terms that are exponentially suppressed in the volume
 
    ∂𝒵^(n)_J M(P,L,α)  /∂α   =    
    -    (-1)^n [ ∑_𝐫ω_k2^⋆/ω_k2 - ∫   d^3 r^⋆]    
      √(4 π)  r^⋆ JY_J M(r̂^⋆)      ,    
       =    
    -    (-1)^n
    (L/2π)^J
    ∑_𝐧≠0∫d^3 k/(2π)^3ω_k2^⋆/ω_k2
       e^i L n·k√(4 π)    k^⋆ JY_J M(k̂^⋆)       =   𝒪( e^-mL)  .

Using this result in an expansion about α = 0 then gives
 
    𝒵^(n)_J M(P,L,α)     =     𝒵^(n)_J M(P,L,0) + 𝒪(α e^-mL )  .



Of course, it is possible to define the 𝒵^(n)_JM functions with a milder cutoff function, for example

    𝒵^(n)_J M(P,L,α) ≡ [ ∑_𝐫ω_k2^⋆/ω_k2 - ∫   d^3 r^⋆]    
      √(4 π)  r^⋆ JY_J M(r̂^⋆)     / (x^2-r^⋆2+ i ϵ)^n   e^- α (r^⋆2 - x^2) .

For n=1 this equivalent to 𝒵^(n)_J M(P,L,α), but for n > 1 it is a less useful prescription, due to an enhancement of the α corrections

 
    𝒵^(n)_J M(P,L,α) -𝒵^(n)_J M(P,L,0) 
       =    
    α  [ ∑_𝐫ω_k2^⋆/ω_k2 - ∫   d^3 r^⋆]    
    √(4 π)  r^⋆ JY_J M(r̂^⋆)     / (x^2-r^⋆2+ i ϵ)^n-1 +𝒪(α^2) ,
     
       =    
    α 𝒵^(n-1)_J M(P,L,0) +𝒪(α^2) .

As 𝒵^(n-1)_J M(P,L,0) is itself a singular function for n > 1, we deduce that the difference between the optimal version, 𝒵^(n)_J M(P,L,α), and the alternative, 𝒵^(n)_J M(P,L,α) can take on arbitrarily large values for any finite α.


 
 In Fig. <ref> we compare the α-dependence of 𝒵^(n)_J M(P,L,α) and 𝒵^(n)_J M(P,L,α) for n=1,2, and show that the large α-dependence of the latter is well described by Eq. (<ref>). In the α→ 0 limit the two prescriptions agree, but to optimize the numerical evaluation we advocate the form of Eq. (<ref>) and use only this definition throughout the remainder of the text.
 
 
 
 









 §.§ Evaluating 𝒵^(n)_J M(P,L,α) 

 
 In this section we discuss the evaluation of 𝒵^(n)_J M(P,L,α), in particular the integral part of this quantity. As already mentioned in Sec. <ref>, the integral entering 𝒵^(n) vanishes for all entries besides J M = 00, implying
 
    𝒵^(n)_J M(P,L,α)     =      ∑_𝐫ω_k2^⋆/ω_k2√(4 π)    (r^⋆)^ JY_J M(r̂^⋆)     / ( x^2-r^⋆2+ i ϵ)^n   e^- α (r^⋆2 - x^2)^n -δ_J0δ_M0Ξ^(n)(x, α) ,

where

    Ξ^(n)(x, α) ≡
    4π∫^∞_0  d r^⋆ r^⋆ 2   e^- α ( r^⋆2- x^2)^n   / (x^2-r^⋆2+ i ϵ  )^n   .

 
For the single-pole function, 𝒵^(1), this can be evaluated analytically and takes the form

    Ξ^(1)(x, α) =
    
    
    4π[ - √(π/4 α)
    e^α  x^2
    +
    π  x/2 Erfi(√(α  x^ 2))
    -i π x/2] ,

 where Erfi(x) is the imaginary error function function, defined by d  Erfi(z)/dz= 2  e^z^2/√(π) and Erfi(0) = 0. 

  
 

We now consider n > 1. Although we only require n=2 for the present work, we find it instructive to consider all values together. 
A consequence of the cutoff function, together with the higher pole factors, is that we are not able to evaluate  Ξ^(n)(x, α) analytically for n > 1. Instead, following our usual trick, we separate the expression into two terms: one that can be evaluated analytically and another that is smooth and converges rapidly under numerical integration


    Ξ^(n)(x, α)
       =
    4π∫^∞_0  d r^⋆ r^⋆ 2   1   / (x^2-r^⋆2+ i ϵ    )^n  
    +
    δΞ^(n)(x, α)  , 
       =
    -i2π^2/(n-1)!(-∂/∂ x^2)^n-1√(x^2)
    +
    δΞ^(n)(x, α)  ,

where

    δΞ^(n)(x, α) ≡ 4π∫^∞_0  d r^⋆ r^⋆ 2    e^- α ( r^⋆2- x^2)^n    - 1    / (x^2-r^⋆2    )^n .


We close with a final remark concerning the n=2 case, of direct relevance for 2 + 𝒥→2 transition amplitudes. Here the relevant integral is

    Ξ^(2)(x, α)
       =
    iπ^2/x
    +
    δΞ^(2)(x, α)  .

Recalling x ∝ q^⋆∝√(s-s_ th), where s = P^2 is the c.m. energy and s_ th = (m_1 + m_2)^2, we deduce that for P_i = P_f, G_σ generically has an inverse square-root singularity at two-particle production threshold. This implies that 𝒲_df, as well as 𝒲, must have the same singularity. This behavior is visible in the values of 𝒵_JM^(2) plotted in Fig. <ref>. In particular, we observe that 𝒵_JM^(1) has milder behavior near threshold.



 

 §.§ Symmetry contraints on 𝒵^(n)_J M

 





To efficiently implement the formalism it is useful to identify, based on symmetry arguments, the values of JM for which 𝒵^(n)_J M(P,L,α) = 0. In this subsection we review these constraints and discuss subtleties that arise for n > 1. 
Our results are summarized in Table <ref>.

We begin with P = [000]. The properties of the zero-momentum zeta function are well-known <cit.>, but we still think it useful to review the arguments here, in order to better understand the  generalization to P≠ [000] and n>1.
Note that the zero-momentum zeta function must be unchanged if we flip r everywhere in the summand. Taking the expression for JM ≠ 00 we write
 
    𝒵^(n)_J M((E, 0),L,α)     =      ∑_𝐫√(4 π)    (| - r|)^ JY_J M(- r̂)     /( x^2- (- r )^2  )^n  e^- α ((- r )^2 - x^2)^n .

Substituting Y_J M(- r̂) = (-1)^J Y_J M( r̂) then gives 𝒵^(n)_J M((E, 0),L,α)  = (-1)^J 𝒵^(n)_J M((E, 0),L,α) implying that the zeta function vanishes for all odd J. We can further rewrite the summand with (r_x,r_y,r_z) ⟶ (-r_x,r_y,r_z) and use Y_JM(θ,  π - ϕ) =  Y_J,-M(θ,ϕ) to show that 𝒵^(n)_J  M((E, 0),L,α)  =  𝒵^(n)_J , -M((E, 0),L,α). Similarly, a π/2 rotation about the z-axis, together with the identity Y_JM(θ,  ϕ - π/2) = e^- i M π/2 Y_JM(θ,ϕ) implies that M must be divisible by 4 for the zeta function to be non-zero. 

A final zero-momentum case worth mentioning is JM=20. To see that this vanishes as well, note that the corresponding spherical harmonic is proportional to 3 r_z^2 - r^2. The sum of this structure (times a function of r^2) over the integer set, r∈ℤ^3, is clearly identical to the same with 3 r_x^2 - r^2 or with 3 r_y^2 - r^2. Thus 𝒵^(n)_20 is equally well defined by averaging the three possibilities. But this gives a summand proportional to 3 ( r_x^2 + r_y^2 + r_z^2 ) - 3 r^2 = 0, implying 𝒵^(n)_20 = 0 as claimed.
The conditions under which 𝒵^(n)_J M = 0 for  P = [000]  are summarized in the first line of Table <ref>.




 






We now turn to non-zero momenta of type P = [00 d_z]. As this momentum type only breaks the symmetry in the z direction, the invariance under (r_x,r_y,r_z) ⟶ (-r_x,r_y,r_z), as well as the π/2 rotation around the z-axis, give the same constraints as for  P = [000].[Of course, it is true in all cases that the sum over r∈ℤ^3 is invariant under any octahedral transformation on r. The relevant question is whether this leads to a useful constraint on 𝒵_JM^(n).] Similarly, for P = [0 d_y d_z], flipping only r_x gives the same relation as above. 

By contrast, parity is broken for any non-zero momentum so that the argument based on r→ - r no longer holds. 
For example, for m_1 = m_2 and P = [00d_z], the summand defining 𝒵^(n)_J M  now depends on 

    r^⋆ =  (r_x,   r_y,   -  2 π d_z /E^⋆ L [ m^2 L^2/4 π^2 + r^2  ]^1/2  + E/E^⋆ r_z )  .

This vector does not transform in a useful way under a flip of r.
Remarkably, in the degenerate-mass case, the single pole functions, 𝒵^(1)_JM, continue to vanish for all odd J.
More precisely, these are smooth functions with exponentially suppressed volume dependence and thus scale as terms that have been dropped in the derivation. As we now explain, this is due to an accidental symmetry  
 inherited from the non-relativistic  system.  

 

The following argument holds for all values of P = (2 π/L) d and so we present the results for the general case. The approach is based on the results of Ref. <cit.> and we begin by recalling Eqs. (62) and (66) from that work 

    r^⋆_∥    = R^⋆_∥ -  P/EL/2 πk^⋆2 - q^⋆2/ω_k^⋆ + E^⋆/2 ,     r^⋆_ =  R^⋆_ , 
    
    ( x^2-r^⋆2+ i ϵ)     = γω_k2^⋆ /ω_k2 (x^2 - R^⋆2 + i ϵ) + 𝒪[(x^2 - R^⋆2)^2]  .

where R^⋆ =  γ̂^-1(r - d/2). Here r^⋆_∥ and r^⋆_ are the vector components parallel and perpendicular to d. We have also introduced the operator γ̂^-1 (r^⋆_∥ , r^⋆_ ) = (   γ^-1r^⋆_∥ , r^⋆_ ). Substituting these results into our definition of 𝒵_JM^(n) we find
 
    𝒵^(n)_J M(P,L,α)     =      ∑_𝐫 [ ω_k2^⋆/ω_k2 (  γω_k2^⋆/ω_k2 )^-n√(4 π)    (R^⋆)^ JY_J M(R̂^⋆)     /( x^2-  R^⋆ 2+ i ϵ  )^n
     + 𝒪((x^2 - r^⋆2 )^1-n)   ] .


Now note that, for n=1, this function exhibits two special features, both unique to the single-pole case. First, the factors of ω_k2/ω_k2^⋆ multiplying the pole exactly cancel; second, the subleading term becomes a smooth function of the summed coordinate, r. We thus reach
 
    𝒵^(1)_J M(P,L,α)     =     1/γ∑_𝐫√(4 π)    (R^⋆)^ JY_J M(R̂^⋆)     /( x^2-  R^⋆ 2+ i ϵ  )
     + 𝒪(e^-m L) .

This simplified form, incidentally the form first derived for the moving frame quantization condition <cit.>, makes the accidental symmetry that we are after manifest. In particular, we can now use that the sum over r∈ℤ^3 is invariant under r⟶d - r. Under this transformation R^⋆→ - R^⋆ leaving R^⋆ = |R^⋆| unchanged. Thus every factor in the summand is invariant except for Y_J M(- R̂^⋆) = (-1)^J Y_J M( R̂^⋆). We deduce that, for odd J and degenerate masses, 𝒵^(1)_J M=𝒪(e^- m L). However, the remaining symmetries do not survive due to the factor of γ̂^-1 in the definition of R^⋆. 


Finally we stress that the vanishing of odd J due to the accidental symmetry only holds for n=1. As is clear from Eq. (<ref>), for all other n values, the ratio of omegas does not cancel, leading to another factor that changes under the r→d - r transformation. Thus, while the odd-J 𝒵^(1)_JM have no poles, for n=2 the functions already exhibit the full double-pole behavior. [See also Fig. <ref>.] In addition, for n>1, the subleading term contains a (n-1)th order pole that also generates an important contribution to the zeta function.
 

At this stage we have completed all details relevant for P_i = P_f and therefore turn to the more complicated case of P_i ≠ P_f, beginning with the numerical integral denoted ℐ_𝒩; σ(α, P_f, P_i).






 §.§ Evaluating ℐ_𝒩; σ(α, P_f, P_i) 


Here we give some further details on the evaluation of ℐ_𝒩; σ(α, P_f, P_i), defined in Eq. (<ref>) of the main text. Recall that this is a 3-dimensional integral with a smooth integrand, to be evaluated numerically. When added to the semi-analytic expression, ℐ_𝒜; σ, it gives the full integral entering the definition of G_σ.

The complication we wish to address here is that the definition in the main text requires evaluating a large number of integrals defined with factors of k_μ (but with no spherical harmonics) in the integrand. These are then contracted with the tensor M^ν_1 ⋯ν_N_σ to reach the final expression. While highly advantageous for the analytic integral, the Lorentz-index-based expressions lead to a costly determination of ℐ_𝒩; σ. For example, for ℓ_f = ℓ_i = 1 and a single current index, the covariant form contains factors of k_ν_1 k_ν_2 k_ν_3 that naively require the evaluation of 64 terms. 



To improve this situation, it is preferable to move M^ν_1 ⋯ν_N_σ back inside the integrals defining ℐ_𝒩; σ, and thereby rewrite the integrands to only carry the final σ index. This procedure is a bit subtle, because the time components of the k^μ are evaluated at various values. To proceed we first recall that ℐ_𝒩, σ generally involves terms evaluated at the physical masses, m_1, m_2, together with regulating integrals evaluated at some higher scale Λ_j. For the physical-mass terms, time components are evaluated at k^0 ∈{ω_k2, E_f + ω_P_fk1, E_i + ω_P_i k1}, and for the regulating integrals at k^0 ∈{√(k^2 + Λ_j^2), E_f + √((P_f-k)^2 + Λ_j^2), E_i + √((P_i-k)^2 + Λ_j^2)}. Some of these four vectors are recombined into the harmonics, 𝒴_ℓ m(k^⋆), and, because the time and space components mix upon boosting k^⋆_μ = [Λ_β]_μ^ν k_ν, we end up with strange spatial components in some of the harmonics.




To give concrete expressions it is convenient to define

    𝒩_σ(Ω_k,  Λ, k) ≡ M^ν_1 ⋯ν_N_σ [   k_ν_1⋯ k_ν_N|_k^0 = Ω_k ]  = 𝒴_ℓ_f m_ f(  k^⋆, Ω_f)  
          [   k_μ_1⋯ k_μ_n|_k^0 = Ω_k ]    𝒴^*_ℓ_i m_ i(  k^⋆, Ω_i)       .

Here Ω_k represents any of the possible choices made for the temporal component of the 4-vectors. 
In each of these three cases, an implicit mass dependence enters and the second argument, Λ, refers to this mass dependence. 
In the following we will use m or Λ_0 to indicate that the ωs are evaluated at their physical masses and Λ_j>0 to indicate evaluation at an unphysical value m_1 = m_2 = Λ_j>0. In short, the first two entries in 𝒩_σ simply serve to indicate the value at which all k^0 are evaluated.
We stress that the tensor M^ν_1 ⋯ν_N_σ does not depend on these parameters but only on the c.m. frame energies E_i^⋆ and E_f^⋆ as well as the boost velocities β_i and β_f. Thus, the only modification to the spherical harmonics is that they now depend on k^⋆, Ω, defined via

    (Ω_k^⋆,  k^⋆, Ω)^μ≡ [Λ_β]^μ_ν  (Ω_k,  k )^ν .

If we set Ω_k = ω_k2 and Λ = m, then we exactly recover the spherical harmonic definitions used everywhere else in this work. 

With our new numerator function  in hand, we are ready to give our final form for ℐ_𝒩; σ

    ℐ_𝒩; σ(α, P_f, P_i)    ≡∫d^3 k/(2π)^3  
     [H(α, k) - 1]
      ∑_j=0^n_j c_j  [ D(Λ_j, k) 𝒩_σ  (ω_k , Λ_j, k) + 𝒦_r;σ(Λ_j, k)   ]  
     
    
       
     -  ∫d^3 k/(2π)^3   H(α, k)  ∑_j=1^n_j  c_j D(Λ_j, k)  𝒩_σ  (ω_k , Λ_j, k)
          
    
       
     -  ∫d^3 k/(2π)^3   H(α, k)   ∑_j=0^n_j c_j     𝒦_r;σ(Λ_j, k)    
      ,

where

    𝒦_r;σ(Λ, k) ≡ D_rf(Λ, k) 𝒩_σ(E_f+ω_P_fk1, Λ, k) +  D_ri(Λ, k) 𝒩_σ(E_i+ω_P_i k1, Λ, k)   .

This is identically equal to the quantity defined in Eq. (<ref>) of the main text. The only difference is that we have absorbed M^ν_1 ⋯ν_N_σ inside the integrands, via the new function 𝒩_σ.

To complete the specification we require explicit expressions for D_rf(m, k) and D_ri(m, k)

    D_rf(m, k)    =
    1/2 ω_P_fk11/ ( E_f+ω_P_fk1)^2 - ω_k2^2   1/  (E_i-E_f-ω_P_fk1)^2 - ω_P_ik1^2    , 
    
    D_ri(m, k)    =
    1/2 ω_P_ik11/ ( E_i+ω_P_ik1)^2 - ω_k2^2    1/  (E_f-E_i-ω_P_ik1)^2 - ω_P_fk1^2    .

The motivation for these quantities is discussed in the text around Eq. (<ref>), where an implicit definition is also given.



We close this subsection with one final  simplification to ℐ_𝒩. We show that it is always possible to simplify the numerical evaluation from a three- down to a two-dimensional integral. To show this we must prove the following key identity

    ℐ^μ_1 ⋯μ_n_𝒩; ℓ_f m_f;ℓ_i m_i(α, P_f,P_i) =
     [Λ_-β]^μ_1_ν_1⋯[Λ_-β]^μ_n_ν_n ℐ^ν_1 ⋯ν_n_𝒩;ℓ_f m_f;ℓ_i m_i(α,Λ_βP_f,Λ_βP_i )  .


Begin the proof by considering the generic integral

    𝒳^μ_1 ⋯μ_n(P_f, P_i) ≡∫d^4 k/(2 π)^4ℱ(P_f, P_i, k) 𝒢(k_i^⋆, k_f^⋆) k^μ_1⋯ k^μ_n ,

assumed to be convergent. Here we have separated the integrand into a Lorentz-scalar function, ℱ, together with an arbitrary function of the c.m. frame momenta, 𝒢. The latter is also Lorentz invariant in the vacuous sense, i.e. because its arguments carry a frame label. Now act on both sides with Λ_β, on each index

    [Λ_β]^μ_1_ν_1⋯[Λ_β]^μ_n_ν_n𝒳^ν_1 ⋯ν_n(P_f, P_i) ≡∫d^4 k/(2 π)^4ℱ(P_f, P_i, k) 𝒢(k_i^⋆, k_f^⋆)  [Λ_β]^μ_1_ν_1⋯[Λ_β]^μ_n_ν_n  k^ν_1⋯ k^ν_n .

To simplify the right-hand side we perform a change of integration variable k'^μ≡[Λ_β]^μ_ν k^ν and also define P_f'^μ≡[Λ_β]^μ_ν P_f^ν and P_i'^μ≡[Λ_β]^μ_ν P_i^ν 

    [Λ_β]^μ_1_ν_1⋯[Λ_β]^μ_n_ν_n𝒳^ν_1 ⋯ν_n(P_f, P_i) ≡∫d^4 k'/(2 π)^4ℱ(P_f', P_i', k') 𝒢(k_i^⋆, k_f^⋆)     k'^ν_1⋯ k'^ν_n .

Here we have used the Lorentz-invariance of the various building blocks, including the fact that k_i^⋆ and k_i^⋆ must be unchanged if we replace k, P_f and P_i with their boosted counterparts. This result, which can be rewritten as

    [Λ_β]^μ_1_ν_1⋯[Λ_β]^μ_n_ν_n𝒳^ν_1 ⋯ν_n(P_f, P_i) ≡𝒳^ν_1 ⋯ν_n(P'_f, P'_i)  ,

is just a statement of Lorentz-covariance for 𝒳.

To conclude our demonstration of Eq. (<ref>) we note that ℐ_𝒩≡ℐ - ℐ_𝒜 and that the two terms on the right-hand side each satisfy the functional form of 𝒳, shown in Eq. (<ref>). In the case of ℐ, the original integral defining G, one takes

    ℱ(P_f, P_i, k) =  Θ(k^0) (2 π) δ(k^2 - m_2^2) 1/  (P_f-k)^2 - m_1^2 + i ϵ1/  (P_i-k)^2 - m_1^2 + i ϵ .

This function is only invariant under orthochronous transformations, as is standard when one discards the anti-particle pole, but this is sufficient for the present argument. For both ℐ and ℐ_𝒜 the spherical harmonics, as well as the cutoff function H(α, k) can be absorbed into the definition of 𝒢. Again the key point is that these objects are frame-independent because they carry a frame label, k^⋆ (k, P) = k^⋆ (k', P'). We deduce that ℐ_𝒩 must satisfy Eq. (<ref>). Multiplying both sides by Λ_-β, we conclude Eq. (<ref>).


 
To see the power of this identity we take β = β_i, implying

    [Λ_β_i]^μ_ν P_f^ν≡ P'_f^ν = (E'_f,𝐏'_f) ,        [Λ_β_i]^μ_ν P_i^ν =(E_i^⋆,0),

and thus that only one external direction enters the integral. In this case the integration coordinate is simply transformed to k_i^⋆. We are then left with

    ℐ^μ_1 ⋯μ_n_𝒩; ℓ_f m_f;ℓ_i m_i(α, P_f,P_i) =
     [Λ_-β_i]^μ_1_ν_1⋯[Λ_-β_i]^μ_n_ν_n
    ×∑_k  ∫d^3 k_i^⋆/(2 π)^3𝒬_k(|k_i^⋆|,  k_i^⋆·P_f' ) 𝒴_ℓ_f m_f(k_f^⋆)    k_i^⋆ν_1⋯    k_i^⋆ν_n𝒴^*_ℓ_i m_i(k_i^⋆) |_k^⋆ 0_i = Ω_k ,

where the sum over k runs over all possible choices for the temporal component, as discussed above. 
The key point here is that, once the spherical harmonics and the factors of k^⋆μ_i are factored out, the remaining integrand can only depend on k_i^⋆ through its magnitude and a single angle. This is because no other direction is defined in the system once we have expressed all coordinates in the rest frame of the incoming state.

Defining cosθ≡k̂_i^⋆·P̂_f' and picking an aribtrary orientation for the azimuthal angle ϕ, we reach

    ℐ^μ_1 ⋯μ_n_𝒩; ℓ_f m_f;ℓ_i m_i(α, P_f,P_i) =
     [Λ_-β_i]^μ_1_ν_1⋯[Λ_-β_i]^μ_n_ν_n
    ×∑_k     ∫ d    k_i^⋆  k_i^⋆ 2  d cosθ/(2 π)^3𝒬_k(     k_i^⋆  ,  cosθ)  ∫_0^2 π d ϕ 𝒴_ℓ_f m_f(k_f^⋆)    k_i^⋆ν_1⋯    k_i^⋆ν_n𝒴^*_ℓ_i m_i(k_i^⋆) |_k^⋆ 0_i = Ω_k .

Here the ϕ integral can be evaluated analytically using the rotation properties of the spherical harmonics. Thus only the k_i^⋆ and θ integrals need to be performed numerically.





 §.§ Evaluating ℐ^χ(P_f, P_i, m, δ) 




In this subsection we show how the the generating functional, ℐ^χ(P_f, P_i, m, δ), defined in Eq. (<ref>) of the main text, reduces to the results given in Eqs. (<ref>)-(<ref>). We begin by inserting Feynman-parameter integrals into the definition to reach

    ℐ^χ(P_f, P_i, m, δ)
       
    =  2 ∫_0^1dx∫_0^1-xdy
    ∫d^4 - δ k_E/(2π)^4 - δe^iχ_E· k_E/[ k^2_E+m^2-2 k_E·(x P_f,E+y P_i,E)+x P_f,E^2+y P_i,E^2]^3 , 
       =
    2 ∫_0^1dx∫_0^1-xdy
      e^iχ_E· (x P_f,E+y P_i,E)∫d^4 - δ k_E/(2π)^4 - δe^iχ_E· k_E/[k_E^2+M(m,x,y)^2]^3 ,

where in the second line we have performed the shift k_E → k_E + (x P_f,E+y P_i,E) and have introduced

    M(m,x,y)^2 ≡ m^2+x(1-x) P_f,E^2+y(1-y) P_i,E^2-2  x y P_f,E· P_i,E .


Next note that the denominator of the integrand of Eq. (<ref>) is invariant under k_E → - k_E, implying that only even terms in k_E contribute to the integral. Expanding the exponential and keeping only the even powers, we find

    ℐ^χ(P_f, P_i, m, δ)
       
    = 
    2 ∫_0^1dx∫_0^1-xdy
      e^iχ_E· (x P_f,E+y P_i,E)∫d^4 - δk_E/(2π)^4 - δ1/[k_E^2+M(m,x,y)^2]^3∑_n=0^∞(iχ_E· k_E)^2n/(2n)! .


 To further reduce the expression note that we can make the substitution (i χ_E · k_E)^2n⟶ A_n (-χ_E^2)^n ( k_E^2)^n, where A_n is a normalization constant, to be determined. This holds because the integral over k_E, μ_1⋯ k_E, μ_2n (multiplied by a function of k_E^2) must be proportional to δ_μ_1 μ_2⋯δ_μ_2n-1μ_2n + ⋯ where the second ellipsis indicates a sum over all possible pairings. Contracting with χ_E, μ_1⋯χ_E, μ_2n then gives the claimed form.

To determine the normalization, first consider the case of n=1, corresponding to k_E,μ k_E,ν⟶ A_1 δ_μν k_E^2. Taking the trace on both sides then gives A_1 = 1/(4 - δ), where we are careful to consistently perform the calculation in 4 - δ dimensions. Similarly, for n=2 one finds 

    k_E,μ_1 k_E,μ_2 k_E,μ_3 k_E,μ_4⟶ A_2 (k_E^2)^2 δ_μ_1 μ_2δ_μ_3 μ_4 + δ_μ_1 μ_3δ_μ_2 μ_4  + δ_μ_1 μ_4δ_μ_2 μ_3/3 .

First summing over μ_1 = μ_2 and then over μ_3 = μ_4 gives A_2^-1 =  (1/3) [(4 - δ)^2 + 2 (4 - δ)]. 

The result for general n can be derived by first writing

    A_n  ∫ d Ω_3-δ   ≡∫ d Ω_3-δ(i χ_E · k_E)^2n/(k_E^2)^n (- χ_E^2)^n  
    =
    ∂_α^2n∫ d Ω_3-δexp[ α (i χ_E · k_E)]/(k_E^2)^n (- χ_E^2)^n|_α=0 ,

where in the second equality we have rewritten the integral with a dummy parameter α, to be set to zero after differentiation.  Next we multiply both sides by exp[- k_E^2] and integrate with respect to d k_E k_E^4-δ to write

    A_n  ∫ d ^4-δk_E
    exp[  - k_E^2   ]
       =
    A_n  π^2-δ/2
    =
    ∂_α^2n∫ d ^4-δk_E
     exp[  - k_E^2  + i αχ_E · k_E ]/(k_E^2)^n (- χ_E^2)^n|_α=0.




Solving for A_n and evaluating the remaining integral using Schwinger parameters, we deduce

    A_n 
       =  1  /π^2 - δ/2  (- χ_E^2)^n     ∂_α^2n∫ d^4 - δ k_E  ∫_0^∞ d ββ^n-1/Γ(n)exp[- β k_E^2]  exp[  - k_E^2  + i αχ_E · k_E ] |_α=0 ,
    
    
       =  1  /π^2 - δ/2  (- χ_E^2)^n     ∂_α^2n∫_0^∞ d ββ^n-1/Γ(n)∫ d^4 - δ k_E  
     exp [ 
      - (1 +β) 
      (
      k_E
      +i αχ_E/2(β +1))^2  
    - 
    α^2 χ_E^2/4(β +1)   ] |_α=0 ,
    
    
    
    
    
    
    
    
    
       =  (2n)!  /Γ(n) 4^nn!∫_0^∞ d ββ^n-1(β+1)^-n-2+δ/2 ,
    
    
       =  (2n)!  /    4^n  n!Γ(2-δ/2)/Γ(2+n-δ/2) .

In the second equality we have completed the square in k_E and in the third we have integrated with respect to k_E and also evaluated the α derivative and set α = 0.
 



 Substituting into Eq. (<ref>) then gives

    ℐ^χ(P_f, P_i, m, δ)
     = 
    2 ∫_0^1dx∫_0^1-xdy
      e^- iχ· (x P_f+y P_i)∑_n=0^∞ 
    A_n   (χ^2 )^n/ (2n)! 𝒥^n(P_f, P_i, m, δ)   ,

  where we have returned to the Minkowski signature for χ, P_f and P_i and have also defined
 
    𝒥^n(P_f, P_i, m, δ)    ≡∫d^4 - δk_E/(2π)^4 - δ k^2n_E/[k_E^2+M(m,x,y)^2]^3 .




To conclude we simplify 𝒥^n(P_f, P_i, m, δ) by evaluating the momentum integral 

    𝒥^n(P_f, P_i, m, δ) 
       =
    ∫dΩ_3 - δ/(2π)^4 - δ ∫_0^∞ dk_E k_E^3-δk^2n_E/[k_E^2+M(m,x,y)^2]^3 , 
    
       =
     2π^2 - δ/2/(2π)^4 - δ Γ(2 - δ/2) ∫_0^1  d ζ M(m,x,y)^2/2 ζ^2  
     M(m,x,y)^2n + 2 - δ  (1/ζ - 1)^n + 1 - δ/2ζ^3/M(m,x,y)^6 .

In the second step we have integrated over dΩ_3 - δ and then changed variables via k^2_E ≡ M^2/ζ - M^2. The measure is modified as 2 k_E d k_E = - M^2 d ζ/ζ^2 with the integral now running from ζ = 0 to ζ = 1. The final factor in the second line is just ζ^3/M^6 = [k_E^2 + M^2]^-3.

Evalauting the ζ integral via

    ∫_0^1 d ζ
    (1-ζ)^α-1 ζ^β-1=Γ(α) Γ(β)/Γ(α+β) ,

we conclude
  
    𝒥^n(P_f, P_i, m, δ)   =  1/2(4π)^2 - δ/2Γ(n+2 - δ/2)
    Γ(1-n + δ/2)
    /Γ(2 - δ/2)  M(m,x,y)^2n - 2 - δ .

This directly gives Eq. (<ref>) in the main text.
 
 






 §.§ Determining c_j and Λ_j 



As discussed in Secs. <ref> and <ref>, an immediate application of the proposed formalism is electromagnetic reactions coupling two-pion states: (π^+ π^0)_I_i + j_μ→ (π^+ π^0)_I_f. For this case we require the function G_σ for σ = [ℓ_f m_f;ℓ_i m_i] (no current index) as well as σ = [μ; ℓ_f m_f;ℓ_i m_i] (one index) for ℓ_i, ℓ_f ≤ 1. This requires evaluating ℐ_𝒜 (with no indices) through ℐ_𝒜; ν_1 ν_2 ν_3. This set depends on only two scalar integrals, 𝒥^0 and 𝒥^1. 
The integral defining 𝒥^0 is convergent so that we only need the c_0=1, Λ_0=m term, i.e. no subtraction is required. The integral defining 𝒥^1, by contrast, has a logarithmic divergence (arising from d^4 k_E k_E^2/k_E^6). This is removed using the subtraction given in Eq. (<ref>), corresponding to c_1=-1 with Λ_1 = Λ equal to any value exceeding 2 m. 








 Though the integrals of direct interest are rendered convergent by (up to) one simple subtraction, we think it useful here to give the recipe for general ℐ_𝒜, ν_1 ⋯ν_N. In a nut shell the approach requires identifying the divergent part of 𝒥^n and, by substituting this into the relation to ℐ_𝒜, ν_1 ⋯ν_N, to identify an expression of the form

    ℐ^div_𝒜; ν_1 ⋯ν_N(  P_f, P_i, δ)   ≡∑_j=0^⌊N/2⌋ c_j ℐ^div_𝒜; ν_1 ⋯ν_N(  P_f, P_i,  Λ_j, δ)   .

The coefficients c_j and the scales Λ_j are then tuned such that this quantity vanishes. Note that ⌊N/2⌋ terms must be tuned to vanish. This is because the integral ℐ_𝒜; ν_1 ⋯ν_N depends on the scalar integrals up to n = ⌊N/2⌋. The integral then scales as d^4 k _E k^2n_E / k_E^6 and generates a series of divergences of the form Λ^2 n - 2, Λ^2n - 4, ⋯, log(Λ), so a total of n =  ⌊N/2⌋ terms. 

It turns out that one does not need to tune both the regularization scales Λ_j and the coefficients c_j. We thus choose the recipe of setting Λ_j = 2^j-1Λ [Λ_1 = Λ ,  Λ_2 = 2 Λ ,  Λ_3 = 4 Λ , ⋯] and tuning only the n different c_j terms. 













The latter  is achieved by studying the 1/δ terms. In particular one can show that 𝒥^n ∼ M(m,x,y)^2n - 2/δ. Thus, for n=1 the divergence is M(m,x,y)-independent and is removed by setting c_1=-1 as explained above. For n=2 the divergence scales as M(m,x,y)^2 leading to the linear combination M(m,x,y)^2 + c_1 M(Λ,x,y)^2 - (1 + c_1) M(2 Λ,x,y)^2 where we have already enforced 1 + c_1 + c_2 = 0. This leads to a cancellation in all terms except for the explicit mass and Λ dependence 

    M(m,x,y)^2 + c_1 M(Λ,x,y)^2 - (1 + c_1) M(2 Λ,x,y)^2 = m^2 + c_1 Λ^2 - 4 (1 + c_1) Λ^2  ,

and requiring this to vanish gives c_1 = [ m^2 - 4 Λ^2   ]/[3 Λ^2].  Here we restrict our attention to the degenerate case, m_1 = m_2 = m. While the subtraction is always mass-independent for integrals up to ℐ_𝒜; ν_1 ν_2 ν_3, beyond this the scheme in the non-degenerate case becomes more complicated. 

Note that the choices of coefficients that remove 𝒥^2 divergences also automatically remove those in 𝒥^1. This is because ∑_j c_j M(Λ_j,x,y)^2 contains both Λ_j-dependent and independent pieces. Thus the vanishing of ∑_j c_j M(Λ_j,x,y)^2 guarantees that the same holds for ∑_j c_j, and the latter is the condition for removing divergences in 𝒥^0. 

As we now show, this pattern continues to all orders, so that it is always sufficient to determine c_j by tuning away the highest 𝒥^n divergences. For the general-n, degenerate case, the system that we need to solve is

    ∑_j=1^n c_j     [ M(0,x,y)^2 + 4^j-1Λ^2  ]^n = -    [ M(0,x,y)^2 + m^2  ]^n  ,

or equivalently

    ∑_j=1^n c_j  ∑_k=0^n C(n,k)   M(0,x,y)^2k   4^(j-1)(n-k)Λ^2(n-k)  = -  ∑_k=0^n C(n,k)     M(0,x,y)^2k   m^2(n-k) ,

where C(n,k)= n!/(n-k)! k! is the binomial coefficient. Defining v_k ≡ C(n,k)   M(0,x,y)^2k and A_kj≡   4^(j-1)(n-k)Λ^2(n-k), the above relation becomes v_k A_kj c_j = v_k (- m^2k) with repeated indices summed. Dropping the v_k from both sides, we conclude that a solution is given by

    [           1           1           1           ⋯           1;         Λ^2     (2 Λ)^2     (4 Λ)^2           ⋯  (2^n-1Λ)^2;         Λ^4     (2 Λ)^4     (4 Λ)^4           ⋯  (2^n-1Λ)^4;           ⋮           ⋮           ⋮           ⋱           ⋮;        Λ^2n    (2 Λ)^2n    (4 Λ)^2n           ⋯ (2^n-1Λ)^2n ][ c_1; c_2; c_3;   ⋮; c_n ]
    =
    [     -1; -  m^2; -  m^4;      ⋮; - m^2n ] .

It is straightforward to invert this matrix and read-off the values of c_i to regulate an integral with any number of indices. In Table <ref> we give the values up to n=5, assuming Λ_j = 2^j-1 (3m).




Returning to the case of different masses, here one must instead solve

    ∑_j=1^n c_j  ∑_k=0^n C(n,k) ⟨⟨      M(0,x,y)^2k⟩⟩  4^(j-1)(n-k)Λ^2(n-k)  = - ⟨⟨ [  M(0,x,y)^2 +(1-x-y)m_2^2+ (x+y)m_1^2   ]^n ⟩⟩ ,

where

    ⟨⟨ f(x,y) ⟩⟩≡ 2  ∫_0^1dx∫_0^1-xdy f(x,y) e^iχ· (x P_f,E+y P_i,E) .

The key distinction is that x and y dependence now appears in both M(0,x,y)^2 and the mass terms on the right-hand side. Thus it is not possible to express the right-hand side in terms of a matrix product in which ⟨⟨      M(0,x,y)^2k⟩⟩ is factored off. The upshot is that, in this case, the c_j depend on the kinematics s_i, s_f, q^2, the masses m_1, m_2 and the generating parameter χ. Again we stress that this is not of immediate concern as it is only relevant for form factors with many indices. 




 §.§ Triangle singularities within ℐ_𝒜(P_f, P_i)


In this appendix we give a more detailed discussion of the singularities that arise away from threshold in ℐ_𝒜(P_f, P_i), as summarized around Eqs. (<ref>) and (<ref>) of the main text and in Figs. <ref> and <ref>. 
The task is to study discontinuities of the integral
  
    ℐ_𝒜( P_f, P_i, m, 0)  
    =   1/(4π)^2∫_0^1dx∫_0^1-xdy  1/M(m,x,y)^ 2 .

These arise when the three particles in the triangle loop of Fig. <ref>(b) can all go on shell. As Landau described in Ref. <cit.>, the on-shell condition is realized at critical points (x_c, y_c) of M(m,x,y)^2 defined by three conditions: ∂ M(m,x_c,y_c)^2/∂ x=0, ∂ M(m,x_c,y_c)^2/∂ y=0, and M(m,x_c,y_c)^2=0.  Since M(m,x,y)^2 is, at most, a second degree polynomial in x and y, solutions to the three conditions can be found analytically.

To see the role of these critical points in practice, we integrate Eq. (<ref>) with respect to y, to reach

    ℐ_𝒜( P_f, P_i, m, 0)  =
    ∫_0^1dx   F^(1)(x)  ,

where
 
    F^(1)(x)    ≡1/(4π)^2 s_iℒ_+ ϵ[y_+(x)] - ℒ_- ϵ[y_-(x)]/y_+(x)-y_-(x) , 
    ℒ_±ϵ[f(x)]    ≡log|1 - x -f(x)/f(x)| + i arctan( 1 - x-Re f(x) /Im f(x)  ±ϵ )  +  i arctan( Re f(x) /Imf(x)  ±ϵ )  ,
 
with y_±(x) ≡ (1/2)(A ±√(A^2+B + i ϵ)). Here A and B are known functions of the kinematic variables and the Feynman parameter x, defined in Eqs. (<ref>) and (<ref>).

We next note that the three Landau conditions given above are satisfied whenever A^2 + B = 0 and, in addition d [A^2 + B]/dx = 0. Noting that A^2 + B is a quadratic polynomial in x we see that the conditions are equivalent to A^2 + B ∝ (x - x_c)^2. Before considering this special case we take our general form and substitute A^2 + B = 2 [(P_i · P_f)^2/s_i^2 - s_f /s_i] (x-x_1)(x-x_2)

    ℐ_𝒜( P_f, P_i, m, 0)  = 1/32 π^2 √(| (P_i · P_f)^2 - s_f  s_i )| [ ∫_0^x_1 + ∫_x_1^x_2  +  ∫_x_2^1    ]  ℒ_+ ϵ[y_+(x)] - ℒ_- ϵ[y_-(x)]/√(   (x-x_1)(x-x_2) ) .

 For concreteness, here we have assumed (P_i · P_f)^2/s_i^2 - s_f /s_i > 0. In addition we have split the integral in x into a sum over regions where Im√(A^2 + B) =0 (x ∈ [0,x_1] ∪ [x_2,1]) as well as the region where it is non-zero (x ∈ [x_1,x_2]). This separation assumes 0 < x_1 < x_2 < 1. If these are instead complex valued, or outside the range of integration, then one can directly evaluate the integral over the entire range. For general values of x_1 and x_2, the integral over each region is well defined even with ϵ = 0 and can be directly evaluated.

The final case to consider is when the external kinematics are tuned to critical values P_f = P_f,c and P_i = P_i,c for which x_1 = x_2 = x_c ∈ [0,1]. This is equivalent to the Landau conditions mentioned above, and corresponds to the apex of the (x - x_1)(x - x_2) parabola sitting on the x-axis. Perturbing the kinematics in one direction shifts the parabola down, opening a finite region of [x_1, x_2 ] that must be integrated in isolation. Perturbing the kinematics in the other direction shifts the parabola upwards, causing the roots to become complex such that we can directly integrate x from 0 to 1. We now demonstrate that, as one approaches P_f,c and P_i,c from the side of real x_1, x_2, the integral ∫_x_1^x_2 dx has a nonzero limit due to a singularity in the integration range. As a consequence, the limit has a different value when approached from opposite sides. This manifests as a step singularity in the real part of ℐ_𝒜. In addition, the imaginary part diverges as log| P_f(i)-P_f(i),c|.

The magnitude of the discontinuity is given by evaluating the integral between x_1 and x_2, for kinematics such that 0 < x_1 < x_2 < 1, and then taking the limit x_1, x_2 → x_c. This can easily be done by noting that, in this region, y_+(x) and y_-(x) are complex conjugates of each other so that the integrand simplifies to

    Disc(ℐ_𝒜) 
       = ∫_x_1^x_2dx F^(1)(x) 
    = 
    1/(4π)^21/s_i∫_x_1^x_2dx 1/iIm y_+
    i(
    arctan[1-x-Re y_+/Im y_+] +
    arctan[Re y_+/Im y_+]
    ) .

Next we note that, as x_1 approaches x_2, Im y_+ goes to zero. Thus it is natural to expand in this quantity

    Disc(ℐ_𝒜)
        = 
    1/(4π)^21/s_i∫_x_1^x_2dx 1/Im y_+(π - Im y_+/1-x-Re y_+ - Im y_+/Re y_+ + 𝒪((Imy_+)^2)
    ) .

We see that only the first term will contribute in the limit x_1 → x_2, and that it only contributes when y_c = y_+(x_c) is in the y integration region, i.e. 0<y_c<1-x_c. Evaluating the remaining integral, we conclude

    Disc(ℐ_𝒜)
       = 1/(4π)^21/s_i∫_x_1^x_2dx π/Im y_+ =
    1/(4π)^21/√((P_i· P_f)^2 -   s_f s_i)∫_x_1^x_2dx π/√((x-x_1)(x_2-x))  ,  
       = 1/16√((P_i· P_f)^2 -   s_f s_i) .








 §.§ Evaluating ℐ_𝒜(P_f, P_i) through ℐ_𝒜; ν_1 ν_2 ν_3 (P_f, P_i) 





Here we provide compact expressions here for the ℐ_𝒜; σ(P_f, P_i) integrals with up to three Lorentz indices. Starting with Eqs. (<ref>) and (<ref>) one can show

    ℐ_𝒜;ν_1(P_f,P_i ,m)  
       =
    P_f,ν_1 ℐ^(1,1)(P_f,P_i ,m)
    + P_i,ν_1ℐ^(1,2)(P_f,P_i ,m)  ,
    
    ℐ_𝒜;ν_1ν_2(P_f,P_i ,m)  
       =
    P_f,ν_1P_f,ν_2ℐ^(2,1)(P_f,P_i ,m)  
    +
    P_i,ν_1P_i,ν_2ℐ^(2,2)(P_f,P_i ,m)  
    
        +
    P_[i,ν_1P_f,ν_2]ℐ^(2,3)(P_f,P_i ,m)  
    -
    g_ν_1ν_2/4ℐ^(2,4)(P_f,P_i ,m)    ,
    
    
    
    ℐ_𝒜;ν_1ν_2ν_3(P_f,P_i ,m)  
       =
    P_f,ν_1P_f,ν_2P_f,ν_3ℐ^(3,1)(P_f,P_i ,m)  
    +
    P_i,ν_1P_i,ν_2P_i,ν_3ℐ^(3,2)(P_f,P_i ,m)  
       
    +
    1/2P_[f,ν_1P_f,ν_2P_i,ν_3]ℐ^(3,3)(P_f,P_i ,m)  
    +
    1/2P_[i,ν_1P_i,ν_2P_f,ν_3]ℐ^(3,4)(P_f,P_i ,m)  
       
    -
    1/8g_[ν_1ν_2P_f,ν_3]ℐ^(3,5)(P_f,P_i ,m)  
    -
    1/8g_[ν_1ν_2P_i,ν_3]ℐ^(3,6)(P_f,P_i ,m)    ,

where the brackets in the indices denote a sum over permutations, even when the indices are identical. The definition is such that, for n indices within a pair of square brackets, the sum runs over n! terms (some of which may vanish). Some examples include

    P_[f,0P_i,0]   =2P_f,0P_i,0,
    
    
    P_[f,0P_f,0P_i,0]   =6P_f,0^2P_i,0,
    
    
    P_[f,1P_i,2]   =P_f,1P_i,2+P_f,2P_i,1,
    g_[00P_i,1]   =
    2g_00P_i,1.


 



Here we have introduced the notation ℐ^(n,m) where the indices just index the integrals needed and do not describe a property of the integrand (i.e. ℐ^(n,m) is just the mth integral needed to evaluate the n-index version of ℐ_𝒜). We now define the set of relevant quantities and also give useful expressions for evaluation. 


To evaluate ℐ_𝒜;ν_1(P_f,P_i ,m) we require

    ℐ^(1,1)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx∫_0^1-xdy
     x   𝒥^0(P_f, P_i, m, 0)  =
     ∫_0^1dx
     x F^(1)(x) ,
    
    
    ℐ^(1,2)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx∫_0^1-xdy
     y   𝒥^0(P_f, P_i, m, 0) 
    =
     ∫_0^1dx
     F^(2)(x).

 
 For ℐ_𝒜;ν_1ν_2(P_f,P_i ,m) we need
 
    ℐ^(2,1)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx∫_0^1-xdy
     x^2   𝒥^0(P_f, P_i, m, 0)
    =
     ∫_0^1dx
     x^2 F^(1)(x) ,
    
    
    ℐ^(2,2)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx∫_0^1-xdy
     y^2   𝒥^0(P_f, P_i, m, 0)
    =
     ∫_0^1dx F^(3)(x) ,
    
    
    
    ℐ^(2,3)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx∫_0^1-xdy
     y x   𝒥^0(P_f, P_i, m, 0)
    =
     ∫_0^1dx
     x F^(2)(x) ,
     
    
    
     ℐ^(2,4)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx∫_0^1-xdy
     𝒥^1(P_f, P_i, m, 0)
    
    
    
    
    
    
    
    
    = ∫_0^1dx
     F^(5)(x)
    +⋯ ,

 where the ellipses denote terms that will be canceled by the Pauli-Villars-like subtractions. These are terms that are independent of m. For example, in the last line above we are ignoring a term proportional to ∫ dx dy y.



  
  Finally, for ℐ_𝒜;ν_1ν_2ν_3(P_f,P_i ,m), six integrals appear
  
    ℐ^(3,1)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx x^3 ∫_0^1-xdy
      𝒥^0(P_f, P_i, m, 0)
     =
     ∫_0^1dx
     x^3 F^(1)(x)
     , 
    ℐ^(3,2)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx∫_0^1-xdy
     y^3   𝒥^0(P_f, P_i, m, 0)
    =  ∫_0^1dx
      F^(4)(x)
    
     , 
    ℐ^(3,3)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx x^2∫_0^1-xdy
     y   𝒥^0(P_f, P_i, m, 0)
    =
     ∫_0^1dx
     x^2 F^(2)(x)
     , 
    ℐ^(3,4)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx x∫_0^1-xdy
     y^2   𝒥^0(P_f, P_i, m, 0)
    =
     ∫_0^1dx
     x F^(3)(x)
     , 
    ℐ^(3,5)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx x∫_0^1-xdy
     𝒥^1(P_f, P_i, m, 0)
    =
     ∫_0^1dx
     x F^(5)(x)
    +⋯
    
     , 
    ℐ^(3,6)(P_f,P_i ,m)
       =Γ(3)∫_0^1dx ∫_0^1-xdy
       y 𝒥^1(P_f, P_i, m, 0)
    =
     ∫_0^1dx
     F^(6)(x)
    +⋯.

 Again, we have ignored terms that cancel after the Pauli-Villars-like subtraction. This includes terms that are proportional to the external momenta but independent of m, e.g.  P_f^ν∫ dx dy x.
 
In the final steps we have evaluated the y-integrals analytically and expressed the remaining x-integral in terms of F^(n)(x). These, in turn, can be written
 
    F^(1)(x) 
       =
     1/(4π)^2 s_iℒ_+ ϵ[y_+(x)] - ℒ_- ϵ[y_-(x)]/y_+(x)-y_-(x) ,
    
     
    
    F^(2)(x) 
       =
     1/(4π)^2 s_iy_+(x)ℒ_+ ϵ[y_+(x)] - y_-(x)ℒ_- ϵ[y_-(x)]/y_+(x)-y_-(x) ,
     
    
     F^(3)(x) 
       
     =
     1/(4π)^2 s_i(   
     y_+(x)^2ℒ_+ ϵ[y_+(x)] - y_-(x)^2ℒ_- ϵ[y_-(x)]/y_+(x)-y_-(x))  
      ,
     
     
    F^(4)(x) 
        
    =
     1/(4π)^2 s_i(   
     (1-x)(y_+(x)+y_-(x))
    + y_+(x)^3ℒ_+ ϵ[y_+(x)] - y_-(x)^3ℒ_- ϵ[y_-(x)] /y_+(x)-y_-(x))
      ,
     
    
     F^(5)(x)
       =
    -1/8π^2(
    (1-x-y_-(x)) log_ϵ(1-x-y_-(x))
    +(1-x-y_+(x)) log_-ϵ(1-x-y_+(x))
    
       +y_-(x) log_ϵ(-y_-(x))
    +y_+(x) log_-ϵ(-y_+(x))
    )  ,
    
    
    
    F^(6)(x)
       
    =
    -1/(4π)^2(
    -(1-x)(y_-(x)+y_+(x))
    +
    ((1-x)^2-y^2_-(x)) log_ϵ(1-x-y_-(x))
    
       
    +
    ((1-x)^2-y^2_+(x)) log_-ϵ(1-x-y_+(x))
    +y^2_-(x) log_ϵ(-y_-(x))
    +y^2_+(x) log_-ϵ(-y_+(x))
    )  ,

where ℒ_±ϵ[f(x)] is defined in Eq. (<ref>). log_ϵ(f(x)) is defined to be the standard log with its branch-cut aligned on the negative real axis. Except if f(x) is purely real and negative, in which case log_±ϵ(f(x))=± i π +log(|f(x)|). As discussed in the main text, these expressions allow one to work identically at ϵ=0 in the numerical evaluation of the integrals with respect to x. The only memory of the non-zero value of ϵ that these functions carry are the sign.




