Prospects for a lattice calculation of the rare decay $\Sigma^+\to p\ell^+\ell^-$

We present a strategy for calculating the rare decay of a $\Sigma^+ (uus)$ baryon to a proton $(uud)$ and di-lepton pair using lattice QCD. To determine this observable one needs to numerically evaluate baryonic two-, three-, and four-point correlation functions related to the target process. In particular, the four-point function arises from the insertion of incoming and outgoing baryons, together with a weak Hamiltonian mediating the $s \to d$ transition and an electromagnetic current creating the outgoing leptons. As is described in previous work in other contexts, this four-point function has a highly non-trivial relation to the physical observable, due to nucleon and nucleon-pion intermediate states. These lead to growing Euclidean time dependence and, in the case of the nucleon-pion states, to power-like volume effects. We discuss how to treat these issues in the context of the $\Sigma^+\rightarrow p\ell^+\ell^-$ decay and, in particular, detail the relation between the finite-volume estimator and the physical, complex-valued amplitude. In doing so, we also make connections between various approaches in the literature.


Introduction
The transition of an sto a d-quark (s → d) requires a flavour changing neutral current, which is only allowed through quantum corrections within the Standard Model of particle physics. Consequently, processes involving such transitions are rare in the Standard Model and could be enhanced by potential new physics that includes a flavour changing neutral current in its Lagrangian. An example for such a process is the rare semi-leptonic hyperon decay Σ + → p + − , for which the muonic mode has been recently measured by the LHCb experiment [1] with a branching ratio of Evidence for this decay had previously been found by the HyperCP collaboration [2] giving B(Σ + → pµ + µ − ) = (8.6 where the first uncertainties are statistical, and the second is systematic. This determination follows from three events, all at nearly the same invariant mass of the µ + µ − pair. However, such a resonant structure in the µ + µ − invariant mass could not be confirmed by the more recent LHCb measurement [1]. The current state-of-the-art Standard Model theory predictions for this process [3][4][5] use a combination of dispersion relations, experimental input, various formulations of Baryon Chiral Perturbation Theory and model estimates (e.g. vector meson dominance) and arrive at a range of results [5] 1.6 × 10 −8 ≤ B Σ + → pµ + µ − ≤ 8.9 × 10 −8 .
(1. 3) More details on the phenomenological background of this decay can be found in section 2. The rare hyperon decay Σ + → p + − can be viewed as the baryonic analogue of the rare kaon decay K → π + − , which has been previously calculated from first principles using lattice simulations by the RBC/UKQCD collaboration, including recent results at physical quark masses [6][7][8]. Taking inspiration from this progress, in this paper we explore prospects for calculating the required form factors for the Σ + → p + − decay using lattice QCD. While other hyperon decays with much higher yields, such as Σ − → n − ν , can be used to make measurements of the CKM matrix element V us in order to test for new physics that would break the CKM unitarity relations, this decay can be sensitive to new physics due to its rarity within the Standard Model.
As already discussed in refs. [6][7][8][9][10], a key challenge in extracting decays such as Σ + → p + − and K → π + − from lattice QCD is that the physical observables (most directly defined in terms of infinite-volume Minkowski-signature correlation functions) contain onshell intermediate states that can propagate between the weak Hamiltonian, effecting the s → d transition, and the electromagnetic current emitting the di-lepton pair. For the case of the rare hyperon decay in particular, intermediate N π states contribute, where N represents the nucleon doublet and π the pion triplet. As we discuss in more detail in the following sections, in practice the states are projected to definite isospin as this is a good quantum number of the numerical calculation, provided the light quarks are degenerate and dynamical electromagnetism is not included, as we assume throughout. Fortunately, three-(or-more)-particle intermediate states are kinematically guaranteed to be off-shell and do not require special treatment.
Complications arise because numerical lattice QCD calculations only allow one to directly determine Euclidean correlation functions in a finite spacetime volume. Specifically, the finite-volume Euclidean correlator that most closely matches the rare hyperon decay is a four-point function, defined with operators to create the incoming Σ + and the outgoing p as well as the weak Hamiltonian and electromagnetic current. Careful examination of this correlator shows that the on-shell N π intermediate states manifest in a number of ways, all of which complicate the calculation.
First, after the baryonic operators are used to project out the Σ + and p states, one finds that the on-shell intermediate states lead to exponentials that grow with the Euclidean-time separation between the weak Hamiltonian and the current [9]. In practice, the number of such exponentials is finite, dictated by the discrete finite-volume spectrum, and thus these states can be removed through various strategies that we detail in the following sections.
However, a consequence of discarding these terms is that the resulting finite-volume estimator has poles at the locations of all finite-volume energies with N π quantum numbers. In addition, the resulting quantity is known to have power-like volume effects away from the poles, and to miss the imaginary part appearing in the physical amplitude due to the long-distance propagation of intermediate states. In short, removing growing exponentials defines a finite-volume object that a priori has no clear relation to the targeted amplitude. Fortunately, as we describe in section 4 following refs. [9,10], the strategy to convert the finite-volume estimator to the physical observable (and thereby cancel the poles and include the imaginary contribution) is known and can be applied in this case.
The remainder of this paper is organized as follows: In section 2 we discuss the currently available phenomenological strategy and predictions to compute the Σ + → p + − branching ratio. Section 3 then outlines our strategy on the lattice, which aims to recover the Σ + → p + − amplitude via carefully chosen, numerically calculable Euclidean correlation functions. Here we discuss various strategies to remove the exponentially growing terms that will appear in the direct lattice result. In section 4 we discuss the removal of finite-volume singularities and make contact with the physical observable. In addition to translating the general formalism of refs. [9,10] to our particular case, we also provide an explicit expression for an expansion that arises when the volume is tuned so that the mass of the Σ + coincides with one of the finite-volume N π energies. We close with a brief summary in section 5. This work also contains five appendices detailing various technical aspects useful for the practical calculation.

Phenomenological background
In the following, we will briefly review the phenomenological determination [3][4][5] of the branching ratio for Σ + → p + − leading to the result quoted in eq. (1.3). Short distance contributions to Σ + → p + − originate from penguin and box diagrams (cf. figure 1) and are found to contribute only at the order of 10 −12 [3] to the branching ratio of the muonic mode Σ + → pµ + µ − , which is much smaller than the experimental measurements (eqs. (1.1) and (1.2)), indicating that this decay is long-distance dominated. Figure 1: Short-distance Standard Model contributions to the s → d transition from penguin and box diagrams.
The matrix element for the long-distance Standard Model contribution to the Σ + → p + − decay can be written as [3,11] in terms of the four form factors a(q 2 ), b(q 2 ), c(q 2 ) and d(q 2 ). Here γ ν , γ 5 and σ µν are (combinations of) Minkowski gamma matrices, with conventions defined in appendix B, and u Σ (k), u p (p), u (p − ), and v (p + ) are the usual Dirac spinors for the incoming Σ + and the outgoing proton, − , and + , respectively. The four-momentum transfer is q = k − p, where k and p are the on-shell four-momenta of the Σ + and proton, respectively. Some information about the form factors a and b can be obtained from the decay Σ + → pγ with a real photon. The respective decay rate 1 can be written as where |q| is the energy of the photon and θ is the angle between the spin of the Σ + and the momentum of the proton. The imaginary parts of the four form factors can be obtained from unitarity using amplitudes for Σ → N π and N π → N γ * . While the amplitude for Σ → N π is known from experimental measurements [12], the authors of ref. [3] calculate the amplitude N π → N γ * from Chiral Perturbation Theory (ChiPT) using either the relativistic baryon ChiPT [13] or the heavy-baryon ChiPT [14,15] formulation. The momentum dependence of the imaginary parts of the form factors is found to be very mild. Once the imaginary parts of the form factors are known, information on the real parts of a(q 2 ) and b(q 2 ) at q 2 = 0 can be obtained from equations (2.2) and (2.3) and experimental data for the decay Σ + → pγ. Since this decay only determines values for |a(0)| 2 +|b(0)| 2 and Re[a(0)b(0) * ], this leads to four possible solutions for Re a(0) , Re b(0) . Motivated by the mild q 2 -dependence of the imaginary parts of the form factors, the authors of ref. [3] assume that Re a(q 2 ) = Re a(0) and Re b(q 2 ) = Re b(0) for their prediction of B(Σ + → p + − ). The real parts of the c(q 2 ) and d(q 2 ) form factors are calculated assuming vector meson dominance in [3] and explicitly calculating vector meson pole contributions to the decay Σ + → p + − . The q 2 -dependence of Re c(q 2 ) and Re d(q 2 ) is found to be mild, just like in the imaginary parts.
Depending on the formulation of baryon ChiPT used and the four possible solutions for Re a(0) and Re b(0) from Σ + → pγ decays, the authors of ref. [3] find the Standard Model prediction for B(Σ + → pµ + µ − ) to be in the range and very similar ranges are found in refs. [4,5].

Extracting the amplitude from Euclidean correlators
In this section and the next, we describe how to extract the Σ + → p + − amplitude from a numerical lattice calculation. The approach closely follows the methods of refs. [6][7][8][9][10], adjusted here to treat issues specific to this system. This section details the Euclidean two-, three-, and four-point correlation functions needed to construct a finite-volume estimator, denoted by F µ (k, p) L . The following section describes how to relate this quantity to the physical rare hyperon amplitude.

Spectral representation
The determination of the long-distance contribution to Σ + → p + − requires a calculation of the Σ + → pγ * amplitude, defined as with r and s labelling the azimuthal spin component of the state. Here we are assuming Minkowski-signature conventions and working in an infinite space-time volume. This amplitude can be re-expressed as a Dirac matrix, A µ (k, p), using the relation with the spinors u p and u Σ of the proton and Σ + , respectively. 2 The effective weak-Hamiltonian density of the qs → qd transition is given by [16] H where the C i are Wilson coefficients, the Q q i are four-quark operators, defined in terms of Dirac spinors for up, down, strange and charm quarks (respectively u, d, s and c) as and it is understood that the two are renormalized in the same scheme, at the same scale, such that the weak Hamiltonian is scheme independent. Here i, j denote colour indices and we define γ L µ ≡ γ µ (1 − γ 5 ). There are additional four-quark operators in (3.3) with Wilson coefficients of order | VtsV td VusV ud | 0.00142 which will be neglected in this work. The electromagnetic current in eq. (3.1) is given by and we make use of translational invariance by fixing the position of the electromagnetic current to y = 0. Including an additional Fourier transform on the current would lead to an overall momentum-conserving Dirac delta function, to be removed at a later step, and we find it more convenient to follow the approach where this is never introduced. The Hamiltonian density decomposes into a parity-positive and a parity-negative component defined via the parity operatorP according tô where P µ ν = diag 1, −1, −1, −1 andP is the Hilbert-space representation of the parity operator. Both parity sectors contribute to the amplitude we are evaluating. Defining A ± µ (k, p) as in eqs. (3.1) and (3.2), but with H W (x) replaced by H ± W (x), one can next decompose each definite-parity amplitude in terms of form factors as follows [3,11]: where we recall that q = k − p is the four-momentum transfer of the virtual photon. This form-factor decomposition is derived in Appendix A. Note also that, while the amplitude is a Dirac matrix and thus depends on individual components of the four-momenta, the form factors are Lorentz scalars and can therefore only depend on q 2 . We will see in the following that the amplitude, and thus also the form factors, are complex-valued due to on-shell intermediate N π states. Since the Euclidean correlators are real-valued, this complexity already signals the fact that it is non-trivial to extract the amplitudes. This turns out to be closely related to the interplay of the Euclidean signature and the finite volume. As we will show in the following, the quantum numbers of the contributing N π states differ for A + µ and A − µ , and thus the finite-volume formalism must be applied independently to the two quantities.
To explain this in more detail we return to eq. (3.1) and insert a complete set of states between the current and the weak Hamiltonian to write 10) where we have introduced the spectral functions, satisfying Note that one must treat the two time orderings separately and this leads to two types of intermediate states encoded in ρ and σ, which have strangeness S = 0 and S = −1 respectively. The sums over α and β represent both sums and phase-space integrals over the multi-particle QCD Fock space for all states that contribute. For example, the sum over α includes N π, N ππ, ∆π and ΛK states. Evaluating the time integrals then gives a compact result The aim of the following sections is to review how this amplitude can be extracted from finite-volume Euclidean-signature correlation functions.

