Consistency checks for two-body finite-volume matrix elements: I. Conserved currents and bound states

Recently, a framework has been developed to study form factors of two-hadron states probed by an external current. The method is based on relating finite-volume matrix elements, computed using numerical lattice QCD, to the corresponding infinite-volume observables. As the formalism is complicated, it is important to provide non-trivial checks on the final results and also to explore limiting cases in which more straightforward predications may be extracted. In this work we provide examples on both fronts. First, we show that, in the case of a conserved vector current, the formalism ensures that the finite-volume matrix element of the conserved charge is volume-independent and equal to the total charge of the two-particle state. Second, we study the implications for a two-particle bound state. We demonstrate that the infinite-volume limit reproduces the expected matrix element and derive the leading finite-volume corrections to this result for a scalar current. Finally, we provide numerical estimates for the expected size of volume effects in future lattice QCD calculations of the deuteron's scalar charge. We find that these effects completely dominate the infinite-volume result for realistic lattice volumes and that applying the present formalism, to analytically remove an infinite-series of leading volume corrections, is crucial to reliably extract the infinite-volume charge of the state.


I. INTRODUCTION
One of the overarching goals of modern-day nuclear physics is the characterization and fundamental understanding of the low-lying strongly interacting spectrum. There is, by now, overwhelming evidence that the detailed properties of all low-lying states are governed by the dynamics of quark and gluon fields in the mathematical framework of quantum chromodynamics (QCD). But still, it remains a significant challenge to extract low-energy predictions from the underlying theory.
The vast majority of QCD states emerge as either bound states or resonances of multihadron configurations. An example is the deuteron, a shallow bound state of the isoscalar proton-neutron channel with a binding energy of m n þ m p − M d ≈ 2.2 MeV. The deuteron has long been hypothesized to be a molecular state of the two nucleons [1], and similar pictures have been proposed for a variety of other QCD states. (See Ref. [2] for a recent review.) However, in many cases a straightforward interpretation is unavailable. For example, the isoscalar f 0 ð980Þ resonance couples strongly to ππ and KK states, and has been postulated to be both a tetraquark [3] and a KK molecule [4]. 1 The challenge of resolving the inner structure of composite hadrons is twofold: First, QCD is nonperturbative, so that systematic low-energy calculations are challenging. This has been addressed with substantial success using low-energy effective theories, methods based in amplitude analysis and numerical calculations using lattice QCD (LQCD). In contrast to the first two methods, LQCD has the unique advantage of relating the fundamental QCD Lagrangian to low-energy predictions. Second, composite states generally manifest as dynamical enhancements of multihadron scattering rates, meaning that the detailed observation depends on the production mechanism and decay channel of the resonance in question. This ambiguity is resolved, at least in principle, by recognizing that across all production and decay channels, a given resonance always leads to the same pole in an analytic continuation of scattering amplitudes to complex energies.
These two points have motivated the community to develop a systematic framework for extracting hadronic scattering amplitudes via LQCD. From the energy dependence of such amplitudes one can then quantitatively describe the bound and resonant states of the theory. In addition, by extracting transition amplitudes involving external currents, one can in principle access structural information of these states. In this work, we focus on an example in the latter class of the amplitudes, namely 2 þ J → 2 transition amplitudes. We consider a method, first introduced in Refs. [7,8], that allows one to determine such quantities from numerical LQCD.
Electroweak interactions involving scattering states can also be accessed using LQCD, via a generalization of the methods described above. The seminal example in this sector is the work of Ref. [65], providing a formal method for determining the electroweak decay, K → ππ. More generally, in processes for which the effects of the electroweak sector can be treated perturbatively, the relevant amplitudes are given via the evaluation of QCD matrix elements, built from the appropriate currents together with multiparticle external states. These ideas have been successfully developed for the case that either the initial or the final state couples strongly to two-particle scattering states [16,19,[65][66][67][68][69] and implemented in number lattice QCD studies, most prominently to determine the K → ππ decay amplitudes [70][71][72][73] as well as the electromagnetic process πγ ⋆ → ππ [74][75][76][77][78]. This progress motivates the consideration of more complicated electroweak transitions, in particular those with two hadrons in both the initial and the final states.
As we discuss in detail in Sec. II, 2 þ J → 2 transition amplitudes allow one to extract elastic form factors of bound states and resonances, thereby providing direct information on the structure of these states and possibly resolving which models are most descriptive [79][80][81]. As compared to the transitions described in the preceding paragraph, the necessary formalism for these quantities is significantly more complicated [7,8,20,82]. 3 Building on previous work, in Ref. [7] two of us derived a modelindependent relation between the corresponding finitevolume matrix elements, schematically denoted h2jJ j2i L (where L indicates the side length of the periodic cubic volume), and the 2 þ J → 2 transition amplitude, W. In Ref. [8] we improved the method by simplifying technical details relating to the on-shell projection of the singleparticle form factor and by using Lorentz covariant poles in the various finite-volume kinematic functions that arise. We stress that the two approaches are equivalent and only differ in the exact definitions of unphysical, intermediate quantities. The results are derived to all orders in the perturbative expansion of a generic relativistic field theory, for any type of two-scalar channels, with generalizations to spin and coupled channels left to future work. Details of this formalism are reviewed in Sec. III.
The purpose of this work is to provide two nontrivial checks on the general relations of Refs. [7,8], and also to demonstrate their predictive power even in simplified special cases. As a first check, in Sec. IV we demonstrate that the method is consistent with the consequences of the conserved vector current. In particular, the formalism predicts that the charge of a two-hadron finite-volume state is exactly equal to the sum of the constituent charges and independent of L. This relies on nontrivial relations between various L-dependent geometric functions, and a relation between the 2 → 2 and 2 þ J → 2 amplitudes that follows from the Ward-Takahashi identity. The second check, presented in Sec. V B, considers the analytic continuation of the formalism below two-particle threshold, for theories with an S-wave bound state. We show that the finite-and infinite-volume matrix elements coincide (once normalization factors are accounted for) up to term scaling as e −κ B L , where κ 2 B ¼ m 2 − M 2 B =4 defines the binding momentum for two constituents of mass m, binding to a mass of M B .
Presently, LQCD calculations of light nuclei properties are being performed at unphysically heavy quark masses, for which the binding momenta exceed their real-world values [41,47,[84][85][86]. In addition the properties of states can be shifted, e.g., the dineutron, in nature a virtual bound state, is found to be a standard bound state for m π ≳ 450 MeV [39,40,[87][88][89]. The increased binding suppresses finite-volume effects, and this has permitted exploratory calculations of matrix elements of these states [90][91][92][93][94], in which volume effects are ignored.
As LQCD calculations of multinucleon systems move toward physical quark masses, the binding momenta of the nuclei decrease, and it is well-known that finite-volume effects of the naively extracted states can become a dominant source of systematic uncertainty [18,95,96]. In the case of spectroscopy, an infinite series of e −κ B L corrections can be removed by applying the Lüscher formalism, as was done in [41,86] as well as in a wide variety of mesonic channels where bound states appear [42,46,53,54,57,58]. The results of this work stress that it is important to pursue the same paradigm for matrix elements of loosely bound states, using the formalism of Refs. [7,8] to nonperturbatively remove binding-momentum-enhanced finite-volume artifacts. To illustrate this point, in Sec. V B we determine the leading e −κ B L corrections and compare these to the full result, which holds up to e −mL . Finally, in Sec. V C we present a numerical example, meant to model the deuteron at physical pion masses, and show that the full formalism is needed to reliably remove the L dependence for box sizes in the region of mL ≈ 4-7. Otherwise the e −κ B L corrections can become comparable in size with the infinite-volume result and thereby dominate the systematic uncertainties.
Though largely addressed above, we close here with a brief summary of the remaining sections. After reviewing basic properties of the infinite-volume 2 → 2 and 2 þ J → 2 amplitudes in Sec. II, in Sec. III we described the corresponding finite-volume formalism for each type of amplitude. Then, in a very compact Sec. IV, we demonstrate that the finite-volume 2 þ J → 2 formalism gives the expected results for matrix elements of a conserved current. Section V is dedicated to volume effects on a twoparticle bound state, including a check that the L → ∞ limit gives the required result, a calculation of the leading Oðe −κ B L Þ corrections, and a numerical exploration intended to guide future LQCD calculations of the deuteron's scalar charge. We briefly conclude in Sec. VI. In addition, this article includes three Appendixes, providing proofs of various technical results used in the main text.

II. INFINITE-VOLUME AMPLITUDES AND BOUND STATES
In this section we review the definitions and key properties of the infinite-volume 2 → 2 and 2 þ J → 2 amplitudes, with particular attention to the expressions relevant for an S-wave bound state. For simplicity, we focus on systems composed of two scalar particles, with degenerate mass m, distinguished by their charge with respect to an external current J μ . One of the particles carries charge Q 0 , while the other is neutral. Here we have in mind a scalar analog of the proton-neutron system.

A. 2 → 2 amplitudes and bound-state poles
In a general Lorentz frame, the two-particle system has a total energy-momentum denoted by P ¼ ðE; PÞ. Boosting to the center-of-momentum frame (CMF) we define P ⋆ ¼ ðE ⋆ ; 0Þ, which is related to the Mandelstam variable s and a generic P by Two-particle scattering is described by s, as well as the back-to-back momentum orientations of the initial and final states in the CMF:k ⋆ i andk ⋆ f , respectively. 4 Using these coordinates we can introduce the scattering amplitude and its partial wave expansion We have used that total angular momentum, l, is conserved, and that the partial-wave amplitude is independent of the projection, m l , both consequences of rotational symmetry. In the following we are interested in the case of a scalar bound state, appearing as a subthreshold pole in M l¼0 ðsÞ. We therefore restrict attention to the S-wave (l ¼ 0) amplitude and do not write the angular momentum index on the partial wave amplitude for the rest of this section. The elastic 2 → 2 scattering amplitude can be represented in terms of the K matrix, which is an on-shell representation that enforces S matrix unitarity explicitly below the inelasticity threshold [97], Here, ρðsÞ is the two-body phase space, encoding the onshell propagation of two particles. It is defined as where q ⋆ is the relative momentum of the two particles in the CMF, q ⋆ ≡ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s=4 − m 2 p . This square root introduces a branch cut in the complex s plane, illustrated in Fig. 1.
Bound states are then defined as subthreshold poles on the first Riemann sheet, the sheet for which Imq ⋆ > 0. 5 The K matrix, KðsÞ, is a real function describing all of the dynamics of the system. It can be written in terms of the scattering phase shift, δðsÞ, via Unlike MðsÞ, the K matrix is an analytic function of s in a domain around s ¼ ð2mÞ 2 set by the nearest left-hand cut. It follows that the effective range expansion has a finite radius of convergence, and gives a useful description of KðsÞ and MðsÞ near threshold. The parameters a and r are called the scattering length and effective range, respectively. Using Eqs. (3)-(5), the condition for a bound state (a real subthreshold pole on the first Riemann sheet) can be expressed as where κ B is the binding momentum, related to the pole position s B via s B ¼ 4ðm 2 − κ 2 B Þ where the mass of the bound state is given as Going beyond the pole, information about the nature and structure of the bound state is also contained in its coupling to two-particle scattering states, g, defined as the residue of the pole As we review in Sec. VA, g governs the prefactor of the bound state's leading finite-volume effects.

B. 2 + J → 2 amplitudes and form factors
We now turn to the less standard 2 þ J → 2 transition amplitude, defined via Here the initial and final states have kinematics as in Sec. II A, and the current, J μ , is a local operator evaluated in position space at the origin. Since the current can inject energy and momentum, the initial and final states carry different total four-momenta, P i and P f , respectively. It is also convenient to define the squared momentum transfer, where the overall minus is included so that Q 2 > 0 for spacelike P f − P i . The amplitude W μ can be defined for local currents with any Lorentz structure, and in Sec. V B 2 we also consider specific results for a scalar current. Here, for concreteness we focus on a conserved vector current J μ ðxÞ satisfying Our first aim is to connect this amplitude to the bound-state form factor, defined via where jP i ; Bi is the bound state, normalized as (11) is related to the S-wave projection of W μ , analytically continued below threshold to the bound-state pole: where s i;f ¼ P 2 i;f . We prove this result in Appendix A. As discussed in some detail in Refs. [7,8], the analytic structure of W μ is significantly more complicated than that of M. One can identify three generic sources of nonanalyticity in the 2 þ J → 2 amplitude: (i) Exactly as for M, diagrams with on-shell two-particle intermediate states lead to factors of ρðsÞ. W μ thus exhibits the same branch cut and multisheet structure as M. (ii) Isolated poles arise due to the subtracted diagrams in Fig. 2(b), in which the current is attached to an external leg. (iii) The triangle diagram shown in Fig. 2(c) induces a new class of singularities, first described by Landau in Ref. [98].
To relate the 2 þ J → 2 amplitude to a physical scattering rate, it is important to recall that the former arises from perturbatively expanding the weakly interacting sector (encoded in J μ ) while keeping the strong dynamics nonperturbative. In fact, both the triangle singularities and the isolated poles can be understood as artifacts, resulting from truncating the expansion at a fixed order. Of course, the more standard single-particle form factors also arise from such an expansion, but happen to exhibit more straightforward analytic structure in the kinematic region considered. Indeed, as we will see below, analytic continuation to the bound-state pole removes all three of the nonanalyticities we have identified. As first explained in Ref. [7], finite-volume matrix elements are more directly related to a subtracted amplitude, denoted W df , from which the isolated poles, item (ii) above, have been removed. Here the subscript "df" stands for "divergence free." The definition, depicted also in Fig. 2 where w μ is the single-particle matrix element, and we have also introduced the corresponding form factor, f. Here we adopt the convention that the charged particles carry momenta P i − k (incoming) and P f − k 0 (outgoing). 6 The overline inM denotes a slight modification to the definition of M to account for the off-shell leg. This is described in Ref. [8] and, since the distinction is irrelevant for the S-wave amplitude, we do not discuss the issue further in this work. The momentum directions within W μ df can be projected to definite angular momentum as done in Eq.
As with the scattering amplitude, for the remainder of this section we restrict attention to the S-wave component of W μ df , i.e., the component containing our bound state. Note that, in contrast to M, W df has off-diagonal elements in angular momentum space, due to the angular momentum injected by the external current. A crucial observation that will guide our later analysis is that the S-wave component of W df exactly satisfies Eq. (17) above, i.e., W μ and W μ df have the same bound-state double pole with the same residue: This equivalence holds because the subtracted terms in Eq. (13) only have a single pole, and thus cannot modify the leading divergence. The bound-state poles within W μ df motivate us to introduce a new object, F μ ðP f ; P i Þ, given by The S-wave scattering amplitudes on each side remove the poles from W μ df ðP f ; P i Þ, implying lim In addition, Eq. (18) factorizes the ρðsÞ branch cuts from W μ df so that F μ does not contain this class of singularities [80,81]. Nonetheless, F μ is in general complex, due to the triangle diagram of Fig. 2(c). This diagram can only contribute complexity (as well as nonanalyticity) when at least one of the two-particle cuts goes on shell, i.e., when such an intermediate state can physically propagate. In particular, for subthreshold energies, and thus for some domain around the bound-state energy, the triangle integral is real and analytic.
At this stage we have argued that each of the three nonanalyticities listed above is irrelevant near the bound-state Representation of (left) the 2 þ J → 2 amplitude with initial and final energy-momentum given by P i and P f , respectively, leading to a current-induced momentum transfer of pole: First, the external-leg poles are removed in the conversion from W to W df ; second, the on-shell threshold cuts in ρðsÞ are removed in the relation between W μ df and F μ ; and finally, the triangle singularity (still contained in F μ ) is avoided by the subthreshold continuation.
We close this section with an important property of W μ df that follows from the Ward-Takahashi identity. As we prove in Appendix B, W μ df satisfies the following simple relation to the scattering amplitude: This identity is crucial for the analysis presented in Sec. V B 1.

III. FINITE-VOLUME FORMALISM FOR TWO-PARTICLE SYSTEMS
Before giving detailed expressions for the finite-volume effects on a two-body bound state, in this section we briefly review the general formalism describing the finitevolume energies and matrix elements of two-particle systems. In the following, we work in a cubic, periodic volume of length L with an infinite temporal extent. The total momentum of the system in the finite-volume frame is allowed to take on any value consistent with the periodicity:

A. Finite-volume energies
In the window of energies for which only two particles can propagate, the finite-volume spectrum is related to the infinite-volume partial-wave amplitudes, defined in Eq. (3), via the Lüscher quantization condition [14][15][16]. Generally, the quantization condition is a determinant over angular momentum space. If we neglect waves higher than l ¼ 0, however, it reduces to a simple algebraic relation where s n ¼ P 2 n ¼ E n ðLÞ 2 − P 2 corresponds to the eigenenergy of the nth finite-volume two-particle state. Here FðP; LÞ is a known finite-volume function, where and ω Pk ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m 2 þ ðP − kÞ 2 p are the on-shell energies of the two particles, and 1 Equation (21) holds up to corrections associated with higher partial waves and only for s n below the first inelastic threshold.

B. Finite-volume matrix elements
Similarly, one can relate finite-volume matrix elements of two-particle systems to infinite-volume 2 þ J → 2 transition amplitudes. Here, the relevant formalism was first derived in Ref. [7] using an all-orders perturbative expansion based in a generic relativistic effective field theory. Recently, in Ref. [8], the formal approach was improved in two ways: First, by rearranging the separation of finite-volume effects, we were able to show that the extracted infinite-volume transition amplitudes are manifestly Lorentz covariant. Second, we reorganized the analysis so that single-particle matrix elements enter via standard form factors (rather than a nonstandard spherical harmonic decomposition used in the first publication). While the two representations are formally equivalent, the work of Ref. [8] is expected to be significantly more convenient in numerical applications going forward. Of course, all expressions used here are taken from the improved approach.
Again, assuming all but the l ¼ 0 partial waves are negligible, the matrix elements of the vector current for two-particle states can be related to W μ df , defined in Eq. (18), as follows: L 3 hP n;f ; LjJ μ jP n;i ; Li ¼ W μ L;df ðP n;f ; P n;i ; LÞ × ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi RðP n;f ; LÞRðP n;i ; LÞ where P n;i ¼ ðE n;i ; P i Þ and P n;f ¼ ðE n;f ; P f Þ. Here W μ L;df is an L-dependent function related to W μ df in a manner detailed in the following paragraph. In addition, R is a generalization of the Lellouch-Lüscher factor [65], first introduced in Ref. [67] RðP n ; LÞ ¼ ∂ where we have given a second form that will be particularly useful for this work. In general, R is a matrix over all twoparticle degrees of freedom, but in the case considered it reduces to a simple derivative of the functions shown. Before defining W μ L;df ðP f ; P i ; LÞ, we need to introduce a second L-dependent kinematic function, G μ 1 ÁÁÁμ n , first introduced in Refs. [7,8]. For the l ¼ 0 truncation it takes the form G μ 1 ÁÁÁμ n ðP f ; P i ; LÞ ¼ 1 In this work we will specifically need the scalar and vector G functions, denoted G and G μ , respectively. These are defined by keeping zero or one factor, respectively, of k μ in the numerator of the integrand. With these in hand, W μ L;df can be defined via its relation to W μ df as follows: Here fðQ 2 Þ is the form factor of the charged particle while the form factor of the neutral particle, which vanishes identically at Q 2 ¼ 0, is assumed negligible for all values of momentum transfer.

IV. MATRIX ELEMENTS OF THE CONSERVED VECTOR CURRENT
Having introduced the general formalism, we proceed to perform the checks outlined in the Introduction. The first check is to show that, for any finite-volume state, the matrix element with respect to the charge operator Here we have also adopted the shorthand GðP; LÞ ≡ GðP; P; LÞ; i.e., we do not repeat the total momentum argument when it is the same for the incoming and outgoing states. Note that, in this subsection, we are considering not only the finitevolume bound state but also excited states. We do continue to restrict attention to the S-wave only. Substituting this result into Eq. (23), and also taking the relation between W μ df and F μ [Eq. (18)], we find L 3 hP n;f ; LjJ μ jP n;i ; Li This result will prove very powerful in the following derivations. To see the consequences of this for the charge operator we set μ ¼ 0 in the vector current and also set the initial and final states to coincide. This yields hP n ; LjQjP n ; Li ¼ where we have used the x independence of the matrix element to replace L 3 J 0 ð0Þ →Q and have defined F 0 ðPÞ ≡ F 0 ðP; PÞ as a convenient shorthand for systems with identical initial and final momenta. This can be further simplified via the identity which immediately follows from Eqs. (18) and (20). Substituting this into the numerator of Eq. (31) and also using fð0Þ ¼ Q 0 , we recover a very satisfying cancellation of all terms to deduce hP n ; LjQjP n ; Li as expected. This is a highly nontrivial verification that the general 2 þ J → 2 finite-volume formalism is consistent with the consequences of current conservation. The derivation relies on two unexpected identities: First, the fact that the energy derivative of FðP; LÞ can be expressed using the G functions, as shown in Eq. (29), and second, that the Ward-Takhashi identity relates F 0 ðPÞ to the scattering amplitude, Eq. (32).

V. BOUND STATE IN A FINITE VOLUME
We now turn to the implications of the general formalism for bound-state matrix elements in a finite volume.

A. Volume effects on the energies
We start by reviewing results for finite-volume effects in the energy level, E P B ðLÞ, defined to coincide with the moving bound state in the infinite-volume limit, Boosting these energies to the rest frame, we also define Note that the finitevolume energies depend on P, even after boosting back to the rest frame. In the following, we give expressions for the volume-induced shift, δs P B ðLÞ, for two values of total momentum. This represents a small subset of the more general expressions derived in Ref. [18]. 7 The quantization condition, Eq. (21), is satisfied only at the finite-volume energies, e.g., at P B ðLÞ ≡ ðE P B ðLÞ; PÞ. We are thus strictly interested in FðP; LÞ only when it is evaluated at these points. However, taking δs B as a small parameter, we note where P B ≡ ðE P B ; PÞ is the infinite-volume bound-state momentum in a moving frame.
As is discussed in detail in Ref. [18] and reviewed in Appendix C, the subthreshold L dependence of the F function is governed by the binding momentum: In particular, from Eqs. (C3) and (C16) we find where and γ ¼ E P B =M B . This result is to be combined with the inverse scattering amplitude, also evaluated at s P B ðLÞ, but then expanded in powers of δs B to yield where we have used M −1 ðs B Þ ¼ 0. Combining Eqs. (37) and (40) then yields the elegant result which shows that the leading shift to the finite-volume bound state is given directly by the F function, evaluated at the infinite-volume bound-state energy.
To close this section we think it useful to unpack Eq. (41) for two specific cases. First, in the case of vanishing momentum in the finite-volume frame, the three universal orders are given by At Oðe −2κ B L Þ higher derivatives of the inverse amplitude enter, requiring information beyond the coupling, g. For nonzero momenta, the expressions are complicated by the relation between m 0 and m, and by the volume dependence entering γ through To compare these results to those in Ref. [18] we note that, in the earlier work, the authors introduce an L-dependent binding momentum, defined via κ B ðLÞ 2 ¼ m 2 − s B ðLÞ=4. Then the finite-volume shift, δκ B ðLÞ ≡ κ B ðLÞ − κ B , satisfies the relation Combining this with the relation between the coupling and the scattering phase yields Eq. (9) of Ref. [18].
In closing we comment that, due to the reduction of rotational symmetry, higher partial waves do induce finite-volume corrections to the scalar bound state and 7 Related results for bound states can be found in Refs. [95,96,99,100]. corresponding matrix elements. In particular, for P ¼ 0, l ¼ 0 mixes with l ¼ 4; 6; …, as can be seen by the fact that the corresponding off-diagonal components of F are nonzero. These additional angular-momentum contributions are, in fact, not volume-suppressed relative to the S-wave contributions, but are suppressed by powers of the binding momentum in units of the scattering-length analogs appearing in higher-partial waves. For example, the l ¼ 4 phase shift satisfies an expansion analogous to Eq. (6), where M 4 has units of energy. In the case of zero spatial momentum in the finite-volume frame, one can show that the first non-S-wave contribution to s ½000 B ðLÞ is suppressed relative to the leading shift by a factor of κ 8 B M B =M 9 4 . Having reproduced the known expansion for the binding energy [18], we now turn to the finite-L corrections of the bound-state matrix element.

Matrix elements in the L → ∞ limit
We begin by confirming that, in the L → ∞ limit, the finite-volume bound-state matrix element (as described by the general 2 þ J → 2 formalism) coincides with its infinite-volume counterpart. Here it is important to stress that the various quantities we consider have a well-defined L → ∞ limit, only because we are considering them at subthreshold kinematics and thus away from a set of finitevolume poles that becomes arbitrarily dense.
We begin with Eq. (27), the relation between W μ L;df and W μ df . For L → ∞, these two quantities coincide because the G function defining their difference vanishes. This is the case because the sum within G [cf. Eq. (26)] is transformed to an integral in the limit and is exactly canceled by the second, subtracted integral. Equation (23)  The next step is to expand R, evaluated at the finitevolume bound-state energy, about large L. Using the form given by Eq. (25), one readily finds We are now in position to evaluate the limit. The only subtlety is that a double zero, arising from the Lellouch-Lüscher factors, is exactly canceled by the double pole in W μ df . Substituting Eqs. (17) and (49) into Eq. (47), we reach This is exactly the desired result, with the extra factors on the left-hand side accounting for the different normalization conventions of finite-and infinite-volume states.
In this derivation we did not make reference to the Lorentz structure of the current, only to the fact that the 2 þ J → 2 amplitude, W, must have a double pole structure associated with the initial and final bound states. As a result, in general the formalism fulfills the expectation that for an arbitrary current J μ 1 ÁÁÁμ n

Large L expansion of the bound-state matrix element
As shown in Sec. IV, the conserved vector current leads to volume-independent matrix elements at zero momentum transfer. Thus, to reach an interesting large-L expansion, in this section we turn to a scalar current J and define where the subscript indicates that this matrix element defines the scalar charge of the bound state. The infinite-volume bound-state scalar charge is recovered in the L → ∞ limit, i.e., g S;B ≡ lim L→∞ g P S;B ðLÞ. In direct analogy to Eq. (31) above, we observe Here F ðsÞ ≡ M −2 ðsÞW df ðP; PÞ with W df given by Eq. (13), in which the vector current is replaced by a scalar. Note that the numerator includes only the scalar G function, reflecting the scalar current considered. However, the denominator remains identical to the vector case since the Lellouch-Lüscher factors are independent of the current. We have also introduced g S as the scalar charge of the single-particle state, g S ≡ fð0Þ, where f is the single particle form factor fðQ 2 Þ ≡ hP f ; g S jJ jP i ; g S i: ð54Þ As above, we take the coupling of the current to the other constituent particle to be negligible.
With these ingredients in hand it is straightforward to expand Eq. (53) about L → ∞ to reach We note that a great deal of structural information enters the leading finite-volume correction. The δs B ðLÞ-dependent term is the correction induced from the energy shift and is thus proportional to energy derivatives of both the inverse amplitude and the 2 þ J → 2 transition amplitude [entering via F ðsÞ]. The second term in Eq. (55) arises due to a mismatch between the scalar charge of the bound state and the summed charges of its constituents. The final term in Eq. (55) is a direct consequence of the triangle diagram, Fig. 2(c). We close this section with a final, more explicit result for the leading-volume correction in the case where the CMF and finite-volume frames coincide, i.e., P ¼ 0. Substituting the leading result for δs ½000 B , and results from Appendix C for the G functions, one finds The leading 1 in the square brackets arises from the triangle diagram, Fig. 2(c), and will be the dominant finite-volume effect provided jg S − g S;B j ≤ jg S;B j=2.

C. Numerical expectations for finite-volume dependence
In this section, we use the full 2 þ J → 2 formalism to explore the finite-volume corrections to a bound-state matrix element in an example with scattering parameters chosen to mimic the deuteron. As above we consider the simplest case of a scalar current and an S-wave bound state and examine the finite-volume corrections given by Eqs. (53) and (55). For the 2 → 2 scattering amplitude, we use the phenomenological values for the pn scattering length (a ¼ 5.425 fm) and effective range (r ¼ 1.749 fm) to describe the scattering amplitude and spectrum. With the nucleon mass at m ¼ 934 MeV, 8 Fig. 3 as a function of q ⋆2 and the location of the bound state is indicated. We make two assumptions to simplify the numerical exercise: First, we assume that the infinite-volume boundstate form factor is constant, i.e., F B ðsÞ ¼ g S;B . As seen in the fourth term of the brackets in Eq. (56), this contribution is suppressed by 1=L, and thus it is reasonable that this approximation will not strongly alter the prediction. Second, we assume that difference g S − g S;B is numerically small and set g S ¼ g S;B .
Within this setup one can numerically evaluate Eqs. (53) and (55) and compare the results. The first step is to determine the finite-volume bound-state energy, using the effective-range description of the pn scattering amplitude in the quantization condition, Eq. (21). Figure 4(a) shows the bound-state energy as a function of L for both jdj ¼ 0 and 1. The solid lines represents the full solution obtained from Eqs. (3) and (21) with the pn-scattering parameters. The dashed lines correspond to the leading-order approximation using Eq. (41) for the same momenta. These results reproduce those of Ref. [18], and we see that for lattice calculations performed at m π L ∼ 4, deviations between the exact and approximated forms are significant.
Turning to the two-particle matrix elements, Fig. 4(b) shows the ratio of the finite-volume bound-state matrix element to the infinite-volume scalar charge, g P S;B ðLÞ=g S;B . Solid lines represent the full solution, using Eq. (53) evaluated at the finite-volume energy, and the dashed lines give the leading-order shift of Eq. (55). Again, significant deviations arise between the full prediction, the leadingorder expansion, and the infinite-volume result. This illustrates that, to reliably extract infinite-volume matrix FIG. 3. Plot of q ⋆ cot δðq ⋆ Þ and − ffiffiffiffiffiffiffiffiffiffi −q ⋆2 p as a function of q ⋆2 , in units of MeV 2 , for the effective range expansion, Eq. (6), using the scattering length and effective range for the pn system in 3 S 1 . The vertical dashed line indicates the deuteron, with binding momentum κ B ∼ 45.58 MeV. 8 We work here with isospin symmetry, approximating m p ¼ m n . elements of shallow bound states such as the deuteron, it is highly beneficial to use the full formalism which removes an infinite series of terms scaling as powers of e −κ B L . In the present example, only at m π L ∼ 8 do corrections scale to the percent level.

VI. CONCLUSION
In this work, we have provided strong consistency checks on, and also explored various consequences of, the formalism derived in Refs. [7,8], which gives a relation between finite-volume matrix elements, schematically denoted h2jJ j2i L , and the corresponding infinite-volume 2 þ J → 2 amplitudes.
First, in the case of the conserved vector current, we have shown that the resulting prediction for the two-particle matrix element of the charge operator, h2jQj2i L , behaves as expected. Specifically, the matrix element is L-independent and equal to the sum of the constituent charges. Though it is clear that this relation must hold, the way it arises in the mapping is highly nontrivial, relying on an identity relating various L-dependent geometric functions [Eq. (29)] as well as a relation between the 2 → 2 and 2 þ J → 2 amplitudes that follows from the Ward-Takahashi identity [Eq. (20)].
Second, for a generic local current, we have demonstrated that the mapping of Refs. [7,8] reproduces the expected behavior in the case of an S-wave two-particle bound state. By analytically continuing the formal relations below threshold to the bound-state pole, we have confirmed that the finite-and infinite-volume matrix elements are equal up to volume corrections scaling as e −κ B L , where κ B is the binding momentum of the state. This is an expected extension of the well-known result for the L dependence of the bound-state energy.
These two checks give confidence that our admittedly complicated formalism correctly describes two-particle finite-volume states and is ready to be implemented in a LQCD calculation, with the first application likely being the ðππÞ I¼1 þ J μ → ðππÞ I¼1 transition amplitude, allowing one to extract the electromagnetic form factors of the ρ.
As an additional example of the utility of the general approach, we have determined the full functional form of the leading, Oðe −κ B L Þ correction to the bound-state matrix element of a local scalar current. The result, Eq. (56), shows that the coefficient of the leading exponential depends on the bound state's coupling to the two-particle asymptotic state, on the scalar charges of both the bound state and its constituents, and also on derivatives of both the 2 → 2 scattering amplitude and the bound-state form factor.
While the structure of this relatively simple prediction is instructive, we stress that in practice it is more useful to use the general relation of Refs. [7,8] to extract the 2 þ J → 2 over a range of energies, including in a neighborhood around the bound-state pole. Doing so removes an infinite series of terms scaling as powers of e −κ B L and, for shallow bound states, allows one to control an otherwise dominant source of systematic uncertainty. To stress this point, as a final exercise, we have presented numerical comparisons of the leading e −κ B L correction with the full finite-volume shift, for a toy setup mimicking a LQCD calculation of the deuteron's scalar charge. For physical pion masses and volumes in the range m π L ∼ 4 to 7, we find that the finite-L correction will dominate the infinite-volume charge and that removing only the leading exponential also does not give a reliable extraction. Thus, we conclude that the full method must be used to gain a reliable result for the form factors of shallow bound states as well as resonances.
For details on two-hadron form factors for general scattering systems, we refer the reader to a future article which will discuss the analytic structure of these amplitudes [101]. An additional check is underway to reproduce analytic expressions for the 1=L expansion presented in Ref. [102], for the threshold-state matrix element of a scalar current in a weakly coupled system. In this appendix we demonstrate that, in theories with a two-particle bound state, the 2 þ J → 2 amplitude satisfies Eq. (12), repeated here for convenience, To show this, it is necessary to return to the matrix element definition of the amplitude, W μ , given in Eq. (9). Inserting a complete set of states on either side of the current, J μ , we reach Z d 3 P 00 ð2πÞ 3 2E P 00 where we have kept only the bound-state sector of the Fock space, as this will be sufficient to identify the pole that we are after. Three additional subtleties arise here: (1) To properly implement the normalization of the bound state, we must integrate over all spatial momenta with the standard Lorentz-invariant factor as shown. (2) Since the spectral decomposition can only be performed on the full matrix element, we have dropped the "conn" subscript that appears in Eq. (9). To preserve the definition we have included the ellipsis on the right-hand side, which is understood to represent all disconnected contributions. These will, however, play no role, since they do not contain the bound-state pole. (3) The expression we are after requires the analytic continuation of P f and P i to the subthreshold region. This is subtle at the level of Fock states and is more easily understood by rewriting the result in terms of operators projected to definite momentum. This, in turn, reveals that the time ordering of the operators must be carefully treated, as we explain in more detail below. The next step is to substitute hP f ;k ⋆ f ; outjP 00 ; Bi ≡ ð2πÞ 4 δ 4 ðP f − P 00 B Þig; ðA4Þ Here the four-dimensional delta function arises in direct analog to the standard relation between the T matrix and scattering amplitude and leads to the definition of the bound-state coupling g. Using the spatial delta functions to evaluate the integrals in Eq. (A2), we reach where it is understood that one must set P 00 → P f and P 00 → P i . Here we have also written the remaining temporal delta functions as integrals over time.
Introducing the integrals over x 00 0 and x 0 0 allows us to address the subtlety mentioned as point (3) above. Studying the correlation functions reveals that the above expression does not correctly treat all time orderings. For the present case, this is resolved by restricting the integral over x 00 0 from 0 to ∞ and similarly that over x 0 0 from −∞ to 0. Doing so, and also including the iϵ prescription required to project the external states in the correlator to the vacuum, one can evaluate both integrals to reach This is the result that we had aimed to prove. Up to the Oððs i;f − s B Þ 0 Þ terms that we neglect, one can replace each pole with the covariant form and also drop the ellipses.
In this appendix, we demonstrate how Eq. (20) follows from the Ward-Takahashi identity. A consequence of current conservation, the Ward-Takahashi identity relates a given n-point Green function, coupled to an external conserved current, to the corresponding (n − 1)-point Green function in which the current is omitted. Let C μ be a five-point function coupling the conserved vector current, J μ , to two neutral and two charged mesons. The Ward-Takahashi identity then reads q μ C μ ðp 0 ; k 0 ; p; kÞ ¼ Q 0 ½Cðp 0 þ q; k 0 ; p; kÞ where q μ ¼ ðp 0 þ k 0 Þ μ − ðp þ kÞ μ ¼ P 0μ − P μ , with the second equality introducing notation for the total momenta of the outgoing and incoming two-meson states. We have also introduced C (with no index) as the four-point function without the current insertion. We further define k and k 0 as the initial-and final-state momenta of the neutral particles, respectively, and p ¼ P − k and p 0 ¼ P 0 − k 0 as the corresponding momenta for the particles carrying the charge, Q 0 . The 2 þ J → 2 and 2 → 2 amplitudes, W μ and M, respectively, are related to the Green functions by amputating the external meson propagators and placing them on the mass shell, i.e., where the amputated Green functions are defined as where the ratios of amputation factors arise since the Ward-Takahashi identity changes the momenta carried by the mesons on the two sides of the equation.
In the limit where p 0 , k 0 , p, and k go on shell, the numerators on the right-hand side of Eq. (B8) vanish but the denominators do not, yielding the well-known Ward identity: q μ W μ ¼ 0. In addition, the long-range pieces that define the difference between W μ and W μ df [see Fig. 2(b)] are proportional to ðP 0 þ PÞ μ and therefore also vanish when contracted with q μ . (Equivalently they are proportional to the single-particle matrix element of J μ and must therefore also satisfy the Ward identity.) It follows that W μ df itself satisfies the identity: q μ W μ df ¼ 0. Returning to the off-shell relation, Eq. (B8), we reexpress all functions in terms of P; P 0 ; k, and k 0 to write q μ C μ amp ðP 0 − k 0 ; k 0 ; P − k; kÞ Applying a P 0 ν derivative on the left-hand side then gives and, applying the same to the right-hand side, one finds Next, before equating the two sides, we take the zero-momentum-transfer limit (P 0 → P) and substitute for the 1 þ J → 1 matrix element at zero momentum transfer. This then gives Setting P ¼ P 0 has greatly simplified the expressions, but care must be taken as the second and third terms on the right-hand side will diverge when we set p 0 , k 0 , p, and k to their on-shell values. Indeed, these are the same divergences that appear in the difference between W μ and W μ df , with the only subtlety that they were first defined in on-shell amplitudes at P 0 − P ≠ 0. Fortunately, in the present case the distinction is unimportant because, when applied to the divergence-free amplitude, the zero-momentum-transfer and on-shell limits commute. We can thus move the second and third terms to the left-hand side and take p 0 , k 0 , p; k on shell to conclude This remarkable result gives a clear interpretation to W μ df in the forward limit. Finally, since the derivative is with respect to total momenta, we can easily project both sides to definite angular momentum. This leads to as claimed in Eq. (20) for the special case of S-wave systems.

APPENDIX C: ANALYTIC CONTINUATION OF FINITE-VOLUME FUNCTIONS BELOW THRESHOLD
In this section we give results for the analytic continuations of the F and G functions below threshold. Specifically we require results for FðP; LÞ, GðP; LÞ, and G μ¼0 ðP; LÞ, where we recall that a single momentum argument within G indicates that the initial-and final-state four-momenta are equal. Each of these can be written in terms of a class of functions naturally extending those defined in Refs. [13][14][15][16] The relations to the finite-volume functions that we require are then given by where the last result also assumes that P is parallel to thê z axis. For P 2 < ð2mÞ 2 , the summand of c ðnÞ is a smooth function of k ⋆ with a finite region of analyticity. As a result, the sum and integral must become exponentially close to each other, with the scale in the exponential given by the grid spacing of the sum (set by L) and the size of the analytic domain (set by 4m 2 − P 2 ). To make this explicit, we apply the Poisson summation formula to c ðnÞ JM , evaluated at a generic subthreshold four-momentum, P κ , satisfying m 2 − P 2 κ =4 ¼ κ 2 . We find 9 The c For the analytic work presented here, the dimensionful versions prove slightly more convenient.
where we have used the fact that the integration measure is a Lorentz invariant, d 3 k=ω k ¼ d 3 k ⋆ =ω ⋆ k . The kinematic variables in the CMF are related to the moving frame variables via standard Lorentz transformations, where k ⊥ ¼ k − k jj , k jj ¼ ðk ·βÞβ, β ¼ P=E is the velocity, and γ ¼ E=E ⋆ is the Lorentz factor. We can then write the phase factor in terms of the CMF momenta, with m 0 defined in Eq. (38). With these relations in hand we can write the integrand solely in terms of k ⋆ , Next, we evaluate the angular piece by introducing spherical Bessel functions and making use of the standard plane wave expansion, where j l ðzÞ is the spherical Bessel function of the first kind. The angular integral in Eq. (C9) becomes For the cases considered here we require only I 00 ðLk ⋆ jm 0 j; mÞ ¼ 4π sin ðLk ⋆ jm 0 jÞ Lk ⋆ jm 0 j ; ðC13Þ I 10 ðLk ⋆ jm 0 j; mÞ ¼ i4π ffiffi ffi 3 p sin ðLk ⋆ jm 0 jÞ ðLk ⋆ jm 0 jÞ 2 − cos ðLk ⋆ jm 0 jÞ Lk ⋆ jm 0 j m 0 · P jm 0 jjPj : To evaluate the remaining integral over k ⋆ , we express the sinusoidal functions in I JM in terms of exponentials and then divide the function I JM into two terms, denoted I ðþÞ JM and I ð−Þ JM : I ðþÞ JM is defined by replacing the sinusoidal functions with the part of their exponential representation that decays as k ⋆ → i∞, e.g., sinðxÞ → e ix =ð2iÞ, and I ð−Þ JM is defined in the same way for the part that decays as k ⋆ → −i∞, e.g., sinðxÞ → −e −ix =ð2iÞ. This leads to a decomposition of c In addition to the k ⋆ ¼ AEiκ pole, the integrand has branch cuts starting at k ⋆ ¼ AEim, associated with the square root in ω k . These lead to exponential corrections of the order of Oðe −mL Þ that are ignored throughout, i.p. already in deriving the formalism considered in this work. Therefore these contributions should also be dropped in the present evaluations. Keeping only the contribution from the κ pole, we deduce When P ¼ 0, these expressions simplify significantly: