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

TL;DR

This paper develops a Lorentz-covariant formalism for relating finite-volume matrix elements to two-hadron transition amplitudes, with detailed implementation guidance and an example involving pion scattering and form factors.

Contribution

It introduces a covariant finite-volume formalism for two-hadron states, including a new finite-volume function and detailed numerical evaluation methods.

Findings

Derived a Lorentz-covariant formalism for finite-volume matrix elements.

Introduced a new finite-volume function, G, for numerical evaluation.

Applied the formalism to pion scattering, incorporating form factors and phase shifts.

Abstract

In this work we develop a Lorentz-covariant version of the previously derived formalism for relating finite-volume matrix elements to $\textbf 2 + \mathcal J \to \textbf 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 $\mathcal J$ 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 $\pi \pi + \mathcal J \to \pi \pi$ amplitude in the $\rho$ channel, for which the single-pion form factor, $F_\pi(Q^2)$, as well as the scattering phase, $\delta_{\pi\pi}$, 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.

Tables

Table 1. Table 1: Summary of the conditions under which 𝒵 J ​ M ( n ) = 0 subscript superscript 𝒵 𝑛 𝐽 𝑀 0 \mathcal{Z}^{(n)}_{JM}=0 .

P	$\left(n\right)$	$m_1,m_2$	${𝒵}_{JM}^{\left(n\right)}=0$	comments
$\left[000\right]$	all	general	for all $M\notin 4\mathbb{Z}$, all $J\notin 2\mathbb{Z}$, $JM=20$	also ${𝒵}_{J,M}^{\left(n\right)}={𝒵}_{J,-M}^{\left(n\right)}$
$\left[00d_z\right]$	all	general	for all $M\notin 4\mathbb{Z}$	also ${𝒵}_{J,M}^{\left(n\right)}={𝒵}_{J,-M}^{\left(n\right)}$
$\left[0d_y{d}_z\right]$	all	general		${𝒵}_{J,M}^{\left(n\right)}={𝒵}_{J,-M}^{\left(n\right)}$
general	1	$m_1=m_2$	for all $J\notin 2\mathbb{Z}$	only up to $𝒪\left(e^{-mL}\right)$

Table 2. Table 2: Values of c j subscript 𝑐 𝑗 c_{j} up to n = 5 𝑛 5 n=5 , assuming Λ j = 2 j − 1 ​ ( 3 ​ m ) subscript Λ 𝑗 superscript 2 𝑗 1 3 𝑚 \Lambda_{j}=2^{j-1}(3m) . This information is sufficient to calculate G [ μ 1 ​ ⋯ ​ μ n ; ℓ f ​ m f ; ℓ i ​ m 1 ] subscript 𝐺 subscript 𝜇 1 ⋯ subscript 𝜇 𝑛 subscript ℓ 𝑓 subscript 𝑚 𝑓 subscript ℓ 𝑖 subscript 𝑚 1 G_{[\mu_{1}\cdots\mu_{n};\ell_{f}m_{f};\ell_{i}m_{1}]} for all indices satisfying n + ℓ f + ℓ i ≤ 10 𝑛 subscript ℓ 𝑓 subscript ℓ 𝑖 10 n+\ell_{f}+\ell_{i}\leq 10 .

$n$	$c_1$	$c_2$	$c_3$	$c_4$	$c_5$
$1$	$-1$
$2$	$-\frac{35}{27}$	$\frac{8}{27}$
$3$	$-\frac{1001}{729}$	$\frac{286}{729}$	$-\frac{14}{729}$
$4$	$-\frac{82225}{59049}$	$\frac{16445}{39366}$	$-\frac{4025}{157464}$	$\frac{143}{472392}$
$5$	$-\frac{37872835}{27103491}$	$\frac{5410405}{12754584}$	$-\frac{1853915}{68024448}$	$\frac{329329}{816293376}$	$-\frac{16445}{13876987392}$

Equations

$$E_i^{\star }\equiv \sqrt{P_i^2}=\sqrt{E_i^2-{\text{P}}_i^2},\text{and}{E}_f^{\star }\equiv \sqrt{P_f^2}=\sqrt{E_f^2-{\text{P}}_f^2}.$$

$$Q^2\equiv -(P_f-P_i{)}^2,$$

$${\omega }_{k2}\equiv \sqrt{k^2+m_2^2},$$

$${\omega }_{Pk1}\equiv \sqrt{(\text{P}-\text{k}{)}^2+m_1^2}.$$

$$k^{\star \mu }\equiv ({\omega }_{k2}^{\star },{\text{k}}^{\star })\equiv {[{\Lambda }_{\mathbf{\beta }}{]}^{\mu }}_{\nu }{k}^{\nu },$$

$${\omega }_{k2}^{\star }=\sqrt{k^{\star 2}+m_2^2},$$

$$P^{\star \mu }-k^{\star \mu }=(E^{\star }-{\omega }_{k2}^{\star },-{\text{k}}^{\star })\equiv {[{\Lambda }_{\mathbf{\beta }}{]}^{\mu }}_{\nu }(P^{\nu }-k^{\nu }),$$

$${\omega }_{k1}^{\star }\equiv \sqrt{k^{\star 2}+m_1^2}.$$

$$E^{\star }\equiv \sqrt{q^{\star 2}+m_1^2}+\sqrt{q^{\star 2}+m_2^2}.$$

$$\det [\mathcal{M}(P^2{)}^{-1}+F(P,L)]=0.$$

$$\mathcal{M}_{\ell m;\ell^{\prime}m^{\prime}}(P^{2})=\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\mathcal{M}^{(\ell)}(P^{2})\,,\ \ \ \text{with}\ \ \ \mathcal{M}^{(\ell)}(s)=\frac{8\pi E^{\star}}{\xi}\frac{1}{q^{\star}\cot\delta^{(\ell)}(q^{\star})-iq^{\star}}\bigg{|}_{E^{\star 2}=s}\,,$$

$$F_{\ell m;\ell^{\prime}m^{\prime}}(P,L)\equiv\xi\bigg{[}\frac{1}{L^{3}}\sum_{\mathbf{k}}-\int\!\!\frac{d^{3}\textbf{k}}{(2\pi)^{3}}\,\bigg{]}\frac{\mathcal{Y}_{\ell m}(\textbf{k}^{\star})\,\,\mathcal{Y}^{*}_{\ell^{\prime}m^{\prime}}(\textbf{k}^{\star})}{2\omega_{Pk1}2\omega_{k2}(E-\omega_{k2}-\omega_{Pk1}+i\epsilon)},$$

$${\mathcal{Y}}_{\ell m}({\text{k}}^{\star })$$

$$\mathcal{R}(E_n,\text{P})\equiv \underset{E\to E_n}{\lim }\left[\frac{(E-E_n)}{F^{-1}(P,L)+\mathcal{M}(P^2)}\right].$$

$${\mathcal{W}}_{{\mu }_1\cdots {\mu }_n}(P_f,k^{\prime };P_i,k)\equiv ⟨P_f,k^{\prime };\text{out}|{\mathcal{J}}_{{\mu }_1\cdots {\mu }_n}(0)|P_i,k;\text{in}{⟩}_{conn.}.$$

$${\mathcal{W}}_{\mathrm{df};{\mu }_1\cdots {\mu }_n}\equiv {\mathcal{W}}_{{\mu }_1\cdots {\mu }_n}-i\overline{\mathcal{M}}(P_f,k^{\prime },k)\frac{i}{(P_f-k{)}^2-m_1^2}{w}_{{\mu }_1\cdots {\mu }_n}-w_{{\mu }_1\cdots {\mu }_n}\frac{i}{(P_i-k^{\prime }{)}^2-m_1^2}i\overline{{\mathcal{M}}^{\prime }}(P_i,k^{\prime },k),$$

$$\overline{\mathcal{M}}(P_f,k^{\prime },k)$$

$$\overline{{\mathcal{M}}^{\prime }}(P_i,k^{\prime },k)$$

$${\mathcal{W}}_{{\rm df};\mu_{1}\cdots\mu_{n}}(P_{f},k^{\prime};P_{i},k)\equiv 4\pi Y^{*}_{\ell^{\prime}m^{\prime}}\big{(}\hat{\textbf{k}}^{\prime\star}_{f}\big{)}\ {\mathcal{W}}_{{\rm df};\mu_{1}\cdots\mu_{n};\ell^{\prime}m^{\prime};\ell m}(s_{f},s_{i},Q^{2})\ Y_{\ell m}\big{(}\hat{\textbf{k}}^{\star}_{i}\big{)}\,,$$

$$w_{{\mu }_1\cdots {\mu }_n}^{\text{on}}(P_f-k,P_i-k)\equiv ⟨P_f-k,m_1|{\mathcal{J}}_{{\mu }_1\cdots {\mu }_n}(0)|P_i-k,m_1⟩.$$

$$w_{{\mu }_1\cdots {\mu }_n}^{\text{off}}(P_f-k,P_i-k)$$

$$w_{{\mu }_1\cdots {\mu }_n}(P_f-k,P_i-k)$$

$${\text{K}}_{{\mu }_1\cdots {\mu }_n}^{(j)}(m,\text{k},P_f,P_i)\equiv \sum _{n^{\prime }=0}^n{K}_{{\mu }_1\cdots {\mu }_{n^{\prime }}}^{\omega }(m,\text{k})C_{{\mu }_{n^{\prime }+1}\dots {\mu }_n}^{(j)}(P_f,P_i),$$

$$K^{\omega}_{\mu_{1}\cdots\mu_{n}}(m,\textbf{k})\equiv k_{\mu_{1}}\cdots k_{\mu_{n}}\bigg{|}_{k_{0}=\omega_{k2}}\,.$$

$$\Big{|}\langle P_{f},L|\mathcal{J}^{\mu_{1}\cdots\mu_{n}}(0)|P_{i},L\rangle\Big{|}^{2}=\frac{1}{L^{6}}{\rm{Tr}}\left[\mathcal{R}(P_{i}){\mathcal{W}}_{L,{\rm df}}^{\mu_{1}\cdots\mu_{n}}(P_{i},P_{f},L)\mathcal{R}(P_{f}){\mathcal{W}}_{L,{\rm df}}^{\mu_{1}\cdots\mu_{n}}(P_{f},P_{i},L)\right]\,,$$

$$\begin{gathered} {\mathcal{W}}_{L,\mathrm{df}}^{{\mu }_1\cdots {\mu }_n}(P_f,P_i,L)-{\mathcal{W}}_{\mathrm{df}}^{{\mu }_1\cdots {\mu }_n}(s_f,s_i,Q^2)\equiv \sum _j\sum _{n^{\prime }=0}^n{C}^{(j),{\mu }_{n^{\prime }+1}\dots {\mu }_n}(P_f,P_i)f^{(j)}(Q^2) \\ \times \mathcal{M}(s_f)G^{{\mu }_1\cdots {\mu }_{n^{\prime }}}(P_f,P_i,L)\mathcal{M}(s_i). \end{gathered}$$

$$G_{{\mu }_1\cdots {\mu }_n;{\ell }_f{m}_f;{\ell }_i{m}_i}(P_f,P_i,L)$$

$$D(m,\text{k})$$

$${\mathcal{W}}_{L,\mathrm{df};{\Lambda }_f{\mu }_f;{\Lambda }_i{\mu }_i}^{\mu }(P_f,P_i,L)=e^{i{\delta }_{\pi \pi }^{I=1}(q_f^{\star })}\frac{⟨E_{f,n^{\prime }},{\text{P}}_f,L;{\Lambda }_f{\mu }_f|j^{\mu }(0)|E_{i,n},{\text{P}}_i,L;{\Lambda }_i{\mu }_i⟩}{|{\mathcal{R}}_{{\Lambda }_f{\mu }_f}^{I=1}(P_f,L){\mathcal{R}}_{{\Lambda }_i{\mu }_i}^{I=1}(P_i,L){|}^{1/2}}{e}^{i{\delta }_{\pi \pi }^{I=1}(q_i^{\star })},$$

$$|\mathcal{R}^{I=1}_{\Lambda\mu}(P,L)|^{1/2}\equiv\frac{1}{L^{3}}\sqrt{\frac{q^{\star}}{16\pi E^{\star}}}\left[\frac{E^{\star}}{E}\frac{\partial}{\partial E^{\star}}\Big{(}\phi_{\Lambda\mu}^{\textbf{d}}(E^{\star},L)+\delta^{I=1}_{\pi\pi}(q^{\star})\Big{)}\right]^{-1/2}\,,$$

$${\mathcal{W}}_{\mathrm{df};{\Lambda }_f{\mu }_f;{\Lambda }_i{\mu }_i}^{\mu }(s_f,s_i,Q^2)$$

Full text

finite volume, lattice QCD

††preprint:

JLAB-THY-18-2878 ††preprint:

CERN-TH-2018-263

I Introduction