Euclidean correlation functions
We now discuss how to extract a finite-volume estimator for the desired Minkowski-space amplitude (3.13) from Euclidean correlation functions that can be calculated on the lattice. All quantities in this section (e.g. Dirac γ-matrices, four-vectors) are defined with Euclidean conventions detailed in Appendix B.

Two-point functions
The two-point function of a baryon B can be written as where ψ B β (t, x) and ψ B α (t, x) create and annihilate, respectively, a baryon B and α and β are Dirac spinor indices. Examples for operators that have overlap with the proton p and Σ + are with Γ A = 1, Γ B = Cγ 5 for spin 1/2 particles and the charge conjugation operator defined by C = γ 0 γ 2 . P + = 1 2 (1+γ 0 ) projects to the positive parity state. Here, and below, Roman and Greek indices refer to colour and Dirac indices, respectively.
For t > 0, the two-point function (3.14) has the spectral representation plus contributions from excited states, which are exponentially suppressed in t by their higher energies. The moving-frame energy of the baryon is given by E B (p) = M 2 B + p 2 and the four-momentum vector is p B ≡ (iE B (p), p). We have introduced the projector The overlap factor Z B is defined by (3.19) Throughout this section and the next, we neglect the finite-volume effects on single-hadron energies and matrix elements. These are known to be exponentially suppressed, scaling as e −MπL . We also include the ∞ label on states that are normalized according to 20) as this differs from the normalization of finite-volume states used below.

Three-point functions
We turn now to the three-point function of the effective weak Hamiltonian H W between a Σ + and a p operator 21) with the effective weak-Hamiltonian density H W given in eq. (3.3). We leave the parity labels implicit throughout this section unless stated explicitly.
Here the weak-Hamiltonian is assumed to be appropriately renormalised. In the case of chirally symmetric fermion discretisations and a massless renormalisation scheme, the 4-quark operators in (3.4) are protected from mixing with other dimension-6 operators. In fact, the V − A structure of our particular 4-quark operator guarantees that all divergences are removed with massless renormalisation factors, even when the matrix elements themselves are evaluated at non-zero quark mass. This is described in detail, in the context of domain-wall-fermion lattice calculations, in refs. [17,18], where the K L − K S mass difference is evaluated using the same operator. The key point for the present work is that no aspects of the renormalisation are affected by the external states, so we can directly adopt the strategy of the earlier publications. See also the text in the paragraph after eq. (3.31).
The three-point function (3.21) gives rise to four different topologies 3 for the Wick contractions, which are shown in Figure 2. The double point labelled with H W shows the position of the weak Hamiltonian, and the points labelled Σ + and p are the positions of the Σ + and proton operator, respectively. The quark lines are labelled by their respective quark flavours. The two diagrams shown on the left-hand side of Figure 2 arise from contractions using the Q 1 operator in the weak Hamiltonian, the two diagrams on the right-hand side arise from contractions using the Q 2 operator. Figure 2: The four different topologies for the Wick contractions of the three-point function (3.21). Two fully connected contractions C sd and C su and two topologies (E and S) containing quark loops. The two diagrams on the left (C sd and E) arise from contractions using the Q 1 operator, the two diagrams on the right (C su and S) from Q 2 .
The spectral representation of the Euclidean three-point function Γ (3) is given by where r and s are the spins of the proton and Σ + , respectively. In the following it will be convenient to define the overall normalization factor where B, B ∈ {p, Σ}. This factor can be constructed using information extracted from Σ + and proton two-point functions (cf. eq. (3.17)). Completing the spin sum, the spectral representation (3.22) can be written as (3.26) The spectral representation of such a three-point function is given by Just as with the weak-Hamiltonian operator, the electromagnetic vector current is also assumed to be renormalised. If the conserved Noether current is used, then the vector Ward identity is exactly obeyed on the lattice and no renormalisation is necessary.
Here we have only considered three-point functions with single-hadron states. However, as discussed in the following subsection, to construct the target finite-volume estimator, one may also require matrix elements involving the finite-volume analogue of a multi-particle excited state. Many details of the construction of A rs H and A rs µ,B (q) also apply for the excited-state analogues, but important differences arise. We describe this in detail in Appendix D.
Finally, we define amputated versions of three-point functions with certain factors removed. For the vector current we takê and for the weak Hamiltonian In both cases we drop time dependence which cancels in the limit that the ground-state dominates.

Four-point functions
We turn now to the four-point function of the time-ordered product of the weak-Hamiltonian density, H W (x), and the electromagnetic current, J µ (0), between a Σ + and a p state Here the subscripts L, T indicate that the quantity is evaluated in a finite space-time volume. This is of particular importance for the four-point function, so we emphasize the fact with our notation. Note also that the lattice path integral will always give the time-ordered product of all four-fields. We restrict attention to the case of t Σ < 0, t H < t p so that the fields can be written as shown.
The weak-Hamiltonian and electromagnetic current operators should be renormalised in the same way as the relevant three-point functions above to remove divergences coming from the operators themselves. However, additional divergences can come from the contact of the two operators x = (t H , x) → 0. Refs. [6,19] describe in detail the origin and cancellation of these divergences, which we summarise here. From power counting, the contact term between H W and J µ can diverge at most quadratically. However, use of the conserved Noether current for the electromagnetic operator gives exact QED gauge invariance, which reduces the degree of divergence by two down to a logarithmic one. Finally, this logarithmic divergence is independent of the quark mass, and therefore cancels in the GIM subtraction between the up and charm quarks . Therefore, so long as a conserved electromagnetic vector current is used, there are no divergences as H W approaches J µ . While this previous work has been performed in the context of the rare Kaon decay K → π + − , since none of the arguments rely on the nature of the external states, no additional modifications are required for use in the baryonic decay Σ + → p + − .
As with the Minkowski amplitude and the three-point function considered above, q = k − p denotes the momentum transfer at the electromagnetic vertex. From the four-point function, one obtains six Wick contractions for each of the four topologies (cf. Figure 2) of the three-point function with the weak Hamiltonian: the electromagnetic current J µ can be inserted on any of the five quark lines or on a disconnected quark-loop. The diagrams corresponding to the in total 24 Wick contractions are shown in Appendix E.
We next remove the overlap factor Z pΣ (t p , t Σ ; p, k), given in eq. (3.24), to defineΓ We drop the dependence on T , t p , and t Σ in this quantity since the ratio is independent of these time coordinates as long as T t p , T |t Σ | 1/∆E Σ and |t p | 1/∆E p , where ∆E Σ , ∆E p are the gaps between the ground and first-excited states for the quantum numbers indicated by the subscript. We view it as a task of the numerical analysis to remove or quantify residual dependence on T, t p , and |t Σ | and omit these coordinates for the remainder of this work.
The amputated four-point function is then equal to the matrix element where we have introduced the finite-volume, Euclidean-time-dependent analogue of eq. (3.1) Inserting a complete set of finite-volume states between the current and weak-Hamiltonian density, one can give a spectral representation of the four-point function in Euclidean spacetime aŝ where ρ µ (ω) L and σ µ (ω) L are defined as in eqs. (3.11) and (3.12) above but here with the Euclidean conventions in the gamma matrices and with the sum now running over the discrete finite-volume spectrum. For example, ρ µ (ω) L can be written as and E n , k L is the nth finite-volume state with the relevant quantum numbers to contribute, normalized as E n , k E n , k L = 2E n (k).