In recent years, interest in hadron spectroscopy has increased significantly, primarily due to various experimental discoveries of unconventional excitations.111For recent reviews of the experimental and theoretical understanding of these states we point the reader to Refs. Guo et al. (2018); Lebed et al. (2017); Liu (2014); Briceño et al. (2016); Chen et al. (2016). 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 Guo et al. (2018); Lebed et al. (2017); Liu (2014); Briceño et al. (2016); Chen et al. (2016). 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 $\textbf{2}+\mathcal{J}\to\textbf{2}$, where each bold-faced 2 counts the QCD-stable hadrons in the state and $\mathcal{J}$ represents a generic, local external current.

The approach discussed here is based on prior formalism developed by two of us in Ref. Briceño and Hansen (2016).222This, in turn, was inspired and guided by the work of Refs. Briceño and Davoudi (2013a); Bernard et al. (2012); Detmold and Savage (2004). In Sec. II 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 $\textbf{1}+\mathcal{J}\to\textbf{1}$, are expressed in terms of standard Lorentz-invariant form factors. After extracting the $\textbf{2}+\mathcal{J}\to\textbf{2}$ transition amplitudes, one can proceed to determine form factors as well as distribution functions333Distribution functions are accessed in lattice calculations via spatially non-local operator insertions Ji (2013). These may suffer from further finite-volume effects associated with the size of the operators as discussed in Ref. Briceño et al. (2018). 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 Winter et al. (2017); Detmold and Shanahan (2016) 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. Briceño et al. (2018a) for a recent review.

The most well-established such relation is that derived by Lüscher in Refs. Lüscher (1986, 1991) 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 Rummukainen and Gottlieb (1995); Feng et al. (2004); He et al. (2005); Bedaque (2004); Liu et al. (2006); Kim et al. (2005); Christ et al. (2005); Lage et al. (2009); Bernard et al. (2011); Fu (2012); Leskovec and Prelovsek (2012); Briceño and Davoudi (2013a); Hansen and Sharpe (2012); Guo et al. (2013); Li and Liu (2013); Briceño et al. (2014); Briceño (2014). These formal ideas, together with significant algorithmic developments, have resulted in a proliferation of scattering amplitudes determined directly from lattice QCD Dudek et al. (2013); Wilson et al. (2015a); Bali et al. (2016); Bulava et al. (2016); Andersen et al. (2018); Alexandrou et al. (2017); Brett et al. (2018); Moir et al. (2016); Woss et al. (2018); Dudek et al. (2014); Wilson et al. (2015b); Dudek et al. (2016); Briceño et al. (2018b, 2017). 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 Briceño and Davoudi (2013b); Hansen and Sharpe (2014); Mai and Döring (2017); Hammer et al. (2017); Briceño et al. (2017, 2018c).

Similar developments have resulted in both perturbative Detmold and Flynn (2015) non-perturbative Lellouch and Lüscher (2001); Kim et al. (2005); Christ et al. (2005); Hansen and Sharpe (2012); Meyer (2011); Feng et al. (2015); Meyer (2012); Agadjanov et al. (2014); Briceño and Davoudi (2013a); Bernard et al. (2012); Briceño and Hansen (2015); Briceño et al. (2015a); Briceño and Hansen (2016) relations between finite-volume matrix elements and electroweak amplitudes. These have already been implemented in a variety of LQCD calculations Bai et al. (2015); Blum et al. (2012a); Boyle et al. (2013); Blum et al. (2011, 2012b); Feng et al. (2015); Briceño et al. (2015b, 2016); Alexandrou et al. (2018). In particular, Refs. Briceño et al. (2015b, 2016); Alexandrou et al. (2018) extracted the $\rho\gamma^{\star}\to\pi$ electromagnetic form factor by determining the energy dependence of the corresponding amplitude, $\pi\pi\gamma^{\star}\to\pi$.

In this work we turn our attention to the prospect of determining $\textbf{2}+\mathcal{J}\to\textbf{2}$ transition amplitudes from finite-volume matrix elements. This was previously considered in Ref. Briceño and Hansen (2016). 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. Briceño and Hansen (2016), 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. Briceño and Hansen (2016) in this work we derive a mapping between finite-volume matrix elements of two-particle states and the infinite-volume $\textbf{2}+\mathcal{J}\to\textbf{2}$ amplitude. The result is summarized by the flow-chart shown in Fig. 1. We find that, given the following quantities:

• •

the two-particle finite-volume spectrum,

• •

the $\textbf{1}+\mathcal{J}\to\textbf{1}$ form factors,

• •

the finite-volume two-particle matrix elements of $\mathcal{J}$,

one can systematically constrain the $\textbf{2}+\mathcal{J}\to\textbf{2}$ amplitude in the kinematic window in which only the accommodated channel contributes. Our relation requires the generalized Lellouch-Lüscher factors Lellouch and Lüscher (2001); Briceño et al. (2015a); Briceño and Hansen (2016), 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. Detmold and Flynn (2015) showed that, in the case of weak interactions, the finite-volume effects on the ground-state $\textbf{2}+\mathcal{J}\to\textbf{2}$ matrix element are given by an expansion in powers of $1/L$.444In fact, the authors of Ref. Detmold and Flynn (2015) consider $\textbf{n}+\mathcal{J}\to\textbf{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. II, where we also discuss three key examples involving $\pi\pi$ states. Then, in Sec. III, 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. IV we conclude and provide an outlook for future studies. Finally we include two appendices: Appendix A gives various details relevant for the derivation of the improved formalism used here. Appendix B includes various technical aspects regarding the evaluation of the finite-volume functions discussed in the main text.

II Covariant representation of the formalism

In this section we revisit the formalism derived in Ref. Briceño and Hansen (2016) 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 A we explain the (minimal) modifications to the original derivation that lead to this new form.

This section is divided into three subsections: In Sec. II.1 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. II.2, we provide a full description of our covariant $\textbf{2}+\mathcal{J}\to\textbf{2}$ formalism. Finally, in Sec. II.3, we consider a handful of specific examples to show how our general expressions reduce for a particular system with specified quantum numbers.

II.1 Kinematics and the quantization condition

We denote the 4-momentum of the incoming state in the finite-volume frame by $P^{\mu}_{i}\equiv(E_{i},\textbf{P}_{i})$ and that of the outgoing state by $P^{\mu}_{f}\equiv(E_{f},\textbf{P}_{f})$. The center-of-momentum (c.m.) energies corresponding to these are then given by

[TABLE]

This also defines the metric used for the Minkowski-signature 4-vector dot products throughout. Following the notation of Ref. Kim et al. (2005) we use $\star$ 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 $\star$, 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. $\textbf{P}_{i}=2\pi\textbf{d}_{i}/L$ and $\textbf{P}_{\!f}=2\pi\textbf{d}_{f}/L$ where $\textbf{d}_{i}$ and $\textbf{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}^{\mu}-P_{i}^{\mu}$ [see Fig. 2(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

[TABLE]

which is positive for spacelike $P_{f}^{\mu}-P_{i}^{\mu}$.

As mentioned above, we restrict our attention here to values of $E_{i}^{\star}$ and $E_{f}^{\star}$ 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, $\mathcal{J}$, is flavor conserving so that the same two particles appear in the initial and final state.555Given the results presented below, implementing the covariant modification to the multi-channel expressions of Ref. Briceño and Hansen (2016) 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^{\mu}\equiv(E,\textbf{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^{\mu}=(\omega_{k2},\textbf{k})$, where

[TABLE]

with $k=|\textbf{k}|$.

Next note that, in order to satisfy the specified total 4-momentum ($P^{\mu}$), particle 1 must carry $P^{\mu}-k^{\mu}\equiv(E-\omega_{k2},\textbf{P}-\textbf{k})$. Thus, for general k, one cannot simultaneously require that the particle momenta sum to $P^{\mu}$ and that particle 1 is on shell. The latter holds only when the temporal component, $E-\omega_{k2}$, is equal to

[TABLE]

To better understand when the on-shell condition ($E-\omega_{k2}=\omega_{Pk1}$) is satisfied, it is useful to introduce ${[\Lambda_{\boldsymbol{\beta}}]^{\mu}}_{\nu}\equiv{\Lambda^{\mu}}_{\nu}(\boldsymbol{\beta})$ as a boost matrix with boost velocity $\boldsymbol{\beta}\equiv\textbf{P}/E$. We then define

[TABLE]

and observe

[TABLE]

where $k^{\star}\equiv|\textbf{k}^{\star}|$.666It is worth emphasizing that the definitions of $\textbf{k}^{\star}$ and $k^{\star}$ depend on $(E,\textbf{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

[TABLE]

where we have used ${[\Lambda_{\boldsymbol{\beta}}]^{\mu}}_{\nu}P^{\nu}=(E^{\star},\textbf{0})$. We deduce that the c.m. frame on-shell condition is $E^{\star}=\omega_{k1}^{\star}+\omega_{k2}^{\star}$ where

[TABLE]

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

[TABLE]

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^{\mu}\to P_{i}^{\mu}$, then the corresponding quantities above become $\omega_{P_{i}k2}$, ${[\Lambda_{\boldsymbol{\beta}_{i}}]^{\mu}}_{\nu}$, $k^{\star\mu}_{i}$, $\textbf{k}^{\star}_{i}$, $k^{\star}_{i}$, $\omega_{k1i}^{\star}$, $\omega_{k2i}^{\star}$, $E^{\star}_{i}$, $q^{\star}_{i}$. The only quantities that do not inherit a frame index are the finite-volume frame momentum, k, and the corresponding energy, $\omega_{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, $\mathcal{M}(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  Kim et al. (2005); Hansen and Sharpe (2012); Briceño (2014)

[TABLE]

This holds in the region $0<E^{\star}<E^{\star}_{\text{th}}$, where $E^{\star}_{\text{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^{-mL}$, where $m$ is the physical mass of the lightest degree of freedom in the theory.

The precise definition of $\mathcal{M}(P^{2})$ is given by

[TABLE]

where $\delta^{(\ell)}(q^{\star})$ is the scattering phase shift, and $\xi=1/2$ if the particles are identical and $\xi=1$ otherwise.

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

[TABLE]

where

[TABLE]

In Appendix B.2 we review an efficient method to evaluate this, based on analytic expressions for the integrals defined using the cutoff functions introduced in Ref. Kim et al. (2005).

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

[TABLE]

This object, introduced in this form in Ref. Briceño et al. (2015a), is the generalization of the Lellouch-Lüscher factor Lellouch and Lüscher (2001) that relates finite-volume matrix elements of two-particle states to the corresponding infinite-volume decay and transition amplitudes.

II.2 Relating finite-volume matrix elements with $\textbf{2}+\mathcal{J}\to\textbf{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. Briceño and Hansen (2016) in two key ways.

First, the separation of singularities, required to disentangle finite-volume effects in the $\textbf{2}+\mathcal{J}\to\textbf{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\omega_{k}(k^{0}-\omega_{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. 2(b) and defined explicitly in Fig. 2(d) and in Eq. (27) below.

Second, we treat the single-particle matrix elements in a simpler way here than we did in Ref. Briceño and Hansen (2016). Our approach requires isolating the $\textbf{1}+\mathcal{J}\to\textbf{1}$ matrix element in order to express the finite-volume effects of the triangle diagram [Fig. 2(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 $\textbf{2}+\mathcal{J}\to\textbf{2}$ matrix element that we are after [see also Fig. 2(c)]

[TABLE]

Here the initial and final states are standard two-particle asymptotic states with the usual relativistic normalization convention and $\mathcal{J}_{\mu_{1}\cdots\mu_{n}}(0)$ is a generic local current insertion. As is shown in Fig. 2(a), the initial state is built from particle 1 [with mass $m_{1}$ and on-shell momentum $(P_{i}-k)^{\mu}$ 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^{\prime}$ and $k^{\prime}$. Following the discussion of the previous sub-section, the c.m. frame 3-momenta are denoted by $\textbf{k}_{i}^{\star}$ and $\textbf{k}_{f}^{\star}$ and have magnitudes equal to $q_{i}^{\star}$ and $q_{f}^{\star}$ respectively. The label “conn.” emphasizes that only fully connected diagrams, those shown in Fig. 2(c), contribute to the definition of the $\textbf{2}+\mathcal{J}\to\textbf{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, $\hat{\textbf{k}}^{\star}_{i}$ and $\hat{\textbf{k}}^{\star}_{f}$. However, in contrast to the scattering amplitude, $\mathcal{M}$, for $\mathcal{W}$ 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. 2(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. Briceño and Hansen (2016) where the quantity $\mathcal{W}_{\mathrm{df}}$ 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. 2(c)]777This 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. (111) as well as Ref. Briceño and Hansen (2016) for explicit expressions.

[TABLE]

where $w_{\mu_{1}\cdots\mu_{n}}$ is the single-particle matrix element of the external current, defined in detail below, and

[TABLE]

Here $P_{\ell}(\cos\theta)$ are the standard Legendre polynomials, satisfying $\sum_{m=-\ell}^{\ell}4\pi Y^{*}_{\ell m}(\hat{\textbf{x}})Y_{\ell m}(\hat{\textbf{y}})=(2\ell+1)P_{\ell}(\hat{\textbf{x}}\cdot\hat{\textbf{y}})$. Unlike ${\mathcal{W}}$, ${\mathcal{W}}_{{\rm df}}$ admits a uniformly convergent decomposition in spherical harmonics

[TABLE]

where $s_{i}\equiv P_{i}^{2}=E_{i}^{\star 2}$, $s_{f}\equiv P_{f}^{2}=E_{f}^{\star 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 $\textbf{2}\to\textbf{2}$ and $\textbf{1}+\mathcal{J}\to\textbf{1}$ transitions. In particular, as we show in Sec. III.2.2 below, ${\mathcal{W}}_{{\rm df}}$ does have so-called triangle singularities as a function of $E_{i}^{\star}$ and $E_{f}^{\star}$ associated with diagram of the kind depicted in Fig. 2(b). These pose no problem to the harmonic decomposition but must, of course, be understood in order to successfully extract and interpret both ${\mathcal{W}}_{{\rm df}}$ and ${\mathcal{W}}$.

The scattering amplitude, $\overline{\mathcal{M}}$, is defined with powers of $(k^{\star}/q^{\star})$, referred to as barrier factors. These must be included due to the manner in which the factors of $w_{\mu_{1}\cdots\mu_{n}}$ and ${\mathcal{M}}$ arise in the triangle diagram of Fig. 2(b). In particular, the loop is summed over all finite-volume momenta in such a way that the current insertion, $w_{\mu_{1}\cdots\mu_{n}}$, as well as the external factors of $\mathcal{M}$, are sampled at off-shell values of momenta (i.e. values for which $k^{2}\neq m_{i}^{2}$). Nonetheless, the power-law finite-volume effects are governed by the on-shell values of $w_{\mu_{1}\cdots\mu_{n}}$ and ${\mathcal{M}}$ 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^{\star}\to q^{\star}$, amounts to replacing $k^{\star\ell}Y_{\ell m}(\hat{{\textbf{k}}^{\star}})$ with $q^{\star\ell}Y_{\ell m}(\hat{{\textbf{k}}^{\star}})$. The latter has spurious singularities near $\textbf{k}^{\star}=\textbf{0}$ and thus more care is needed. The inclusion of barrier factors resolves the issue.

The single-particle matrix element, $w_{\mu_{1}\cdots\mu_{n}}$, appearing in Eq. (16), is a function of $(P_{f}-k,P_{i}-k)$ in the first term and $(P_{f}-k^{\prime},P_{i}-k^{\prime})$ 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

[TABLE]

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

[TABLE]

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 $w_{\mu_{1}\cdots\mu_{n}}$ that appears in Eq. (16). In this step we set the virtualities within the form factors to be on shell [$(P_{f}-k)^{2},(P_{i}-k)^{2}\longrightarrow m_{1}^{2}$] and also set $k^{0}=\omega_{k2}$ everywhere. We reach

[TABLE]

This definition is not completely on shell because, within $\textbf{K}^{(j)}_{\mu_{1}\cdots\mu_{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 $\textbf{K}^{(j)}_{\mu_{1}\cdots\mu_{n}}$ as a sum of terms that isolate the dependence on this quantity

[TABLE]

where $C^{(j)}_{\mu_{n^{\prime}+1}\ldots\mu_{n}}(P_{f},P_{i})$ has no indices for $n^{\prime}=n$ and

[TABLE]

Here the superscript $\omega$ indicates that the 4-momenta in $K^{\omega}$ are on shell. This will be important in Sec. III 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 $\mathcal{J}$ is sandwiched between two finite-volume states, each of which has the quantum numbers of the two-particle system. As we demonstrate in Appendix A, following the derivation of Ref. Briceño and Hansen (2016), this LQCD observable is related to the infinite-volume $\textbf{2}+\mathcal{J}\to\textbf{2}$ transition amplitude via

[TABLE]

where ${\mathcal{W}}_{L,{\rm df}}$ directly determines ${\mathcal{W}}_{\rm df}$ via

[TABLE]

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. (114) as well as Ref. Briceño and Hansen (2016) for explicit expressions.

In these expressions we have suppressed angular-momentum indices on $\mathcal{R}$, ${\mathcal{W}}_{L,{\rm df}}^{\mu_{1}\cdots\mu_{n}}$, ${\mathcal{W}}_{\rm df}^{\mu_{1}\cdots\mu_{n}}$, $\mathcal{M}$ and $G^{\mu_{1}\cdots\mu_{n^{\prime}}}$. Each object carries the set $\ell^{\prime}m^{\prime};\ell m$ and these are contracted between adjacent factors in the usual matrix multiplication. The trace in Eq. (25) is also over this angular-momentum space. Note, by contrast, that the index set $\mu_{1}\cdots\mu_{n}$ is not summed but rather fixed to common values for all objects appearing in these equations.

Finally, the matrix $G_{\mu_{1}\cdots\mu_{n^{\prime}}}(P_{f},P_{i},L)$ is defined diagrammatically in Fig. 2(d) and has the explicit form

[TABLE]

where

[TABLE]

This differs from the form presented in Ref. Briceño and Hansen (2016) due to the aforementioned modifications: The poles are Lorentz invariant and the $\textbf{1}+\mathcal{J}\to\textbf{1}$ matrix element is expressed in terms of tensor structures leading to $K^{\omega}$. Note that the modifications to $G$ are directly connected to those in the definition of ${\mathcal{W}}_{\rm df}$, Eq. (16). We have altered these two intermediate objects in such a way that $\mathcal{W}$ is unchanged.

II.3 Examples

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

II.3.1 $(\pi^{+}\pi^{0})_{J=1}+j_{\mu}\to(\pi^{+}\pi^{0})_{J=1}$

We begin with the electro-magnetic form factors of a charged $\rho$ meson. The $\rho$ decays predominantly to the vector-isovector $\pi\pi$ state. Indeed for heavier than physical light-quark masses (such that $4M_{\pi}>M_{\rho}$) 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 $\rho$ mass ($2M_{\pi}>M_{\rho}$) and the latter becomes a stable particle. In this case, one can extract the form factors of the $\rho$ directly from finite-volume matrix elements. See, for example, Ref. Shultz et al. (2015).

To determine the analogous observable at quark masses for which the $\rho\to\pi\pi$ decay occurs, it is necessary to first extract the $(\pi^{+}\pi^{0})_{J=1}+j_{\mu}\to(\pi^{+}\pi^{0})_{J=1}$ amplitude for a wide range of kinematic points. As depicted in Fig. 1, 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 $\rho$ form factors. Detailing the steps of this continuation will be the focus of future work. Here we focus on the extraction of the $(\pi^{+}\pi^{0})_{J=1}+j_{\mu}\to(\pi^{+}\pi^{0})_{J=1}$ amplitude for real $\pi\pi$ energies.

By interpolating isospin-one initial and final states ($I=1,m_{I}=\pm 1$), we project the system to a sector where all even angular momenta vanish. Thus, regardless of the values of $\textbf{P}_{f}$ and $\textbf{P}_{i}$, we will always have a finite-volume irrep that contains $J^{P}=1^{-}$, with the next contamination coming from $J\geq 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 Thomas et al. (2012).

Considering first the incoming state, we restrict attention to a specified set of finite-volume quantum numbers: $\Lambda_{i}$ and $\mu_{i}$, labeling the irrep and row respectively, of the little group defined by $\textbf{P}_{i}$. In addition we assume that, within this irrep, the interactions are dominated by $\ell_{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. (25) becomes888Note, the procedure for subducing the matrix elements onto the appropriate symmetry group is discussed in detail in Ref. Briceño and Hansen (2016). Although some of the details of the formalism has changed, this aspect remains the same.

[TABLE]

where $\delta^{I=1}_{\pi\pi}$ is the elastic scattering phase shift for $\pi\pi\to\pi\pi$ in the $\rho$ channel and

[TABLE]

with $\phi_{\Lambda\mu}^{\textbf{d}}(E^{\star},L)$ defined for the irreps that couple to $\ell=1$ $\pi\pi$ in, for example, Ref. Briceño et al. (2015a). 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^{\prime}}(L)$ and $E_{i,n}(L)$ respectively, satisfying Eq. (10). In Eq. (29), $\mathcal{W}_{L,{\rm df};\Lambda_{f}\mu_{f};\Lambda_{i}\mu_{i}}^{\mu}$ is the subduced version of $\mathcal{W}_{L,{\rm df}}^{\mu}$. As discussed in Ref. Briceño and Hansen (2016), this can be obtained from $\mathcal{W}_{L,{\rm df};\ell_{f}m_{f};\ell_{i}m_{i}}^{\mu}$ by rotating these into the helicity basis and then using the subduction matrices Thomas et al. (2012). 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. (25), in Eq. (29) we have dropped the trace, since we are ignoring all but one partial wave, and have solved for ${\mathcal{W}}_{L,{\rm df}}^{\mu}$. For the latter step there is a potential sign ambiguity that one must address. Note that $\mathcal{R}^{I=1}_{\Lambda\mu}(P,L)=|\mathcal{R}^{I=1}_{\Lambda\mu}(P,L)|e^{-2i\delta^{I=1}_{\pi\pi}(q^{\star})}$ as is shown, for example, in Eqs. (132)-(134) of Ref. Briceño and Hansen (2015). The phases in $\mathcal{R}^{I=1}_{\Lambda\mu}(P_{f},L)$ and $\mathcal{R}^{I=1}_{\Lambda\mu}(P_{i},L)$ precisely generate the Watson phases within ${\mathcal{W}}_{L,{\rm df}}^{\mu}(P_{f},P_{i},L)$ as they must, since the finite-volume matrix element is real. The remaining $\pm$ ambiguity is constrained by the known value in the $Q^{2}\to 0$ limit, but in certain cases a remaining ambiguity may survive.

The final step is to convert this to the infinite-volume quantity ${\mathcal{W}}_{\rm df}^{\mu}$ via

[TABLE]

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

[TABLE]

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}\neq 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 $\mu$-indexed equations by a common factor.

II.3.2 $(\pi^{+}\pi^{0})_{I_{i}}+j_{\mu}\to(\pi^{+}\pi^{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 $\pi\pi$ states must vanish. If we take the angular momentum to be unconstrained then the incoming $\pi\pi$ state may, in general carry isospin $I_{i}=0,1,2$. The iso-vector current then couples the fixed incoming isospin as follows: $0\to 1$, $1\to 0,1,2$ and $2\to 1,2$. If we further restrict attention to $|m_{I}|=1$ states (i.e. $\pi^{+}\pi^{0}$ states) then this reduces to $1\to 1,2$ and $2\to 1,2$. An alternative way to distinguish these possibilities is by fixing orbital angular momentum: $J=0\Leftrightarrow I=2$ and $J=1\Leftrightarrow I=1$.

In this way we identify four possible transitions involving $\pi^{+}\pi^{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. (29) to (II.3.1) provided that we can neglect the effects of finite-volume mixing with $J\geq 2$ states. The only modifications are that one must generally project to different irreps $\Lambda$ to access the $J=0$ components and that the scattering amplitudes must correspond to the isospin of the state

[TABLE]

with

[TABLE]

II.3.3 Gluonic structure

Thirty years ago, Jaffe and Manohar identified a structure function that provides a measure of the gluon distribution within hadrons Jaffe and Manohar (1989). This has since lead to lattice calculations of the leading moments of these distributions, for example within the $\phi$-meson. Thus far, the calculations are restricted to heavy quark masses where the $\phi$ is stable within QCD Detmold et al. (2017); Detmold and Shanahan (2016). 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 Winter et al. (2017). 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 [$\mathcal{O}_{\mu\nu}\sim\mathcal{G}_{\mu\beta}\,{\mathcal{G}}_{\beta\nu}$]. 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 $\rho$, discussed in Sec. II.3.1 above. In this case, the relation between the finite-volume matrix elements and the transition amplitude is very similar to Eq. (29). The only distinction arises in relating $\mathcal{W}^{\mu\nu}_{{\rm df}}$ and $\mathcal{W}^{\mu\nu}_{L,{\rm df}}$. To do so one must determine the scalar ($G$), vector ($G^{\mu}$), and tensor ($G^{\mu\nu}$) 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) $\textbf{2}\to\textbf{2}$ scattering amplitudes and (iii) finite-volume kinematic functions, can be combined with finite-volume two-particle matrix elements to extract the infinite-volume $\textbf{2}+\mathcal{J}\to\textbf{2}$ transition amplitudes. In this recipe the ingredient that remains most obscure is the new finite-volume function $G^{\mu_{1}\cdots\mu_{n}}(P_{f},P_{i},L)$, defined in Eq. (27). Thus, in the next section, we give a detailed description of how this can be efficiently evaluated numerically.

III Evaluating $G(P_{f},P_{i},L)$

Our aim is to evaluate

[TABLE]

where we have introduced the collective index $\sigma\equiv[\mu_{1}\cdots\mu_{n};\ \ell_{f}m_{f};\ \ell_{i}m_{i}]$. We also take the convention that if $\sigma$ 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^{\sigma}$ is a universal quantity with no scheme dependence.

Evaluating the integral part of $G^{\sigma}$ 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\epsilon$ 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}\neq 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^{\sigma}$ in the case of $P_{i}\neq P_{f}$ case is detailed in Sec. III.2 below.

III.1 $P_{i}=P_{f}$

The function $G^{\sigma}$ simplifies considerably when the initial and final momenta coincide, $P_{i}=P_{f}\equiv P$, i.e. when the momentum transfer vanishes, $P_{f}^{\mu}-P_{i}^{\mu}=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,\textbf{k})$ becomes a double pole of the form

[TABLE]

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

[TABLE]

where $k^{\star}$ and $q^{\star}$ are defined with respect to the $P$ rest frame.

If $D(m,\textbf{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

[TABLE]

where

[TABLE]

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

[TABLE]

As we explain in Appendix B.1, $\mathcal{C}_{\sigma,JM}(\boldsymbol{\beta},k^{\star})$ can be efficiently calculated by writing the factors within $K^{\omega}$ as boosted c.m. frame vectors, $k^{\mu}={[\Lambda_{-\boldsymbol{\beta}}]^{\mu}}_{\nu}k^{\star\nu}$. Such factors can then be written as spherical harmonics and, using Clebsch-Gordon coefficients, these can be combined with the external factors of $\mathcal{Y}_{\ell_{f}m_{f}}(\textbf{k}^{\star})$ and $\mathcal{Y}^{*}_{\ell_{i}m_{i}}(\textbf{k}^{\star})$ to identify the a final harmonic basis. As a final note, we stress that it is possible to unambiguously separate the dependence on $\boldsymbol{\beta}=\textbf{P}/E$ and $k^{\star}$ within $\mathcal{C}_{\sigma}$, i.e. one can vary $k^{\star}$ while holding $\boldsymbol{\beta}$ constant. This will be important for the manipulations performed below.

The construction of $\mathcal{C}_{\sigma,JM}(\boldsymbol{\beta},k^{\star})$ is discussed in detail in Appendix B.1. As a specific simple example, here we consider the case $\sigma=[\mu;\,10;\,10]$. The numerator within $G_{\sigma}$ is then

[TABLE]

The current insertion, $k_{\mu}$, can be written as a combination of $\ell=0$ and $\ell=1$ spherical harmonics. Combining this with the two $\ell=1$ harmonics from the external states, one finds that $\mathcal{C}_{\sigma,JM}$ is zero for $J>3$.

The $JM=00$ component only has a non-zero contribution in the $k_{0}^{\star}$ component. Isolating this contribution and substituting the definition of the boost matrix, we reach

[TABLE]

The remaining nonzero coefficients, arising for $J\leq 3$, take a more complicated form in general, but simplify for $\textbf{P}\propto\hat{\textbf{z}}$. If we additionally focus on the $\mu=z$ component, then only three coefficients survive

[TABLE]

We revisit this case below to show how it enters our final construction for $G^{\mu}_{10;10}$.

Returning to the main line of argument, we now substitute Eq. (40) into the expression for $G_{\sigma}(P,P,L)$, Eq. (38), to reach

[TABLE]

The next step is to expand the $k^{\star}$ dependence in the first factor within the summand about the pole location

[TABLE]

In the second term we have rewritten the derivative with respect to $k^{\star 2}$ (to be evaluated afterwards at $k^{\star 2}=q^{\star 2}$) directly as a derivative with respect to $q^{\star 2}$. This is possible because $\mathcal{C}_{\sigma,JM}(\boldsymbol{\beta},k^{\star})$ has no implicit $q^{\star 2}$ dependence and it is formally possible to vary $q^{\star}$, and thus $E^{\star}$, while holding $\boldsymbol{\beta}$ constant.

Combining this with the definition of $\overline{D}(m,k^{\star},q^{\star})$ one can show

[TABLE]

Remarkably, the operator in parenthesis vanishes when acting on $1/\omega_{q2}^{\star}$ so that this reduces to

[TABLE]

It follows that Eq. (44) can be rewritten as

[TABLE]

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 $\mathbf{r}^{\star}=\mathbf{k}^{\star}L/(2\pi)$ and $x=q^{\star}L/(2\pi)$. Then, shuffling around terms and introducing a new geometric function, we conclude

[TABLE]

where

[TABLE]

These two equations give the main result of this sub-section. In Appendix B.4 we give some details about the evaluation of $\mathcal{Z}^{(n)}_{JM}(P,L,\alpha)$. In Fig. 3 we plot $\mathcal{Z}^{(n)}_{JM}(P,L,\alpha)$ for $\textbf{P}=(2\pi/L)[001]$, for various values of $J$ and two different volumes.

We give two final comments concerning the new kinematic function, $\mathcal{Z}^{(n)}$. First we note an advantage of the decomposition over a single spherical harmonic performed in Eq. (40). It is now straightforward to use the symmetries of the finite-volume system to identify for which values of $JM$, $\mathcal{Z}^{(n)}_{JM}$ will be nonzero. This is discussed in detail in Appendix B.5.

Second, we remark that the cutoff function used here is designed with the property that the $\mathcal{O}(\alpha)$ correction cancels the pole, and thus generates a smooth quantity within the sum-integral difference. If, for example, one were to instead use $e^{-\alpha(r^{\star 2}-x^{2})}$ for all $n$ values, this would still be formally valid, but would lead to corrections of the form $\alpha/L^{k}$, making it more difficult to estimate the $\alpha\to 0$ limit. In fact these can be systematically subtracted and, as we shown in Appendix B.3, this proves to be an efficient alternative approach for evaluating $\lim_{\alpha\to 0}\mathcal{Z}^{(n)}_{JM}$.

We close this sub-section by returning to the specific case discussed above, $\sigma=[\mu=z;\,10;\,10]$ and $\textbf{P}=(2\pi/L)[00d_{z}]$. Suppressing the arguments of $\mathcal{Z}^{(n)}$, the final result for $G^{\sigma}$ can be written as

[TABLE]

Note that, in the case of mass-degenerate particles, $\mathcal{Z}^{(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 B.5.] In Fig. 4 we plot the real and imaginary parts of $G^{\mu}_{10;10}(P,P,L)$ for $\textbf{P}=(2\pi/L)[001]$.

III.2 $P_{i}\neq P_{f}$

We now turn our attention to the more challenging general case of $P_{i}\neq P_{f}$. Note that this is realized if any of the four components of the 4-vectors differ, in particular also for $\textbf{P}_{i}=\textbf{P}_{f}$ but $E_{i}\neq 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^{\sigma}$, 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}\textbf{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 $\mathcal{Z}_{JM}^{(n)}$ in the previous subsection

[TABLE]

where $\bar{\alpha}=L^{4}/(2\pi)^{4}\alpha$ and $r^{\star 2}_{i}(2\pi/L)^{2}=k^{\star 2}_{i}$, $x^{2}_{i}(2\pi/L)^{2}=q^{\star 2}_{i}$, etc. We then write

[TABLE]

where

[TABLE]

The sum can be evaluated directly as written, and thus we make no further modifications to $\mathcal{S}_{\sigma}$. The remainder of this section is dedicated to $\mathcal{I}_{\sigma}$.

III.2.1 Separation into $\mathcal{I}_{\mathcal{A};\sigma}(P_{f},P_{i})$ and $\mathcal{I}_{\mathcal{N};\sigma}(\bar{\alpha},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 $\mathcal{I}_{\mathcal{A};\sigma}(P_{f},P_{i})$, and a smooth remainder that is well-suited to numerical evaluation, denoted $\mathcal{I}_{\mathcal{N};\sigma}(\bar{\alpha},P_{f},P_{i})$. To proceed we define

[TABLE]

and then introduce $D_{r}(m,\textbf{k})$ via

[TABLE]

In Appendix B.6 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_{\mu}$ and the spherical harmonics multiplying $D(m,\textbf{k})$. To include these, we first introduce a tensor, $M$, defined such that

[TABLE]

where $N=n+\ell_{f}+\ell_{i}$. A more explicit definition, together with various examples, is given in Appendix B.1. This simply amounts to recasting the factors of $\textbf{k}^{\star}$ within the harmonics as boosted factors of $k^{\mu}$.

To incorporate $K^{\omega}_{\nu_{1}\cdots\nu_{N}}$, note that $D_{r}$ receives two contributions, one each from the poles at $k^{0}=E_{f}+\omega_{P_{f}k1}-i\epsilon$ and $k^{0}=E_{i}+\omega_{P_{i}k1}-i\epsilon$. (The third pole, at $k^{0}=\omega_{k2}-i\epsilon$, generates the term we are after, $D(m,\textbf{k})$.) We thus define

[TABLE]

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\to i$ everywhere, and explicit expressions for $D_{ri}(m,\textbf{k})$ and $D_{rf}(m,\textbf{k})$ are given in App. B.6. With these quantities in hand, Eq. (57) generalizes to

[TABLE]

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

[TABLE]

The next step is to address the issue of ultraviolet convergence for these integrals. Equation (56) is manifestly convergent, due to the inclusion of the cutoff function $H(\bar{\alpha},\textbf{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. (62) and subtracting from this the same equation defined with $m_{1,2}\to\Lambda$

[TABLE]

The regularization scale, $\Lambda$, is chosen so that the integrands that depend on it are smooth functions of $k$ with no need for an $i\epsilon$ prescription. This holds for any $\Lambda$ satisfying $2\Lambda>\text{max}[E^{\star}_{f},E^{\star}_{i}]$ though in practice it is useful to take the cutoff well above this minimum value. On the left-hand side of Eq. (64) 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/|\textbf{k}|^{5}$ in the limit $|\textbf{k}|\to\infty$. These thus give convergent integrals with respect to $d^{3}\textbf{k}$.

This approach can be extended to any number of $k_{\mu}$ factors, simply by forming more complicated linear combinations to cancel all divergent powers

[TABLE]

where we have introduced $c_{0}=1$ and $\Lambda_{0}=\{m_{1},m_{2}\}$. As above, for $j>0$ we take $\Lambda_{j}$ such that the corresponding integrands are smooth functions of $k$ ($\Lambda_{j>0}>\text{max}[E^{\star}_{f},E^{\star}_{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}\textbf{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 B.8.

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

[TABLE]

where

[TABLE]

Equations (66)-(68) are the main results of this subsection.

We emphasize here that the integrands in the definition of $\mathcal{I}_{\mathcal{N}}$ 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(\bar{\alpha},\textbf{k})$ decays exponentially, together with the observation that $\mathcal{K}_{r}$ is a smooth function, as is $D(\Lambda_{j},\textbf{k})$ for $j>0$. For the first term, smoothness is guaranteed because $[H(\bar{\alpha},\textbf{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, $\mathcal{I}_{\mathcal{N};\sigma}(\bar{\alpha},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$, $\ell_{i}=\ell_{f}=0$, in which no factors of $k_{\mu}$ appear in the numerator of $\mathcal{I}_{\sigma}$. In this case the integrals are directly convergent, without any need for additional subtraction terms (i.e. the sum over $j$ in Eq. (65) reduces to the $j=0$ term). Equations (67) and (68) then reduce to

[TABLE]

This no-index version of $\mathcal{I}_{\mathcal{N}}$ is plotted in the left panel of Fig. 5.

We emphasize here that $\mathcal{I}_{\mathcal{N}}$, and thus also the original integral, diverges in the $\bar{\alpha}\to 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/|\textbf{k}|$ and the integrand as $d^{3}\textbf{k}/|\textbf{k}|^{3}$. In other words the convergence of $\mathcal{I}_{\mathcal{A}}$ is always better than that of the original integral by two powers of $k$, resulting from the off-shell integration of $k^{\mu}$. For fixing the subtraction scheme in Eqs. (67) and (68), it is only necessary that $\mathcal{I}_{\mathcal{A}}$ be rendered finite. [See also Appendix B.8.]

We close with one final example: $\sigma=[\mu=0;10;10]$, corresponding to a factor of $\mathcal{Y}_{10}k_{0}\mathcal{Y}^{*}_{10}$ in the numerator. This leads to an $\mathcal{I}_{\mathcal{A};\sigma}$ integral with an integrand scaling as $d^{4}kk^{2}/k^{6}$, i.e. diverging as $\log\Lambda$. Performing a single subtraction of the same integral with $m\to\Lambda$ is therefore sufficient to render the result finite. In fact, to improve the numerical evaluation of $\mathcal{I}_{\mathcal{N};\sigma}$, and to test our general method, here we choose to make two subtractions. As explained in Appendix B.8, one possible choice is to add an integral evaluated at $\Lambda=3m$ (with coefficient $-35/27$) and a second at $\Lambda=6m$ (with coefficient $8/27$). Implementing this in Eqs. (67) and (68) leads to convergent forms of $\mathcal{I}_{\mathcal{A};\sigma}$ and $\mathcal{I}_{\mathcal{N};\sigma}$ respectively, with integrands scaling as $dk/k^{5}$. $\mathcal{I}_{\mathcal{N};\sigma}$ in this scheme is plotted in the right panel of Fig. 5.

As we include additional factors of $k_{\mu}$ 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. (68) with an efficient numerical integrator. We give some details on our approach in Appendix B.6, but consider $\mathcal{I}_{\mathcal{N};\sigma}(\bar{\alpha},P_{f},P_{i})$ as a numerically known function for the remainder of the main text.

Thus it remains only to evaluate $\mathcal{I}_{\mathcal{A};\nu_{1}\cdots\nu_{N}}(P_{f},P_{i})$, to which we now turn.

III.2.2 Evaluating $\mathcal{I}_{\mathcal{A};\sigma}(P_{f},P_{i})$

We first use the tensor $M^{\nu_{1}\cdots\nu_{N}}_{\sigma}$, introduced in Eq. (59) above, to define a version of $\mathcal{I}_{\mathcal{A}}$ with no spherical-harmonic indices

[TABLE]

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

[TABLE]

where

[TABLE]

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

To evaluate the integral in Eq. (73), 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

[TABLE]

to reach

[TABLE]

where $\xi_{\nu_{1}\cdots\nu_{N}}\equiv(i)^{\delta_{{\nu_{1}},0}+\cdots+\delta_{{\nu_{N}},0}}$. Note here that the indices are not contracted between $\xi_{\nu_{1}\cdots\nu_{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

[TABLE]

As we show in Appendix B.7, this reduces to

[TABLE]

where

[TABLE]

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. (76) is given by

[TABLE]

Taken together, Eqs. (78)-(80) give the main result of this subsection.

As above, it is instructive to consider $\mathcal{I}_{\mathcal{A}}$ for $N=0$, i.e. with no factors of $k_{\mu}$ in the numerator. In this case Eq. (80) is evaluated with no $\chi$-derivatives and gives

[TABLE]

As noted at the end of the previous subsection, this integral is convergent in the $\delta\to 0$ limit, and simplifies to

[TABLE]

where we have substituted

[TABLE]

with

[TABLE]

The final analytic step is to evaluate the integral with respect to $y$. We do so via the identity

[TABLE]

where in the second equality we have introduced a short-hand for the result of the integral. Note that, as long as $\text{Im}f(x)$ is nonzero, we can safely set $\epsilon=0$ in these expressions. In the case that the imaginary part does vanish, then we use the relation $\lim_{\epsilon\to 0}\arctan(a/\epsilon)=\text{sign}[a]\pi/2$.

Applying this to Eq. (82) we deduce

[TABLE]

where $F^{(1)}(x)$ [also given in Eq. (212)] is

[TABLE]

In this case $\text{Im}y_{+}(x)=-\text{Im}y_{-}(x)=\text{Im}\sqrt{A^{2}+B}/2$. Thus the substitution $\lim_{\epsilon\to 0}\arctan(a/\epsilon)=\text{sign}[a]\pi/2$ is only needed when the argument of the square root is positive. As we explain in detail in Appendix B.9, evaluating the remaining integral over $x$ reveals that $\mathcal{I}_{\mathcal{A}}(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. 2(b) can go on shell.

More precisely, we show in Appendix B.9 that the singularity locations are governed by the discriminant of the polynomial $A^{2}+B\equiv ax^{2}+bx+c$, given by

[TABLE]

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}\equiv 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 Landau (1959).

At values of $P_{f}$ and $P_{i}$ satisfying these conditions, the real part of $\mathcal{I}_{\mathcal{A}}$ has a step-function discontinuity of height

[TABLE]

and the imaginary part shows a logarithmic divergence. In Fig. 6 we illustrate how these singularities form in the $\epsilon\to 0$ limit of the $i\epsilon$ pole prescription. In Fig. 7 we show the singularity structure as a function of $E_{f}^{\star}$ for various fixed values of $E_{i}^{\star}$. In particular one sees that, for sub-thresold $E_{i}^{\star}$, $\mathcal{I}_{\mathcal{A}}$ is a smooth function away from the two-particle production threshold. As $E_{i}^{\star}$ approaches $2m$ a step forms in $\text{Re}\mathcal{I}_{\mathcal{A}}$ and $\text{Im}\mathcal{I}_{\mathcal{A}}$ develops a log divergence. Then, as $E_{i}^{\star}$ is further increased, the location of the singularity in $E_{f}^{\star}$ 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. $\textbf{P}_{i}=\textbf{P}_{f}$, $\mathcal{I}_{\mathcal{A}}$ 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. 8. 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

[TABLE]

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 2(b), going on-shell. (See also Refs. Bayar et al. (2016); Coleman and Norton (1965).)

This concludes our discussion of $\mathcal{I}_{\mathcal{A}}$ within the main text. In Appendix B.10 we extend the results here by explicitly evaluating $\mathcal{I}_{\mathcal{A};\nu_{1}}$, $\mathcal{I}_{\mathcal{A};\nu_{1}\nu_{2}}$ and $\mathcal{I}_{\mathcal{A};\nu_{1}\nu_{3}\nu_{3}}$. For these integrals we find that the above-threshold discontinuities persist, but are milder when factors of $k_{\mu}$ appear in the numerator.

III.2.3 Examples of $G^{\sigma}\!(P_{f},P_{i},L)$

Having discussed the integral entering $G^{\sigma}\!(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. 8, we consider the case of $\textbf{P}_{i}=\textbf{P}_{f}=(2\pi/L)[001]$ with $\sigma=[\mu=0;10;10]$. Then, in Fig. 9 we take $\textbf{P}_{i}=[000]$ and $\textbf{P}_{f}=(2\pi/L)[001]$ with $\sigma=[00;00]$, i.e. with no numerator factors. In both cases we fix $mL=6$ and plot $\text{Re}\,G^{\sigma}$ for all values of $E_{i}^{\star}$ and $E_{f}^{\star}$ in the region of interest.

As explained in the figure captions, each example illustrates important issues and features that arise. In Fig. 8 we consider various diagonal slices of the $E_{i}^{\star}$, $E_{f}^{\star}$ plane. We find that the result for $E_{i}^{\star}=E_{f}^{\star}$ is in perfect agreement with the $P_{i}=P_{f}$ result determined by combining the various $\mathcal{Z}_{JM}$ functions. This provides a strong check on the two different methods. Figure 8 also illustrates that double poles arise along the $E_{i}^{\star}=E_{f}^{\star}$ line, but these split to single poles as the slice is rotated away from this singular choice.

In addition, Fig. 8 illustrates the results of using the Pauli-Villars-like regulator to separate the integral into $\mathcal{I}_{\mathcal{A}}$ and $\mathcal{I}_{\mathcal{N}}$. As discussed towards the end of Sec. III.2.1, the original integral contains a factor scaling as $\mathcal{Y}_{10}k^{\mu}\mathcal{Y}_{10}\sim|\textbf{k}|^{3}$ in the numerator leading to a $\log\Lambda$ divergence in $\mathcal{I}_{\mathcal{A}}$. Following Appendix B.8 we handle this by evaluating Eqs. (67) and (68) with $\{\Lambda_{1},\Lambda_{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 $\mathcal{I}_{\mathcal{N}}$. We already gave the result for $\mathcal{I}_{\mathcal{N}}$ in this prescription in the right panel of Fig. 5. Here the result is combined with the sum and $\mathcal{I}_{\mathcal{A}}$ to reach $G^{\sigma}$. With all building blocks summed together, the dependence on the Pauli-Villars parameters, $\Lambda_{i}$, cancels (between $\mathcal{I}_{\mathcal{A}}$ and $\mathcal{I}_{\mathcal{N}}$) as does the dependence on the smooth cutoff parameter, $\bar{\alpha}$, (between $\mathcal{I}_{\mathcal{N}}$ and the sum).

Turning now to Fig. 9, this result displays two additional features of $G^{\sigma}$. First we see that, when $\textbf{P}_{i}$ and $\textbf{P}_{f}$ differ, the non-interacting two-particle poles appear in different locations for $E_{i}^{\star}$ and $E_{f}^{\star}$. Interactions shift the finite-volume energies away from these singularities so that $G^{\sigma}$, 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 $\mathcal{I}_{\mathcal{A}}$. Such features are simply part of the correct definition of $G^{\sigma}$. Indeed, because the singularity structure is directly induced by the infinite-volume diagram of Fig. 2(b), it also appears within the infinite-volume $\textbf{2}+\mathcal{J}\to\textbf{2}$ transition amplitude itself. The steps, cusps and log divergences of $G^{\sigma}$ are present in both $\mathcal{W}$ and $\mathcal{W}_{\text{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.

IV Conclusion

In this work we have presented a modified version of the finite-volume formalism for studying $\textbf{2}+\mathcal{J}\to\textbf{2}$ transition amplitudes. This is closely related to the approach of Ref. Briceño and Hansen (2016), but differs in that all infinite-volume quantities are Lorentz covariant and the $\textbf{1}+\mathcal{J}\to\textbf{1}$ matrix elements have been reformulated in terms of standard form factors. As explained in Sec. II.2 and Appendix A, the new result is reached by making minor adjustments to the derivation presented in Ref. Briceño and Hansen (2016). For example, in that work finite-volume effects are expressed as sums over poles of the form $1/[2\omega_{k}(k^{0}-\omega_{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. III.2.2, 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. III.2.2 and Appendix B.9, and illustrated in Figs. 6, 7 and 9. We also recall that the form of $G$ presented in Ref. Briceño and Hansen (2016) carries four sets of spherical harmonic indices, resulting from a cumbersome description of the $\textbf{1}+\mathcal{J}\to\textbf{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. Briceño and Hansen (2016). 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 Chang et al. (2018); Savage et al. (2017). From our previous experience with spinning particles Briceño (2014); Briceño and Hansen (2016), 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. Ji (2013), and would open the possibility for lattice QCD calculations of two-body contributions to double-beta decays Tiburzi et al. (2017); Nicholson et al. (2018). 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. Christ et al. (2015), or in a spatial direction, as considered for example in Ref. Briceño et al. (2018). 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 Briceño and Davoudi (2013b); Hansen and Sharpe (2014); Mai and Döring (2017); Hammer et al. (2017); Briceño et al. (2017, 2018c).

Returning to the present formalism, several open questions persist that we plan to address in future work. For example, ${\mathcal{W}}_{\rm df}$ 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 ${\mathcal{W}}_{\rm df}$. For example, in the $\rho$ 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 Wilson et al. (2015b); Dudek et al. (2016); Moir et al. (2016); Woss et al. (2018); Dudek et al. (2014); Briceño et al. (2018b), it will be necessary to perform global fits of the matrix elements using Eq. (25). This requires a detailed understanding of the analytic structure of these amplitudes in which triangle singularities play a crucial role.

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

Appendix A Derivation of the covariant formalism

In this appendix we describe how the derivation of Ref. Briceño and Hansen (2016) must be modified to give the covariant version of the formalism, presented in Sec. II.2 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\to 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. 10(a). For $E_{i}^{\star},E_{f}^{\star}$ 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, $\textbf{k}=2\pi\textbf{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\to 2}(P_{f},P_{i})$ are encoded in the skeleton expansion of Fig. 10(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^{-mL}$.

The set-up here is identical to that of Ref. Briceño and Hansen (2016). 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. 10(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. Briceño and Hansen (2016)

[TABLE]

Here $\mathcal{L}(P_{f},k)$ and $\mathcal{R}^{\dagger}(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. Briceño and Hansen (2016), we define $\mathcal{G}_{L,\mu_{1}\cdots\mu_{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 $\textbf{1}+\mathcal{J}\to\textbf{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}\neq m_{1}^{2}$. [See also Eq. (21) 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. Briceño and Hansen (2016).

We now express $\mathcal{L}(P_{f},k)$, $\mathcal{R}^{\dagger}(P_{f},k)$ and $w^{\text{off}}_{\mu_{1}\cdots\mu_{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. Briceño and Hansen (2016), by first defining

[TABLE]

and then recombining the components with the $\mathcal{Y}_{\ell m}$ harmonics defined in Eq. (13)

[TABLE]

With these in hand we introduce the $\delta$ operator as follows

[TABLE]

These results match Eqs. (41) and (42) of Ref. Briceño and Hansen (2016), up to the minor notational differences discussed above. A key property that we will use below is that $[\mathcal{L}\delta](P_{f},k)$ vanishes like $(P_{f}-k)^{2}-m^{2}_{1}$ in the on-shell limit, and $[\delta\mathcal{R}^{\dagger}](P_{i},k)$ vanishes like $(P_{i}-k)^{2}-m^{2}_{1}$.

We now imitate this separation with the $\textbf{1}+\mathcal{J}\to\textbf{1}$ insertion, and in doing so introduce the first major difference as compared to our earlier work. Beginning with Eq. (21) of the main text, above, we introduce the shorthand $k_{f}\equiv P_{f}-k$ and $k_{i}\equiv P_{i}-k$ to write

[TABLE]

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

[TABLE]

where $\omega^{\star f}_{k2}$ is the temporal component of $k^{\mu\star f}$, given by boosting $k^{\mu}=(\omega_{k2},\textbf{k})$ to the $\textbf{P}_{f}=\textbf{0}$ frame, and $\omega^{\star i}_{k2}$ is the $f\to i$ analog.

The key point is that, when expressed in this way, $f^{(j)}$ has no dependence on $\hat{\textbf{k}}_{f}^{\star}$ or $\hat{\textbf{k}}_{i}^{\star}$. Thus there is no need to decompose in harmonics, nor to include barrier factors. The on-shell projection is simply

[TABLE]

To avoid clutter we have suppressed the arguments on the right-hand side, identical to those of Eq. (99). We stress here that the on-shell projections are subtle in that k is used to define two separate variables $k^{\star}_{f}$ and $k^{\star}_{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}^{\star}$ and $k_{i}^{\star}$ dependence, it is possible to separately project the initial and final states on shell via $k^{\star}_{i}\to q^{\star}_{i}$ and $k^{\star}_{f}\to q^{\star}_{f}$ respectively.

We next form linear combinations of the on- and partially-off-shell form factors in direct analogy to Eq. (37) of Ref. Briceño and Hansen (2016) to reach

[TABLE]

The individual terms on the right-hand side are defined by analogy to Eqs. (38)-(40) of Ref. Briceño and Hansen (2016) and have the explicit form

[TABLE]

Note that functions with a $\delta$ 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. Briceño and Hansen (2016), the key point here is that $\delta$ 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 $\delta$ 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

[TABLE]

where

[TABLE]

is the non-interacting covariant propagator. Our conventions are such that $\Delta(k^{2})$ has unit residue at the single particle pole, implying $\mathcal{S}(k^{2})$ is smooth and finite near $k^{2}=m_{1}^{2}$. This form of separation, in which $\mathcal{D}(k^{2})$ remains Lorentz invariant, is the second key distinction relative to our earlier formalism.

At this stage we have separated the endcaps [Eqs. (96), (97)] the single-particle form factor [Eq. (103)] and the covariant propagator [Eq. (107)] into on-shell terms plus corrections. Substituting these four identities into Eq. (93) then gives the analog of Eq. (43) of Ref. Briceño and Hansen (2016). 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\epsilon)$ singularity, 8 with $1/(k_{i}^{2}+m^{2}-i\epsilon)$, and a single maximally singular term that leads to the appearance of $G^{\sigma}\!(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 $\delta_{\text{df}}$ via

[TABLE]

and similarly for $\delta_{\text{df}}$ acting on the left.

We can make use of this to rewrite Eq. (93) as

[TABLE]

where $D(m,\textbf{k})$ is defined in Eq. (28). 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. (12). We have included the inverted symmetry factor, $\xi^{-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 $\mathcal{W}_{\text{df}}$ is defined with four subtraction terms, given by coupling the current to each of the four external propagators. When $\mathcal{W}_{\text{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 $\Delta(k_{f}^{2})$ as well as the alternative mass assignments within $D(m,\textbf{k})$, i.e. $\{m_{1},m_{2}\}\to\{m_{2},m_{1}\}$.

Using Eq. (23), we can write the first term on the right-hand side in terms of $G_{\mu_{1}\cdots\mu_{n};\ell_{f}m_{f};\ell_{i}m_{i}}(P_{f},P_{i},L)$, defined in Eq. (27). Following the steps outlined in Ref. Briceño and Hansen (2016), one finally arrives at Eq. (25) – the relationship between the finite-volume matrix element and $i{\mathcal{W}}_{{\rm df};\mu_{1}\cdots\mu_{n}}$, defined in Eq. (16). The latter emerges through the identity $\delta_{\text{df}}\mathcal{W}\delta_{\text{df}}\equiv\mathcal{W}_{\text{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\neq 2$, each admit a $\textbf{1}+\mathcal{J}\rightarrow\textbf{1}$ transition, then the definition of ${\mathcal{W}}_{{\rm df}}$, Eq. (16), is replaced by an expression with four subtractions

[TABLE]

Here the $w_{2;\mu_{1}\cdots\mu_{n}}$ in the bottom line have arguments $(k^{\prime},P_{i}-P_{f}+k^{\prime})$ and $(k^{\prime},P_{i}-P_{f}+k^{\prime})$ 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 ${\mathcal{W}}_{L,{\rm df}}$ and ${\mathcal{W}}_{\rm df}$. 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 $w_{2;\mu_{1}\cdots\mu_{n}}$ as a function of $(P_{f}-k,P_{i}-k)$ and decompose exactly as in the main text

[TABLE]

This then allows us to write

[TABLE]

where $G_{12}$ is exactly equal to the quantity defined in Eqs. (27) and (28) of the main text and $G_{21}$ is the same but with $m_{1}\leftrightarrow m_{2}$ everywhere.

Turning to the case of identical particles here the relation between between ${\mathcal{W}}_{{\rm df}}$ and ${\mathcal{W}}$ again has four subtractions, exactly as in Eq. (111). The four terms continue to be distinct due to the four different momentum assignments. However, the relation between ${\mathcal{W}}_{L,{\rm df}}$ and ${\mathcal{W}}_{\rm df}$ is exactly as Eq. (26) 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. 10(b), and that this has no symmetry factor, even in the case that the three hadrons in the loop are identical.

Appendix B 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_{\ell m;\ell^{\prime}m^{\prime}}(P,L)$ and $G^{\sigma}\!(P_{f},P_{i},L)$.

B.1 Index gymnastics

We first discuss various identities for rearranging spherical-harmonic and Lorentz indices in the evaluation of $G_{\sigma}(P_{f},P_{i},L)$. Begin with the case of $P_{i}=P_{f}$, in particular with Eq. (40) of the main text. Multiplying both sides by $(q^{\star})^{\ell_{i}+\ell_{f}}$ and substituting the definition of $\mathcal{Y}_{\ell m}(\textbf{k}^{\star})$, this reduces to

[TABLE]

The first aim of this subsection is to use this result to derive a useful expression for evaluating $\mathcal{C}_{\sigma,JM}(\boldsymbol{\beta},k^{\star})$.

Substituting the relation $k_{\mu}={[\Lambda_{-\boldsymbol{\beta}}]_{\mu}}^{\nu}k^{\star}_{\nu}$, we reach

[TABLE]

At this stage, the factor in square brackets is a simple polynomial in the coordinates $k_{x}^{\star}$, $k_{y}^{\star}$ and $k_{z}^{\star}$ with additional dependence on the magnitude entering through $\omega_{k2}^{\star}$. Thus, for a given set of indices, one can readily determine an explicit expression, and then use the orthogonality of the $Y_{\ell m}$s, to deduce

[TABLE]

As a specific example, we return to the case of $\sigma=[\mu=z;10;10]$ and $\textbf{P}=[00d_{z}]$, already discussed in the main text. For these indices, Eq. (117) reduces to

[TABLE]

Here we have used the fact that, with $\mu$ fixed to $z$, we only need to sum over one row of the boost matrix, ${[\Lambda_{-\boldsymbol{\beta}}]^{z}}_{\nu}=(\beta^{z}\gamma,0,0,\gamma)$. The integrals are now trivial to evaluate. For example, the $JM=00$ component reduces to

[TABLE]

where we have used $\beta^{z}=P^{z}/E$ and $\gamma=E/E^{\star}$. This matches Eq. (42) of the main text.

An alternative approach for determining $\mathcal{C}_{\sigma,JM}(\boldsymbol{\beta},k^{\star})$ is to note that the product of two spherical harmonics can be written in terms of Clebsch-Gordon coefficients

[TABLE]

and that each factor of $k^{\star}_{\nu}$ can be rewritten using the identity

[TABLE]

where the last equality defines ${T}^{\nu}_{\ell m}(k^{\star})$. Substituting this into the left-hand side of Eq. (116) 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 $\mathcal{C}_{\sigma,JM}(\boldsymbol{\beta},k^{\star})$.

We next consider the case of $P_{i}\neq 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_{\sigma}^{\nu_{1}\cdots\nu_{N}}$, in Eq. (71) above. The definition can be re-expressed as

[TABLE]

The second, and final aim of this subsection, is to derive a useful expression for $M_{\sigma}^{\nu_{1}\cdots\nu_{N}}$.

First we introduce a new set of tensors, denoted $\mathcal{T}$, that allow us to express the right-hand side in terms of 4-vectors

[TABLE]

The exact definition of the $\mathcal{T}^{\alpha_{1}\cdots\alpha_{\ell_{f}}}_{[\ell m]}$ can be inferred by comparing Eqs. (122) and (123). Note that they contain the $\sqrt{4\pi}/(q^{\star})^{\ell}$ prefactors and also encode the combinations of $\textbf{k}^{\star}$ components needed to form the various spherical harmonics. For example $\sqrt{4\pi}k^{\star}Y^{*}_{10}(\hat{\textbf{k}}^{\star})/q^{\star}=-\sqrt{3}k_{\mu=3}^{\star}/q^{\star}$ implies $\mathcal{T}^{\mu}_{[10]}=\sqrt{3}(0,0,0,-1)/q^{\star}$, since $k_{\mu}^{\star}=(\omega^{\star}_{k2},-\textbf{k}^{\star})$. The single-index tensor, $\mathcal{T}^{\alpha}_{[\ell m]}$, is closely related to $T^{\nu}_{\ell m}$, introduced above.

The final step is to boost all 4-vectors to the finite-volume frame. We deduce

[TABLE]

B.2 Evaluating $F_{\ell m;\ell^{\prime}m^{\prime}}(P,L)$

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

[TABLE]

Our aim here is to rewrite this in terms of $\mathcal{Z}^{(1)}$. With this in mind, we first observe

[TABLE]

where we have used the fact that $2\omega_{Pk1}(E-\omega_{Pk1}-\omega_{k2}+i\epsilon)$ can be replaced with $2E^{\star}(q^{\star 2}-k^{\star 2}+i\epsilon)/2\omega_{k2}^{\star}$ up to smooth terms in the integrand that lead to exponentially suppressed volume dependence. The definition of $B_{JM}^{[\ell m;\ell^{\prime}m^{\prime}]}$ can be inferred by comparing Eqs. (125) and (126) and is given more explicitly by

[TABLE]

where we have used Eq. (120) to reach the last equality.

The final step here is to note that the additional barrier factor, $(k^{\star}/q^{\star})^{\ell+\ell^{\prime}-J}$, appearing in Eq. (126), can in fact be set to 1. This is justified because $(k^{\star})^{J}Y_{JM}(\hat{\textbf{k}}^{\star})$ is an analytic function and because $\ell+\ell^{\prime}-J\geq 0$ for all nonzero $B_{JM}^{[\ell m;\ell^{\prime}m^{\prime}]}$. The latter point directly follows from the explicit expression of the $B_{JM}^{[\ell m;\ell^{\prime}m^{\prime}]}$ in terms of the Clebsch-Gordon coefficients. As a result, the difference $(k^{\star}/q^{\star})^{\ell+\ell^{\prime}-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

[TABLE]

B.3 $\alpha$-dependence of $\mathcal{Z}^{(n)}_{JM}$

Here we explain our choice of cut-off function used in the definition of $\mathcal{Z}^{(n)}_{JM}$, Eq. (51), recalled here for convenience

[TABLE]

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

By including the power of $n$ in the cutoff function, $e^{-\alpha(r^{\star 2}-x^{2})^{n}}$, we ensure that differentiating with respect to $\alpha$ gives a smooth summand. This, in turn, implies that $\partial_{\alpha}\mathcal{Z}^{(n)}$ vanishes up to terms that are exponentially suppressed in the volume

[TABLE]

Using this result in an expansion about $\alpha=0$ then gives

[TABLE]

Of course, it is possible to define the $\mathcal{Z}^{(n)}_{JM}$ functions with a milder cutoff function, for example

[TABLE]

For $n=1$ this equivalent to ${\mathcal{Z}}^{(n)}_{JM}(P,L,\alpha)$, but for $n>1$ it is a less useful prescription, due to an enhancement of the $\alpha$ corrections

[TABLE]

As ${\mathcal{Z}}^{(n-1)}_{JM}(P,L,0)$ is itself a singular function for $n>1$, we deduce that the difference between the optimal version, ${\mathcal{Z}}^{(n)}_{JM}(P,L,\alpha)$, and the alternative, $\widetilde{\mathcal{Z}}^{(n)}_{JM}(P,L,\alpha)$ can take on arbitrarily large values for any finite $\alpha$.

In Fig. 11 we compare the $\alpha$-dependence of ${\mathcal{Z}}^{(n)}_{JM}(P,L,\alpha)$ and $\widetilde{\mathcal{Z}}^{(n)}_{JM}(P,L,\alpha)$ for $n=1,2$, and show that the large $\alpha$-dependence of the latter is well described by Eq. (134). In the $\alpha\to 0$ limit the two prescriptions agree, but to optimize the numerical evaluation we advocate the form of Eq. (129) and use only this definition throughout the remainder of the text.

B.4 Evaluating $\mathcal{Z}^{(n)}_{JM}(P,L,\alpha)$

In this section we discuss the evaluation of $\mathcal{Z}^{(n)}_{JM}(P,L,\alpha)$, in particular the integral part of this quantity. As already mentioned in Sec. III.1, the integral entering $\mathcal{Z}^{(n)}$ vanishes for all entries besides $JM=00$, implying

[TABLE]

where

[TABLE]

For the single-pole function, $\mathcal{Z}^{(1)}$, this can be evaluated analytically and takes the form

[TABLE]

where ${\rm Erfi}(x)$ is the imaginary error function function, defined by $d{\rm Erfi}(z)/dz=2e^{z^{2}}/\sqrt{\pi}$ and ${\rm 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 $\Xi^{(n)}(x,\alpha)$ 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

[TABLE]

where

[TABLE]

We close with a final remark concerning the $n=2$ case, of direct relevance for $\textbf{2}+\mathcal{J}\to\textbf{2}$ transition amplitudes. Here the relevant integral is

[TABLE]

Recalling $x\propto q^{\star}\propto\sqrt{s-s_{\rm th}}$, where $s=P^{2}$ is the c.m. energy and $s_{\rm th}=(m_{1}+m_{2})^{2}$, we deduce that for $P_{i}=P_{f}$, $G_{\sigma}$ generically has an inverse square-root singularity at two-particle production threshold. This implies that $\mathcal{W}_{\text{df}}$, as well as $\mathcal{W}$, must have the same singularity. This behavior is visible in the values of $\mathcal{Z}_{JM}^{(2)}$ plotted in Fig. 3. In particular, we observe that $\mathcal{Z}_{JM}^{(1)}$ has milder behavior near threshold.

B.5 Symmetry contraints on $\mathcal{Z}^{(n)}_{JM}$

To efficiently implement the formalism it is useful to identify, based on symmetry arguments, the values of $JM$ for which $\mathcal{Z}^{(n)}_{JM}(P,L,\alpha)=0$. In this subsection we review these constraints and discuss subtleties that arise for $n>1$. Our results are summarized in Table 1.

We begin with $\textbf{P}=[000]$. The properties of the zero-momentum zeta function are well-known Lüscher (1991); Rummukainen and Gottlieb (1995), but we still think it useful to review the arguments here, in order to better understand the generalization to $\textbf{P}\neq[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\neq 00$ we write

[TABLE]

Substituting ${Y}_{JM}(-\hat{\textbf{r}})=(-1)^{J}{Y}_{JM}(\hat{\textbf{r}})$ then gives $\mathcal{Z}^{(n)}_{JM}((E,\textbf{0}),L,\alpha)=(-1)^{J}\mathcal{Z}^{(n)}_{JM}((E,\textbf{0}),L,\alpha)$ implying that the zeta function vanishes for all odd $J$. We can further rewrite the summand with $(r_{x},r_{y},r_{z})\longrightarrow(-r_{x},r_{y},r_{z})$ and use $Y_{JM}(\theta,\pi-\phi)=Y_{J,-M}(\theta,\phi)$ to show that $\mathcal{Z}^{(n)}_{JM}((E,\textbf{0}),L,\alpha)=\mathcal{Z}^{(n)}_{J,-M}((E,\textbf{0}),L,\alpha)$. Similarly, a $\pi/2$ rotation about the $z$-axis, together with the identity $Y_{JM}(\theta,\phi-\pi/2)=e^{-iM\pi/2}Y_{JM}(\theta,\phi)$ 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 $3r_{z}^{2}-r^{2}$. The sum of this structure (times a function of $\textbf{r}^{2}$) over the integer set, $\textbf{r}\in\mathbb{Z}^{3}$, is clearly identical to the same with $3r_{x}^{2}-r^{2}$ or with $3r_{y}^{2}-r^{2}$. Thus $\mathcal{Z}^{(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})-3r^{2}=0$, implying $\mathcal{Z}^{(n)}_{20}=0$ as claimed. The conditions under which $\mathcal{Z}^{(n)}_{JM}=0$ for $\textbf{P}=[000]$ are summarized in the first line of Table 1.

We now turn to non-zero momenta of type $\textbf{P}=[00d_{z}]$. As this momentum type only breaks the symmetry in the $z$ direction, the invariance under $(r_{x},r_{y},r_{z})\longrightarrow(-r_{x},r_{y},r_{z})$, as well as the $\pi/2$ rotation around the $z$-axis, give the same constraints as for $\textbf{P}=[000]$.999Of course, it is true in all cases that the sum over $\textbf{r}\in\mathbb{Z}^{3}$ is invariant under any octahedral transformation on r. The relevant question is whether this leads to a useful constraint on $\mathcal{Z}_{JM}^{(n)}$. Similarly, for $\textbf{P}=[0d_{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 $\textbf{r}\to-\textbf{r}$ no longer holds. For example, for $m_{1}=m_{2}$ and $\textbf{P}=[00d_{z}]$, the summand defining $\mathcal{Z}^{(n)}_{JM}$ now depends on

[TABLE]

This vector does not transform in a useful way under a flip of r. Remarkably, in the degenerate-mass case, the single pole functions, $\mathcal{Z}^{(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 $\textbf{P}=(2\pi/L)\textbf{d}$ and so we present the results for the general case. The approach is based on the results of Ref. Kim et al. (2005) and we begin by recalling Eqs. (62) and (66) from that work

[TABLE]

where $\textbf{R}^{\star}=\hat{\gamma}^{-1}(\textbf{r}-\textbf{d}/2)$. Here $\textbf{r}^{\star}_{\parallel}$ and $\textbf{r}^{\star}_{\bot}$ are the vector components parallel and perpendicular to d. We have also introduced the operator $\hat{\gamma}^{-1}(\textbf{r}^{\star}_{\parallel},\textbf{r}^{\star}_{\bot})=(\gamma^{-1}\textbf{r}^{\star}_{\parallel},\textbf{r}^{\star}_{\bot})$. Substituting these results into our definition of $\mathcal{Z}_{JM}^{(n)}$ we find

[TABLE]

Now note that, for $n=1$, this function exhibits two special features, both unique to the single-pole case. First, the factors of $\omega_{k2}/\omega_{k2}^{\star}$ multiplying the pole exactly cancel; second, the subleading term becomes a smooth function of the summed coordinate, r. We thus reach

[TABLE]

This simplified form, incidentally the form first derived for the moving frame quantization condition Rummukainen and Gottlieb (1995), makes the accidental symmetry that we are after manifest. In particular, we can now use that the sum over $\textbf{r}\in\mathbb{Z}^{3}$ is invariant under $\textbf{r}\longrightarrow\textbf{d}-\textbf{r}$. Under this transformation $\textbf{R}^{\star}\to-\textbf{R}^{\star}$ leaving $R^{\star}=|\textbf{R}^{\star}|$ unchanged. Thus every factor in the summand is invariant except for ${Y}_{JM}(-\hat{\textbf{R}}^{\star})=(-1)^{J}{Y}_{JM}(\hat{\textbf{R}}^{\star})$. We deduce that, for odd $J$ and degenerate masses, $\mathcal{Z}^{(1)}_{JM}=\mathcal{O}(e^{-mL})$. However, the remaining symmetries do not survive due to the factor of $\hat{\gamma}^{-1}$ in the definition of $\textbf{R}^{\star}$.

Finally we stress that the vanishing of odd $J$ due to the accidental symmetry only holds for $n=1$. As is clear from Eq. (146), for all other $n$ values, the ratio of omegas does not cancel, leading to another factor that changes under the $\textbf{r}\to\textbf{d}-\textbf{r}$ transformation. Thus, while the odd-$J$ $\mathcal{Z}^{(1)}_{JM}$ have no poles, for $n=2$ the functions already exhibit the full double-pole behavior. [See also Fig. 3.] 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}\neq P_{f}$, beginning with the numerical integral denoted $\mathcal{I}_{\mathcal{N};\sigma}(\alpha,P_{f},P_{i})$.

B.6 Evaluating $\mathcal{I}_{\mathcal{N};\sigma}(\alpha,P_{f},P_{i})$

Here we give some further details on the evaluation of $\mathcal{I}_{\mathcal{N};\sigma}(\alpha,P_{f},P_{i})$, defined in Eq. (68) 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, $\mathcal{I}_{\mathcal{A};\sigma}$, it gives the full integral entering the definition of $G_{\sigma}$.

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_{\mu}$ (but with no spherical harmonics) in the integrand. These are then contracted with the tensor $M^{\nu_{1}\cdots\nu_{N}}_{\sigma}$ to reach the final expression. While highly advantageous for the analytic integral, the Lorentz-index-based expressions lead to a costly determination of $\mathcal{I}_{\mathcal{N};\sigma}$. For example, for $\ell_{f}=\ell_{i}=1$ and a single current index, the covariant form contains factors of $k_{\nu_{1}}k_{\nu_{2}}k_{\nu_{3}}$ that naively require the evaluation of $64$ terms.

To improve this situation, it is preferable to move $M^{\nu_{1}\cdots\nu_{N}}_{\sigma}$ back inside the integrals defining $\mathcal{I}_{\mathcal{N};\sigma}$, and thereby rewrite the integrands to only carry the final $\sigma$ index. This procedure is a bit subtle, because the time components of the $k^{\mu}$ are evaluated at various values. To proceed we first recall that $\mathcal{I}_{\mathcal{N},\sigma}$ generally involves terms evaluated at the physical masses, $m_{1},m_{2}$, together with regulating integrals evaluated at some higher scale $\Lambda_{j}$. For the physical-mass terms, time components are evaluated at k^{0}\in\big{\{}\omega_{k2},E_{f}+\omega_{P_{f}k1},E_{i}+\omega_{P_{i}k1}\big{\}}, and for the regulating integrals at k^{0}\in\big{\{}\sqrt{k^{2}+\Lambda_{j}^{2}},E_{f}+\sqrt{(\textbf{P}_{f}-\textbf{k})^{2}+\Lambda_{j}^{2}},E_{i}+\sqrt{(\textbf{P}_{i}-\textbf{k})^{2}+\Lambda_{j}^{2}}\big{\}}. Some of these four vectors are recombined into the harmonics, $\mathcal{Y}_{\ell m}(\textbf{k}^{\star})$, and, because the time and space components mix upon boosting $k^{\star}_{\mu}={{[\Lambda_{\boldsymbol{\beta}}]_{\mu}}^{\nu}}k_{\nu}$, we end up with strange spatial components in some of the harmonics.

To give concrete expressions it is convenient to define

[TABLE]

Here $\Omega_{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, $\Lambda$, refers to this mass dependence. In the following we will use $m$ or $\Lambda_{0}$ to indicate that the $\omega$s are evaluated at their physical masses and $\Lambda_{j>0}$ to indicate evaluation at an unphysical value $m_{1}=m_{2}=\Lambda_{j>0}$. In short, the first two entries in $\mathcal{N}_{\sigma}$ simply serve to indicate the value at which all $k^{0}$ are evaluated. We stress that the tensor $M^{\nu_{1}\cdots\nu_{N}}_{\sigma}$ does not depend on these parameters but only on the c.m. frame energies $E_{i}^{\star}$ and $E_{f}^{\star}$ as well as the boost velocities $\boldsymbol{\beta}_{i}$ and $\boldsymbol{\beta}_{f}$. Thus, the only modification to the spherical harmonics is that they now depend on ${\textbf{k}}^{\star,\Omega}$, defined via

[TABLE]

If we set $\Omega_{k}=\omega_{k2}$ and $\Lambda=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 $\mathcal{I}_{\mathcal{N};\sigma}$

[TABLE]

where

[TABLE]

This is identically equal to the quantity defined in Eq. (68) of the main text. The only difference is that we have absorbed $M^{\nu_{1}\cdots\nu_{N}}_{\sigma}$ inside the integrands, via the new function $\mathcal{N}_{\sigma}$.

To complete the specification we require explicit expressions for $D_{rf}(m,\textbf{k})$ and $D_{ri}(m,\textbf{k})$

[TABLE]

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

We close this subsection with one final simplification to $\mathcal{I}_{\mathcal{N}}$. 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

[TABLE]

Begin the proof by considering the generic integral

[TABLE]

assumed to be convergent. Here we have separated the integrand into a Lorentz-scalar function, $\mathcal{F}$, together with an arbitrary function of the c.m. frame momenta, $\mathcal{G}$. 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 $\Lambda_{\boldsymbol{\beta}}$, on each index

[TABLE]

To simplify the right-hand side we perform a change of integration variable $k^{\prime\mu}\equiv{[\Lambda_{\boldsymbol{\beta}}]^{\mu}}_{\nu}k^{\nu}$ and also define $P_{f}^{\prime\mu}\equiv{[\Lambda_{\boldsymbol{\beta}}]^{\mu}}_{\nu}P_{f}^{\nu}$ and $P_{i}^{\prime\mu}\equiv{[\Lambda_{\boldsymbol{\beta}}]^{\mu}}_{\nu}P_{i}^{\nu}$

[TABLE]

Here we have used the Lorentz-invariance of the various building blocks, including the fact that $\textbf{k}_{i}^{\star}$ and $\textbf{k}_{i}^{\star}$ must be unchanged if we replace $k$, $P_{f}$ and $P_{i}$ with their boosted counterparts. This result, which can be rewritten as

[TABLE]

is just a statement of Lorentz-covariance for $\mathcal{X}$.

To conclude our demonstration of Eq. (154) we note that $\mathcal{I}_{\mathcal{N}}\equiv\mathcal{I}-\mathcal{I}_{\mathcal{A}}$ and that the two terms on the right-hand side each satisfy the functional form of $\mathcal{X}$, shown in Eq. (155). In the case of $\mathcal{I}$, the original integral defining $G$, one takes

[TABLE]

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 $\mathcal{I}$ and $\mathcal{I}_{\mathcal{A}}$ the spherical harmonics, as well as the cutoff function $H(\alpha,\textbf{k})$ can be absorbed into the definition of $\mathcal{G}$. Again the key point is that these objects are frame-independent because they carry a frame label, $\textbf{k}^{\star}(k,P)=\textbf{k}^{\star}(k^{\prime},P^{\prime})$. We deduce that $\mathcal{I}_{\mathcal{N}}$ must satisfy Eq. (158). Multiplying both sides by $\Lambda_{-\boldsymbol{\beta}}$, we conclude Eq. (154).

To see the power of this identity we take $\boldsymbol{\beta}=\boldsymbol{\beta}_{i}$, implying

[TABLE]

and thus that only one external direction enters the integral. In this case the integration coordinate is simply transformed to $\textbf{k}_{i}^{\star}$. We are then left with

[TABLE]

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^{\star\mu}_{i}$ are factored out, the remaining integrand can only depend on $\textbf{k}_{i}^{\star}$ 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\theta\equiv\hat{\textbf{k}}_{i}^{\star}\cdot\hat{\textbf{P}}_{f}^{\prime}$ and picking an aribtrary orientation for the azimuthal angle $\phi$, we reach

[TABLE]

Here the $\phi$ integral can be evaluated analytically using the rotation properties of the spherical harmonics. Thus only the $k_{i}^{\star}$ and $\theta$ integrals need to be performed numerically.

B.7 Evaluating $\mathcal{I}^{\chi}(P_{f},P_{i},m,\delta)$

In this subsection we show how the the generating functional, $\mathcal{I}^{\chi}(P_{f},P_{i},m,\delta)$, defined in Eq. (77) of the main text, reduces to the results given in Eqs. (78)-(79). We begin by inserting Feynman-parameter integrals into the definition to reach

[TABLE]

where in the second line we have performed the shift $k_{E}\to k_{E}+(xP_{f,E}+yP_{i,E})$ and have introduced

[TABLE]

Next note that the denominator of the integrand of Eq. (164) is invariant under $k_{E}\to-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

[TABLE]

To further reduce the expression note that we can make the substitution $(i\chi_{E}\cdot k_{E})^{2n}\longrightarrow A_{n}(-\chi_{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,\mu_{1}}\cdots k_{E,\mu_{2n}}$ (multiplied by a function of $k_{E}^{2}$) must be proportional to $\delta_{\mu_{1}\mu_{2}}\cdots\delta_{\mu_{2n-1}\mu_{2n}}+\cdots$ where the second ellipsis indicates a sum over all possible pairings. Contracting with $\chi_{E,\mu_{1}}\cdots\chi_{E,\mu_{2n}}$ then gives the claimed form.

To determine the normalization, first consider the case of $n=1$, corresponding to $k_{E,\mu}k_{E,\nu}\longrightarrow A_{1}\delta_{\mu\nu}k_{E}^{2}$. Taking the trace on both sides then gives $A_{1}=1/(4-\delta)$, where we are careful to consistently perform the calculation in $4-\delta$ dimensions. Similarly, for $n=2$ one finds

[TABLE]

First summing over $\mu_{1}=\mu_{2}$ and then over $\mu_{3}=\mu_{4}$ gives $A_{2}^{-1}=(1/3)[(4-\delta)^{2}+2(4-\delta)]$.

The result for general $n$ can be derived by first writing

[TABLE]

where in the second equality we have rewritten the integral with a dummy parameter $\alpha$, to be set to zero after differentiation. Next we multiply both sides by $\exp[-k_{E}^{2}]$ and integrate with respect to $dk_{E}k_{E}^{4-\delta}$ to write

[TABLE]

Solving for $A_{n}$ and evaluating the remaining integral using Schwinger parameters, we deduce

[TABLE]

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 $\alpha$ derivative and set $\alpha=0$.

Substituting into Eq. (166) then gives

[TABLE]

where we have returned to the Minkowski signature for $\chi$, $P_{f}$ and $P_{i}$ and have also defined

[TABLE]

To conclude we simplify $\mathcal{J}^{n}(P_{f},P_{i},m,\delta)$ by evaluating the momentum integral

[TABLE]

In the second step we have integrated over $d\Omega_{3-\delta}$ and then changed variables via $k^{2}_{E}\equiv M^{2}/\zeta-M^{2}$. The measure is modified as $2k_{E}dk_{E}=-M^{2}d\zeta/\zeta^{2}$ with the integral now running from $\zeta=0$ to $\zeta=1$. The final factor in the second line is just $\zeta^{3}/M^{6}=[k_{E}^{2}+M^{2}]^{-3}$.

Evalauting the $\zeta$ integral via

[TABLE]

we conclude

[TABLE]

This directly gives Eq. (78) in the main text.

B.8 Determining $c_{j}$ and $\Lambda_{j}$

As discussed in Secs. II.3.1 and II.3.2, an immediate application of the proposed formalism is electromagnetic reactions coupling two-pion states: $(\pi^{+}\pi^{0})_{I_{i}}+j_{\mu}\to(\pi^{+}\pi^{0})_{I_{f}}$. For this case we require the function $G_{\sigma}$ for $\sigma=[\ell_{f}m_{f};\ell_{i}m_{i}]$ (no current index) as well as $\sigma={[\mu;\ell_{f}m_{f};\ell_{i}m_{i}]}$ (one index) for $\ell_{i},\ell_{f}\leq 1$. This requires evaluating $\mathcal{I}_{\mathcal{A}}$ (with no indices) through $\mathcal{I}_{\mathcal{A};\nu_{1}\nu_{2}\nu_{3}}$. This set depends on only two scalar integrals, $\mathcal{J}^{0}$ and $\mathcal{J}^{1}$. The integral defining $\mathcal{J}^{0}$ is convergent so that we only need the $c_{0}=1$, $\Lambda_{0}=m$ term, i.e. no subtraction is required. The integral defining $\mathcal{J}^{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. (64), corresponding to $c_{1}=-1$ with $\Lambda_{1}=\Lambda$ equal to any value exceeding $2m$.

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 $\mathcal{I}_{\mathcal{A},\nu_{1}\cdots\nu_{N}}$. In a nut shell the approach requires identifying the divergent part of $\mathcal{J}^{n}$ and, by substituting this into the relation to $\mathcal{I}_{\mathcal{A},\nu_{1}\cdots\nu_{N}}$, to identify an expression of the form

[TABLE]

The coefficients $c_{j}$ and the scales $\Lambda_{j}$ are then tuned such that this quantity vanishes. Note that $\left\lfloor{N/2}\right\rfloor$ terms must be tuned to vanish. This is because the integral $\mathcal{I}_{\mathcal{A};\nu_{1}\cdots\nu_{N}}$ depends on the scalar integrals up to $n=\left\lfloor{N/2}\right\rfloor$. The integral then scales as $d^{4}k_{E}k^{2n}_{E}/k_{E}^{6}$ and generates a series of divergences of the form $\Lambda^{2n-2},\Lambda^{2n-4},\cdots,\log(\Lambda)$, so a total of $n=\left\lfloor{N/2}\right\rfloor$ terms.

It turns out that one does not need to tune both the regularization scales $\Lambda_{j}$ and the coefficients $c_{j}$. We thus choose the recipe of setting $\Lambda_{j}=2^{j-1}\Lambda$ [$\Lambda_{1}=\Lambda\,,\ \Lambda_{2}=2\Lambda\,,\ \Lambda_{3}=4\Lambda\,,\cdots$] and tuning only the $n$ different $c_{j}$ terms.

The latter is achieved by studying the $1/\delta$ terms. In particular one can show that $\mathcal{J}^{n}\sim M(m,x,y)^{2n-2}/\delta$. 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(\Lambda,x,y)^{2}-(1+c_{1})M(2\Lambda,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 $\Lambda$ dependence

[TABLE]

and requiring this to vanish gives $c_{1}=[m^{2}-4\Lambda^{2}]/[3\Lambda^{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 $\mathcal{I}_{\mathcal{A};\nu_{1}\nu_{2}\nu_{3}}$, beyond this the scheme in the non-degenerate case becomes more complicated.

Note that the choices of coefficients that remove $\mathcal{J}^{2}$ divergences also automatically remove those in $\mathcal{J}^{1}$. This is because $\sum_{j}c_{j}M(\Lambda_{j},x,y)^{2}$ contains both $\Lambda_{j}$-dependent and independent pieces. Thus the vanishing of $\sum_{j}c_{j}M(\Lambda_{j},x,y)^{2}$ guarantees that the same holds for $\sum_{j}c_{j}$, and the latter is the condition for removing divergences in $\mathcal{J}^{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 $\mathcal{J}^{n}$ divergences. For the general-$n$, degenerate case, the system that we need to solve is

[TABLE]

or equivalently

[TABLE]

where $C(n,k)=\frac{n!}{(n-k)!k!}$ is the binomial coefficient. Defining $v_{k}\equiv C(n,k)M(0,x,y)^{2k}$ and $A_{kj}\equiv 4^{(j-1)(n-k)}\Lambda^{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

[TABLE]

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 2 we give the values up to $n=5$, assuming $\Lambda_{j}=2^{j-1}(3m)$.

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

[TABLE]

where

[TABLE]

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 $\langle\langle M(0,x,y)^{2k}\rangle\rangle$ 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 $\chi$. Again we stress that this is not of immediate concern as it is only relevant for form factors with many indices.

B.9 Triangle singularities within $\mathcal{I}_{\mathcal{A}}(P_{f},P_{i})$

In this appendix we give a more detailed discussion of the singularities that arise away from threshold in $\mathcal{I}_{\mathcal{A}}(P_{f},P_{i})$, as summarized around Eqs. (90) and (91) of the main text and in Figs. 6 and 7. The task is to study discontinuities of the integral

[TABLE]

These arise when the three particles in the triangle loop of Fig. 2(b) can all go on shell. As Landau described in Ref. Landau (1959), the on-shell condition is realized at critical points $(x_{c},y_{c})$ of $M(m,x,y)^{2}$ defined by three conditions: ${\partial M(m,x_{c},y_{c})^{2}}/{\partial x}=0$, ${\partial M(m,x_{c},y_{c})^{2}}/{\partial 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. (187) with respect to $y$, to reach

[TABLE]

where

[TABLE]

with $y_{\pm}(x)\equiv(1/2)(A\pm\sqrt{A^{2}+B+i\epsilon})$. Here $A$ and $B$ are known functions of the kinematic variables and the Feynman parameter $x$, defined in Eqs. (85) and (86).

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\propto(x-x_{c})^{2}$. Before considering this special case we take our general form and substitute $A^{2}+B=2[(P_{i}\cdot P_{f})^{2}/s_{i}^{2}-s_{f}/s_{i}](x-x_{1})(x-x_{2})$

[TABLE]

For concreteness, here we have assumed $(P_{i}\cdot 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 $\text{Im}\sqrt{A^{2}+B}=0$ ($x\in[0,x_{1}]\cup[x_{2},1]$) as well as the region where it is non-zero ($x\in[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 $\epsilon=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}\in[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 [math] 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 $\int_{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 $\mathcal{I}_{\mathcal{A}}$. In addition, the imaginary part diverges as $\log\lvert P_{f(i)}-P_{f(i),c}\rvert$.

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}\to 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

[TABLE]

Next we note that, as $x_{1}$ approaches $x_{2}$, $\text{Im}y_{+}$ goes to zero. Thus it is natural to expand in this quantity

[TABLE]

We see that only the first term will contribute in the limit $x_{1}\to 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

[TABLE]

B.10 Evaluating $\mathcal{I}_{\mathcal{A}}(P_{f},P_{i})$ through $\mathcal{I}_{\mathcal{A};\nu_{1}\nu_{2}\nu_{3}}(P_{f},P_{i})$

Here we provide compact expressions here for the $\mathcal{I}_{\mathcal{A};\sigma}(P_{f},P_{i})$ integrals with up to three Lorentz indices. Starting with Eqs. (78) and (80) one can show

[TABLE]

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

[TABLE]

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

To evaluate $\mathcal{I}_{\mathcal{A};\nu_{1}}(P_{f},P_{i},m)$ we require

[TABLE]

For $\mathcal{I}_{\mathcal{A};\nu_{1}\nu_{2}}(P_{f},P_{i},m)$ we need

[TABLE]

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 $\int\,dx\,dy\,y$.

Finally, for $\mathcal{I}_{\mathcal{A};\nu_{1}\nu_{2}\nu_{3}}(P_{f},P_{i},m)$, six integrals appear

[TABLE]

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}^{\nu}\int\,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

[TABLE]

where $\mathcal{L}_{\pm\epsilon}[f(x)]$ is defined in Eq. (87). $\log_{\epsilon}(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_{\pm\epsilon}(f(x))=\pm i\,\pi\,+\log(|f(x)|)$. As discussed in the main text, these expressions allow one to work identically at $\epsilon=0$ in the numerical evaluation of the integrals with respect to $x$. The only memory of the non-zero value of $\epsilon$ that these functions carry are the sign.