Finite-volume estimator for the decay amplitude
Our aim now is to extract the infinite-volume, Minkowski-signature amplitude (with spectral representation given in eq. (3.13)) from the finite-volume Euclidean four-point function, decomposed above in eq. (3.36). To do so, one needs to separately treat the issues of Euclidean time and finite volume, and we find it most instructive to address the first point in this subsection and the second in the section following. To this end we define a physicalenergy finite-volume estimator of eq. (3.13) as follows: .
This definition looks similar to eq. (3.13) but with the key difference that finite-volume spectral functions have been substituted. As a result, the i pole prescription has also been discarded, since this has no effect in the finite volume. To see this more explicitly we write out the first term by substituting the definition of .
The sum over n runs over the discrete set of finite-volume states including the proton-like ground state and multi-hadron excited states that can be related to proton-pion and other scattering amplitudes. A subtlety of this analysis is that, for non-zero k, parity is no longer a good quantum number for the finite-volume multi-particle states. We will avoid this issue by restricting attention to k = 0 in the following section.
Various strategies are possible for extracting this finite-volume estimator from the amputated four-point function. One technical issue affecting all methods is that, becausê Γ (4) µ (t H ; p, k) L only depends on the projected spectral functions, P p (p) σ µ (ω) L P Σ (k) and P p (p) ρ µ (ω) L P Σ (k), it is only possible to extract a similarly projected version of F µ (k, p) L . Rather than carrying these projectors in all subsequent equations, we find it most convenient to change to spin indices at this stage, defining and similar for all other quantities defined as Dirac matrices above.
To extract F rs µ (k, p) L fromΓ (4)rs µ (t H ; p, k) L , a crucial issue that any method must address is that certain intermediate states lead to exponentially growing Euclidean time dependence. This arises because, in eq. (3.32), one is multiplying by growing exponentials depending on the energies of the incoming Σ + and outgoing proton. If the intermediate energies in the sum over n are sufficiently large, then these contribute decaying exponentials that outweigh the growth. However, due to low-lying finite-volume states,Γ To understand the point in more detail, we consider the unphysical quantity Here ω is chosen such that the integral over t H is convergent, and one finds a result that is very similar to the targeted finite-volume estimator, F rs µ (k, p) L . In fact, the right-hand side gives this desired quantity in the ω → 0 limit, but this is not useful as the integral on the left-hand side is divergent if evaluated at ω = 0. Physically, this expression corresponds to allowing the weak Hamiltonian to carry away energy from the system such that no on-shell intermediate states occur. It thus solves the problem of growing exponentials, but at the unacceptable cost of giving an unphysical quantity.
Two closely related options are available to reach the desired expression at ω = 0. The original proposal, introduced in refs. [6,9], is to integrate t H over a finite range of times, i.e. over the range t H ∈ [−T a , T b ]. One can then remove growing exponentials as a function of T a and T b in order to extract the desired finite-volume quantity. A closely related alternative, described in ref. [10], is to explicitly remove the exponentials as a function of t H before integrating and then to re-introduce the missing poles in a second step.
Here we focus on the method of refs. [6,9], defining The summed correlator has a spectral representation given by in which the growing exponentials are displayed explicitly. From these expressions one sees that contributions growing as T a → ∞ will arise if ρ rs µ (ω) L includes finite-volume energies for which E n (k) < E Σ (k) and similarly contributions growing as T b → ∞ will arise if σ rs µ (ω) L includes finite-volume energies for which E n (p) < E p (p). By studying the contributions in this specific system, we deduce that the limit T b → ∞ can be taken without any difficulties, since all possible baryonic intermediate states with strangeness S = −1 and momentum p have energies E n (p) > E p (p). This will be true so long as the strange quark mass is greater than the down quark mass. As a result, the term with e −T b (ω−Ep(p)) is exponentially suppressed for large T b . However, any intermediate state with S = 0 that has an energy smaller than the initial state energy E Σ (k) will lead to an exponentially growing term when T a → ∞. For simulations at physical or close-to-physical quark masses, such intermediate states can be either a single proton state with momentum k or a nucleon-pion state with total momentum k and an energy smaller than E Σ (k). For sufficiently large light-quark masses, the energy of the finite-volume nucleon-pion states will be greater than the Σ + energy and become decaying exponentials in T a , leaving only the single proton state to grow as T a → ∞.
To extract F rs µ (k, p) L from I rs µ (T a , T b ; p, k), all growing terms need to be removed. It is additionally possible to remove decaying states such that the T a → ∞ limit is saturated for smaller values of T a . The removal of slowly decaying states was already applied in ref. [7], in the context of rare kaon decays. To express this compactly, it is convenient to introduce a modification of ρ rs µ that is cut off to only include low-lying states where C n,µ (k) is defined in eq. (3.38) above, and N must satisfy the condition that E n (k) > E Σ (k) for n ≥ N . We stress here that n = 0 refers to the finite-volume single-proton state. Then one can write Note that I rs µ (T a , T b ; p, k) then has the desired large T a,b limits In contrast to I rs µ , the separate quantities I rs µ and ∆I rs µ , as well as F rs µ (k, p) L , have poles as a function of L for any fixed kinematics. The distinction arises because, in the original expression, 1 − e −(En(k)−E Σ (k))Ta vanishes whenever E n and E Σ coincide so that the combination has a finite limit This must be the case since the manifestly finite correlatorΓ (4)rs µ , integrated over a finite range of times, cannot diverge for any L. By contrast, ∆I rs µ , I rs µ and F rs µ are divergent if E Σ (k) coincides with a finite-volume energy. We stress that there is no problem here, this is simply part of the correct definition of the finite-volume estimator. These poles will be removed in the final relation between F rs µ and the infinite-volume Minkowski amplitude A rs µ . A slight variation in extracting F rs µ (k, p) L , discussed in ref. [10], is to remove the growing exponentials before integrating. In the present context it leads one to define This object now has a well-defined T a,b → ∞ limit, but it is not the desired estimator as the poles from the subtracted states are completely absent. These are then re-introduced by the relation where we have introduced .
This method is analogous to the approach of using low-lying states to estimate the T → ∞ integral for the hadronic-vacuum-polarization contribution to the magnetic moment of the muon [20]. Whether the removal of growing (and slowly decaying) exponentials is performed before or after t H integration, it requires determination of the overlaps C rs n,µ (k) and energies E n (k) of all states to be removed. We discuss the detailed approach for determining this information in the following subsections. Having completed this, it is equally important to understand how to relate F rs µ (k, p) L to the physical amplitude A rs µ (k, p). This requires treating the multi-hadron finite-volume effects and understanding how to include the N π branch cut that is part of the physical amplitude's definition. The method is discussed in detail in section 4. In addition, the general method for extracting the form factors from the physical amplitude A rs µ (k, p) is given in Appendix C, as well as an example for a specific kinematic setup.

Removal of the single-proton state
In this subsection, we describe two methods for removing the growing exponential arising from the single-proton state. Recalling the definitions denoting the single-proton state with a p subscript (i.e. C rs p,µ (k) = C rs n=0,µ (k)) one can show where in the last line we have used the hatted three-point functions defined in eqs. (3.29) and (3.30) and have also applied the identity P p (p) 2 = P p (p). This expression for C rs p,µ (k) can then be used to remove the single-particle state. One can define as the single-proton contribution to ∆I rs µ (T a ; p, k), defined in eq. (3.48). In fact, it is instructive here to consider the case where only the single proton leads to a growing exponential, as would be the case for sufficiently large pion mass calculations. Then the method for extracting F rs µ (p, k) L , with the proton removed after summation, can be summarized succinctly via The estimator with the proton removed before summation instead gives where one can directly make use of the time-dependence arising within eq. (3.30) One can readily show that these two expressions are mathematically equivalent. They may, however, lead to statistical differences depending on how exactly first term in eq. (3.60) and the first and last terms in eq. (3.61) are estimated.
Note also that it follows from eqs. (3.60) and (3.61) that, in the case where C rs p,µ (k) = 0, there is no need to explicitly treat the single-proton state. As was shown in refs. [6][7][8] for the single pion intermediate state in K → π + − , it is in fact possible to define a modified weak Hamiltonian, denoted H W (0), such that F rs µ (p, k) L is invariant under H W → H W , but p(k), r| H W (0) Σ + (k), s = 0 =⇒ C rs p,µ (k) = 0 , where the prime on C rs p,µ indicates the H W → H W replacement. The redefined Hamiltonian density is given explicitly by where c S and c P are constants to be determined and are flavour non-singlet scalar and pseudo-scalar densities. In the following we first prove that F rs µ (p, k) L is invariant under H W → H W and then explain how one fixes c S and c P to set the single-proton contribution to vanish.
Begin by recalling that the conserved and partially conserved vector and axial currents, Vq q µ and Aq q µ respectively, exactly satisfy the chiral Ward identities: Thus, inserting the scalar and pseudo-scalar densities between generic final and initial states with matching momenta, E f , k| L and |E i , k L respectively, one finds The crucial point is that Sd s (x) and Pd s (x) have no effect on F rs µ (p, k) L . To demonstrate this we define F rs µ (p, k) L as the result of replacing H W → H W in F rs µ , and then use the c S and c P dependence to unambiguously decompose as F rs µ (p, k) L = F rs µ (p, k) L − c S S rs µ (p, k) L − c P P rs µ (p, k) L , (3.69) thereby defining S rs µ (p, k) L and P rs µ (p, k) L . Taking the scalar for concreteness note that this can then be written explicitly as where the sum over n runs over finite-volume states with proton quantum numbers and that over n runs over states with strangeness S = −1. Using eq. (3.67) then gives where the energy differences in the chiral Ward identity have cancelled the poles. As a result each term now contains an insertion of the identity that can be collapsed to reach Similarly, for the pseudoscalar one finds Since the electromagnetic current J µ is a singlet in flavour space and the flavour changing axial and vector currents are not, they will commute causing the shifts to the estimator to vanish. We can therefore shift the weak Hamiltonian by any amount of the scalar and pseudoscalar operators without affecting the value of the estimator: It remains only to fix the values of c S and c P . This can be achieved by demanding where we have restored the ± superscript as it is relevant here that one can study the definite parity sectors separately. As shown in appendix A, a generic Lorentz (pseudo)scalar operator O can be decomposed into scalar and pseudoscalar form factors a O and b O respectively, giving the decompositions where we have used that only Sd s contributes to H + W (and only Pd s to H − W ). We deduce c S = a S a H and c P = b P b H are required for the single proton intermediate state to vanish. Comparing this to the scalar shift in the rare kaon decay [6], we note that due to the additional spin degree of freedom, there is an additional pseudoscalar form factor describing the weak Hamiltonian matrix element. In general, it is therefore not sufficient to perform only a scalar shift to remove the single proton intermediate state from both parity sectors. There is however a kinematic point where the pseudoscalar shift is no longer required. This is when the spinor contractionū r p (k)γ 5 u s Σ (k) vanishes at k = 0, corresponding to the Σ + at rest.

Removal of multi-hadron states
At close-to-physical quark masses, the other growing exponentials come from the lowestlying N π states, which have energy smaller than the mass of the Σ + . Because these are excited states, their subtraction is more involved. In general, all lattice interpolating operators with the correct quantum numbers will overlap the states of interest, but in practice one can only reliably extract excited state energies and matrix elements by solving a generalized eigenvalue problem (GEVP) with a diverse set of multi-hadron operators. In the present case, the low-lying states are expected to be N π like states together with resonances. Therefore, operators built from a nucleon and pion that are individually momentum projected are required to reliably determine the excited spectrum. Depending on the finite-volume box size, three-particle states could also become important, requiring an even more complicated operator basis.
It is also important to break up the intermediate excited states according to irreducible representations (irreps) of the octahedral symmetry group (including parity) or, in the case that the Σ + baryon has non-zero spatial momentum k, to a little group that leaves the latter invariant. The GEVP analysis is then performed separately within each irrep and the final result is constructed from the separately determined spectra. Details are given in Appendix D.
Once an optimized operator for a given excited state is determined, the removal of that state proceeds as in eq. (3.38) for the single proton. The energy E n (k) is extracted from the two-point function and the product of matrix elements from the three-point functions. It should be noted that this procedure can be applied to remove the exponentials associated with arbitrarily many intermediate states (whether growing or decaying), and for states with higher numbers of hadrons, so long as a sufficient operator basis can be obtained.
This completes our discussion of the construction of F rs µ (k, p) L . We now turn to the formalism required to relate this object to the physical amplitude for the weak decay Σ + → p + − .

Finite-volume effects
For calculations with sufficiently heavy pions, the energy of the lowest N π states will lie above the E Σ (k) threshold, and therefore only exponentially suppressed finite volume effects will be present. At the physical point however, the low-lying N π states will be below this threshold, inducing additional power-like finite volume effects.
In this section, we detail the correction of these power-like finite volume effects from such N π states which can be accounted for via a simple additive term, denoted by ∆F rs µ (k, p) L . This allows one to determine the physical amplitude, up to exponentially suppressed L effects, using the relation At this stage, we also switch from labelling states by definite individual quark flavour content (e.g. p and Σ + ) to isospin state labels (N for the neutron-proton doublet and N π for two-particle states with possible isospin values I = 1/2, 3/2). This is necessary to introduce the finite-volume formalism for multi-particle states below.
Since A rs µ (k, p) does not contain poles associated with the finite-volume N π intermediate states, these must cancel between the two terms on the right-hand side. As discussed in detail in ref. [10], this means that the numerical steps that introduce poles in F rs µ (k, p) L must match those in the construction of ∆F rs µ (k, p) L , such that the singularities exactly cancel. To keep track of this, it is useful to group the cancelling poles via the definitions An important subtlety here is that I rs µ (T a , T b ; p, k) and δF rs µ (T a ; k, p) L each diverge as T a → ∞, in such a way that the combination in eq. (4.2) remains finite.
An alternative to eq. (4.2) can be written by directly using the quantity I ≥N,rs µ (T a , T b ; p, k), defined in eq. (3.51). In particular, combining the various definitions given above, one finds where we have introduced a T a -independent analogue of δF rs µ (T a ; k, p) L , defined as This is the unique combination in this work that is (i) T a -independent, (ii) free of finitevolume singularities, and (iii) depends only on the states explicitly removed and not on the full sum over all states in the spectral decomposition. For this reason it will be useful in our detailed discussion of the single-channel case below.
We are now ready to give the full definition for ∆F rs µ (k, p) L . The finite-volume correction can be written as [6,10,21] ∆F rs where we have introduced .
Here F (E, P , L) is a known geometric function, reviewed in this work in the following paragraphs, and A r Jµ (k, p), A s H W (k), and M(E cm ) are three types of infinite-volume amplitudes, each involving N π states.
This construction and the following discussion is similar to that given in refs. [9,22], in that case of ππ intermediate states in K −K mixing and the decay K → πνν. The formalism of that work was generalized to particles with any spin, including the present case, in ref. [10]. The purpose of the following is to explain the relevant formalism in the specific context of N π states, including the role of parity, spin and non-degenerate masses. The formalism for long-range matrix elements also draws on that used for scattering and transition amplitudes [21,[23][24][25][26][27][28].
We begin by specifying the index space used in the definition of ∆F rs µ (k, p) L . The product on the right-hand side of eq. (4.6) is understood as a matrix (F) contracted with two vectors (A r Jµ and A s H W ) such that no hanging indices remain. The index space of this contraction is also required for the exact definitions of F (E, P , L) and the other quantities appearing above. The space is denoted by J, , µ, referring to the total angular momentum J, orbital angular momentum, , total spin (in this case fixed to 1/2), and azimuthal component of total angular momentum µ. For example, 11) and the combined labels of total energy E, total momentum k, and J, , µ are both necessary and sufficient to exactly specify the N π state. We define the geometric matrix, F (E, P , L), for the N π system, by first considering the quantity for two non-identical scalar particles with masses M π , M N [29]. Focusing on the P = 0 case, the definition reads We have also introduced p as the magnitude of back-to-back momentum, satisfying As discussed in Refs. [21,23], ultraviolet divergences in the sum-integral difference cancel so that any smooth regulator can be used in the evaluation of each. The explicit exponential included gives one option that is also convenient for numerical evaluation. As discussed in ref. [27], a straightforward combination of Clebsch-Gordon coefficients can then be used to promote this scalar version to the final F function on the full space for particles with spin: This is the quantity entering the definition of F in eq. (4.9). The second matrix entering the F matrix is the N π → N π scattering amplitude, denoted M. This can be represented on the J µ index space via where defines the centre-of-mass energy and δ J, (p) is the scattering phase shift with quantum numbers as indicated. This concludes our explanation of all quantities entering ∆F rs µ (k, p) L . The content of eq. (4.6) can be summarized as follows: By combining a determination of the N π → N π scattering amplitude with the Σ H W −→ N π and N π Jµ −→ N transition amplitudes, one can calculate the correction that relates the finite-volume estimator F rs µ (k, p) L to the amplitude A rs µ (k, p).

Determining M, A s H W and A r Jµ
So far we have made reference to poles cancelling between F rs µ (k, p) L and ∆F rs µ (k, p) L , but without explicitly explaining why the poles are the same between the two terms. This follows from the observation that the condition that F(E, P , L) diverges is equivalent to the Lüscher quantization condition [21,23] This relation allows one to constrain the amplitude M(E cm ) from a numerical lattice calculation by computing many finite-volume energies. The energies, E n , directly correspond to the poles in ∆F rs µ (k, p) L , as can be readily seen from its definition in eq. (4.6). In the limit of infinite statistics, the formalism guarantees a perfect cancellation between the poles in F rs µ (k, p) L and those in ∆F rs µ (k, p) L . However, for realistic data, with statistical uncertainties, care must be taken. In particular, given the best fit for the lattice energies E n (L) from Euclidean correlators, one can use a given parametrization in the quantization condition to determine the best fit for M(E cm ). If the lattice energies (rather than the best-fit quantization energies) are used to construct F rs µ (k, p) L singularities will fail to cancel in the amplitude. As discussed in ref. [10], the most straightforward solution here is to directly use the best-fit quantization energies everywhere.
The Lellouch-Lüscher formalism [25] and its extensions can be used to determine the remaining amplitudes A s H W and A r Jµ . This follows from the eigenvectors of the quantization condition matrix which is known to be rank one near the finite-volume energy where the ⊗ indicates an outer-product in the same index space, so that E (n),in and E (n),out each carry the J, , µ indices. Then the amplitudes are related to finite-volume matrix elements via up to an inherent phase ambiguity that will cancel between E (n),out and E (n),in in the full finite-volume correction.

Single-channel case
We now show that, in the limit where only a single channel contributes, that this formalism is equivalent to that of ref. [9]. We first consider the general case in which L is not tuned such that E n (L) = M Σ and then examine the finely tuned case discussed in ref. [9]. We also extend the latter by expanding to higher orders about the tuned volume point. This allows one to propagate uncertainties, given that E n (L) = M Σ will never hold exactly in estimators of the variance. For the purpose of this discussion, we focus attention on the case of zero momentum in the finite-volume frame (P = 0). We specifically consider the case of a single N π channel in the P -wave ( = 1). In this case the matrices M(E) and F (E, L) are one-dimensional. It is standard to describe these in terms of the scattering phase shift δ =1 (E) and the so-called pseudophase φ(E, L): Substituting this into eq. (4. 19) gives This can be used to express the relation between finite-and infinite-volume matrix elements as where p(E) is given by eq. (4.13). Then, in general, one can envision parametrizing the unknown function A r Jµ (E, 0, p)A s H W (E, 0) as a function of E and constraining the unknown parameters using eq. (4.26).
Taking this function as known allows us to construct δF rs µ (k, p) L as defined in eq. (4.5). This is a useful quantity to focus on as it only depends on the states that have been removed explicitly, together with the finite-volume correction.
where p Σ = p(M Σ ). This can alternatively be written as where, with the exception of the single-proton state, we have expressed also the finite-volume matrix elements within C n,µ (0) in terms of the infinite-volume transition amplitudes. Here we have also introduced the shorthand ∂ E = ∂/∂E to express the derivative factor in the first term more compactly. The result summarized in eq. (4.28) matches eq. (35) of ref. [9] up to some superficial differences in notation, motivated by the fact that we are describing a different system. For example in contrast to the case of the K L -K S mass difference described in ref. [9], here the two infinite-volume amplitudes entering the correction, A r Jµ (E n , 0, p) and A s H W (E n , 0), are different. For this reason we can not express the result as a magnitudesquared but instead explicitly include the Watson phase e −2iδ such that the combination A r Jµ e −2iδ A s H W is real-valued (but not necessarily positive).
This concludes our discussion of the general finite-volume formalism. Before moving to the next section, which describes the case of a volume tuned such that M Σ coincides with a finite-volume energy, we close here by summarizing the necessary steps to obtain the full amplitude A rs µ (k, p) with all power-like finite-volume corrections removed. First, one needs to compute the finite-volume estimator P B (p) F µ (k, p) L P B (k) on the lattice, giving access to its spin-projected counterpart, F rs µ (k, p) L , as described in eqs. (3.49) and (3.52). The next step is to perform an N π scattering analysis as has been done in refs. [30,31] and as we outline in appendix D. This involves obtaining the finite-volume energies, which via the generalized Lüscher formalism lead to parametrizations of the phase shift δ(E) as well as its derivative. Similarly, a lattice computation of the three-point functions on the right-hand side of eq. (4.26), together with the generalized Lellouch-Lüscher relation, gives access to the infinite-volume transition amplitudes A r Jµ (E n , 0, p) and A s H W (E n , 0). Putting everything together, these ingredients can now be put into eq. (4.28) to obtain δF rs µ (0, p) L and the final amplitude via eqs. (4.2), (4.3), and (4.5). We stress that the final +i contribution in eq. (4.28) allows direct access to the imaginary part of the amplitude A rs µ (k, p).

Expanding about the pole
The result of the previous section gives the finite-volume correction for any value of L for the case of E = M Σ . In this section we denote byL andn particular values such that En(L) = M Σ for one of the finite-volume states. The assumption that one is exactly tuned to such a volume is built into ref. [6] and performing the expansion here allows us to make contact to the earlier work.
Cancellation of leading-order singularity Defining δL = L −L we first note that the expression for δF rs µ (0, p) L has a pole that must cancel. The leading order expansion explicitly gives

(4.29)
To see that the term in square brackets vanishes, first use the identity where we have introduced the shorthand ∂ L = ∂/(∂L). In the first line it is understood that the total derivative is being evaluated along the solution. This allows us to rewrite the L derivative of the energy to reach Finally, using and replacing φ(M Σ ,L) with −δ =1 (M Σ ), we see the desired cancellation to deduce that δF µ remains finite as δL → 0.
Result for L =L To derive the (δL) 0 term simply requires carefully expanding to the next order. Performing the expansion about δL = 0 leads to numerous terms, depending on both first and second derivatives of En(L) with respect to L. The result can be written as where we emphasize that n = 0 (corresponding to the single-proton state) is included in the sum in the first term. We have reached this result both by expanding about δL (equivalently taking the L →L limit) with E = M Σ fixed and, alternatively, by setting L =L exactly and taking the E → M Σ limit. As expected the same result is recovered in both approaches.

O(δL) correction
Finally, one can perform the tedious exercise to push this expansion to one higher order in δL. The result can be written as where the square brackets and superscript on the left-hand side indicate that we are reporting only the linear order, to be added to the constant order above. Here we have introduced the lengthy expressions together with the shorthand A Mathematica version of this can be provided upon request to the authors. Although this result looks quite cumbersome, we argue that it is in fact both conceptually straightforward and quite useful in practice. Given the value of the phase shift and its first and second derivatives at E = M Σ and L =L together with the amplitudes A r Jµ (M Σ , 0, p) and A s H W (M Σ , 0), all information is available to express the above result as a single numerical coefficient multiplying (M Σ δL). If one is in the regime where the volume is sufficiently well tuned so that δL 2 terms can be neglected, then this gives the relevant correction from the finite-volume formalism due to inevitable slight mistuning. In particular, in a jackknife or bootstrap error estimate, the above expression gives a straightforward way to evaluate δF rs µ (0, p) L on each sample and thus to propagate uncertainties into the physical amplitude.

Summary
In this paper, we described the different steps necessary to extract the Σ + → p + − rare decay amplitude from Euclidean lattice correlation functions. The main challenge is the radically different contributions from on-shell intermediate states between Minkowski and Euclidean finite-volume correlation functions. In Euclidean space-time, these states create spurious contributions in the time-integrated correlation function that grow exponentially with the time range. In this paper we explain in detail how these contributions can be reconstructed and subtracted, which generalizes for the hyperon decays what was previously discussed in the kaon sector [6,22], and more generally in [10]. Even after subtracting growing exponentials, the finite-volume estimator for the amplitude still suffers from finitesize effects following a power-law in the volume extent L. We presented the asymptotic volume behaviour of the finite-volume estimator in different kinematics, and compared our result to previous approaches. We also expand on previous work by providing the leading order correction of the finite-volume effects for lattice volumes tuned to have a finite-volume energy level near the Σ mass.
This work paves the way for predicting this amplitude in future lattice calculations, and we are in parallel conducting exploratory numerical simulations to that effect [32]. Beyond the theoretical problems presented in this paper, one can expect important numerical challenges to be addressed, in particular related to the generally poor statistical signal of baryonic correlation functions in lattice calculation. However, such lattice calculation will generate numerous interesting results beyond the exciting perspective of studying this flavour-changing neutral current decay observed at LHCb. Indeed, tackling the growing exponential issue will generate as a side-product first-principle determinations of N π scattering parameters. More generally, this type of process is a perfect ground to extend our general understanding on how to predict non-local hadronic matrix elements from lattice simulations.

A Form-factor decomposition (Minkowski)
Generally the matrix element of an operator O between two baryonic states |B(p), s and |B (p ), r , with momenta p and p and spin projections s and r can be split into the external spinors u and a 4 × 4 Dirac-matrix amplitude A where we now parametrize in terms of the 4-momenta q = p − p and k = p + p . We then write A in terms of all of the available Dirac-matrix structures Γ ∈ {1, γ 5 , γ µ , γ µ γ 5 , σ µν } in such a way as the Lorentz structure of the matrix element is recovered. These terms are associated with a form factor which is a scalar coefficient, and can only be a function of the Lorentz and spin scalar objects q 2 and k 2 . However, k 2 = 2m 2 + 2m 2 − q 2 and therefore k 2 and q 2 are not independent, so the forms factors can be written as only a function of q 2 . The base set of Dirac-matrix Lorentz scalar objects is: The full list of structures is the infinite set of all combinations of these base elements. These additional structures can always be decomposed into a linear combination of the base elements using / p / p = p 2 and / p / p = 2p · p − / p / p for any p, p . The base set of spin-matrix Lorentz vectors objects is: {q µ , q µ γ 5 , k µ , k µ γ 5 , γ µ , γ µ γ 5 , σ µν q ν , σ µν q ν γ 5 , σ µν k ν , σ µν k ν γ 5 } , where again the infinite combinations can be decomposed in terms of the base set.

A.1 Explicit form-factor decompositions
Using the base set for a generic Lorentz (pseudo-)scalar operator S as given above yields the form factor decomposition which is relevant for the matrix element of the weak Hamiltonian in the rare hyperon decay. Analogously, using the base set for a generic Lorentz (axial-)vector operator J µ one obtains +g 1 (q 2 )γ µ γ 5 + ig 2 (q 2 )σ νµ q ν γ 5 + g 3 (q 2 )q µ γ 5 u s B (p) , where we have used the generalisations of the Gordon decomposition identity with initial and final states of different mass Finally, the relevant amplitude for the rare hyperon decay Σ + → p + − decomposes like where we have used the Ward-Takahashi identity q µ A rs µ = 0. Note the different definitions of momenta in this appendix and in the main text, where the p carries momentum p and the Σ + carries momentum k.

B Minkowski and Euclidean definitions and conventions
In the following we give details on our conventions in Minkowski and Euclidean metric. Quantities with a sub-or super-script "M " and "E" denote quantities in Minkowski and Euclidean metric respectively 4 .

B.1 Kinematic variables
For the Minkowski metric we use the mostly-minus convention, i.e. g µν = g µν = diag(1, −1, −1, −1). Then the four-momentum of a state with energy E and spatial momentum p is given by leading to Of course, the limitation of the Euclidean signature in numerical lattice calculations is that one can often only access Euclidean momenta with real components, so that some effort is needed to access extract observables for which p 2 E < 0.

B.2 Gamma matrices
We define the Euclidean γ-matrices from their Minkowski-counterparts as follows: In total, for the rare-hyperon amplitude we obtain As derived above, written with Minkowski quantities the form factor decomposition for the rare hyperon decay is given as For µ = 0 one finds that this can be written with Euclidean quantities as follows: Similarly, for µ = j can be written as In equations (B.17) and (B.18) we have defined with the rare-hyperon decay matrix element written with Euclidean quantities.

C (Euclidean) traces for extracting the form factors
All quantities (Γ-matrices etc) in this section are understood to be the Euclidean versions, and we drop sub-/superscripts "E" for better readability. Details of our conventions for Euclidean metric can be found in Appendix B above. In the following we give examples of traces and combinations that can be used to extract the form factors a(q 2 ), b(q 2 ), c(q 2 ) and d(q 2 ). For convenience we define tr Γ µ ≡ Tr Γ(−i / p p + M p )Ã µ (−i/ k Σ + + M Σ + ) (C.1) withÃ µ = σ νµ q ν [a − γ 5 b] + −q 2 γ µ + q µ/ q [c − γ 5 d] . (C.2) For example, one finds ] .

(C.4)
The goal is to find suitable combinations of the traces tr Γ µ to extract the four form factors a, b, c and d. As an example we will give results for the kinematics where the Σ + is at rest k = 0 and the proton is moving along the x-direction p = (p p , 0, 0). Table 1: Traces for µ = 0 and k = 0, p = (p p , 0, 0). All traces that are not explicitly given are 0. Table 2: Traces for µ = 2 and k = 0, p = (p p , 0, 0). All traces that are not explicitly given are 0.

D Details on extracting finite-volume excited states
We assume throughout that the Σ + is at rest in the finite-volume frame. Then the infinitevolume quantum numbers, J P = 1/2 + , of the |Σ + (0) state, imply that it transforms in the G + 1 irrep [33,34] of the double-cover of the octahedral group (with parity). In addition to the J = 1/2 channel, this irrep couples to J = 7/2, which corresponds to orbital angular momentum values of = 3, 4. Neglecting this mixing, by formally assuming that the N π → N π scattering amplitude vanishes for ≥ 3, the only sectors we have to consider are the = 0 (P = −) and = 1 (P = +) states with total angular momentum J = 1/2.
As discussed in the main text, the weak Hamiltonian density is a sum of a scalar (H + W ) and pseudo-scalar (H − W ) operators. The corresponding states H + W |Σ + (0) and H − W |Σ + (0) transform in the G + 1 and G − 1 irreps respectively and therefore couple to multi-hadron excited states with these quantum numbers. The latter correspond to different infinite-volume 5 Using other combinations is possible as well; here we give only one of many possibilities. With this information we can now construct the relevant states to remove from the amplitude. In both these irreps, the calculation is analogous to a phase-shift analysis for J = 1/2 N π scattering, similar to what has been done in [30,31]. The operator basis is different in the two irreps and includes both N π interpolators with different back-toback momenta (for G − 1 , this includes an operator with both operators projected to zero spatial momentum, but for G + 1 , at least one unit of back-to-back momentum is required). Additionally, the Roper N (1440) resonance is present in the G + 1 irrep and a local quark bilinear corresponding to this might be included. In a more ambitious calculation threehadron N ππ operators might also be incorporated.
Once the basis is established one can apply a suitable technique (e.g. the method of distillation [35,36]) to form all possible two-point correlation functions from the operators, allowing one to construct a correlator matrix where Λ ∈ {G − 1 , G + 1 } denotes the irrep. In a practical implementation, this basis should include more operators than the number of finite-volume states one plans to reliably extract. The next step is to apply the variational method, which requires one to solve a GEVP of the form [37][38][39] C Λ (t)v(t, t 0 ) = λ(t, t 0 )C Λ (t 0 )v(t, t 0 ) .

(D.4)
From the eigenvalues one can then extract the energy spectrum of the system via λ n (t, t 0 ) = e −E Λ n (t−t 0 ) + Be −∆E Λ n (t−t 0 ) + · · · , (D.5) where, for a suitable choice of t 0 , ∆E Λ n is an energy gap to the lowest level not covered by the operator basis [40]. With the eigenvectors arising from the GEVP, an optimised interpolator can then be formed which couples particularly well to the n th state in the system. These can then be used exactly as the proton creation operators in the main text, in order to extract the relevant energies and matrix elements both for determining F µ (p, k) L and the finite-volume corrections that define the amplitude.