P-Wave Two-Particle Bound and Scattering States in a Finite Volume including QED

The mass shifts for two-fermion bound and scattering P-wave states subject to the long-range interactions due to QED in the non-relativistic regime are derived. Introducing a short range force coupling the spinless fermions to one unit of angular momentum in the framework of pionless EFT, we first calculate both perturbatively and non-perturbatively the Coulomb corrections to fermion-fermion scattering in the continuum and infinite volume context. Motivated by the research on particle-antiparticle bound states, we extend the results to fermions of identical mass and opposite charge. Second, we transpose the system onto a cubic lattice with periodic boundary conditions and we calculate the finite volume corrections to the energy of the lowest bound and unbound $T_1^{\pm}$ eigenstates. In particular, power law corrections proportional to the fine structure constant and resembling the recent results for S-wave states are found. Higher order contributions in $\alpha$ are neglected, since the gapped nature of the momentum operator in the lattice environnement allows for a perturbative treatment of the QED interactions.


Preamble
Effective Field Theories [1,2,3,4,5,6,7,8] nowadays play a fundamental role in the description of many-body systems in nuclear and subnuclear physics, employing the quantum fields which can be excited in a given regime of energy. Once the breakdown scale Λ of the EFT is set, the scattering amplitudes are usually expressed in power series of p/Λ, where p represents the characteristic momentum of the processes under consideration. The Lagrangian density is typically written in terms of local operators of increasing dimensions obeying pertinent symmetry constraints.
Moreover, power counting rules establish a hierarchy among the interaction terms to include in the Lagrangian, thus permitting to filter out the contributions that become relevant only at higher energy scales [7].
In the case of systems of stable baryons at energies lower than the pion mass, the Lagrangian density contains only the nucleon fields and their Hermitian conjugates, often combined toghether with differential operators. The corresponding theory, the so-called pionless EFT [2,9,10,11,12] counts a number of successes in the description of nucleon-nucleon scattering and structure properties of few-nucleon systems. Despite the original difficulties in the reproduction of S-wave scattering lengths, that were solved via the introduction of the Power Divergence Subtraction (PDS) as a regularization scheme [10,13,14], the theory has permitted so far to reproduce the 1 S 0 np phase shift [15,16], structure properties of the triton as a dn S-wave compound [12,17,18] and the scattering length [19,20] and the phase shift [21,22,23] of the elastic dn scattering process.
In the first applications of QED in pionless EFT, the electromagnetic interactions were treated perturbatively, as in the case of the electromagnetic form factor [24] and electromagnetic polarizability [25] for the deuteron or the inelastic process of radiative neutron capture on protons [26].
Afterwards, a non-perturbative treatment of electromagnetic (Coulomb) interactions on top of the same EFT was set up, in the context of proton-proton S-wave elastic [27] and inelastic [28,29] scattering.
Inspired by the P-wave interactions presented in refs. [30,31,32], we generalize in the first part of the present paper the analysis in ref. [27] to fermion-fermion low-energy elastic scattering ruled by the interplay between the Coulomb and the strong forces transforming as the = 1 representation of the rotation group (cf. ref. [33] for the empirical S-and P-wave phase shifts in the pp case). As in ref. [27], we treat the Coulomb photon exchanges both in a perturbative and in a non-perturbative fashion. During the derivation of the T-matrix elements, we observe that at sufficiently low energy the repulsion effects from the Coulomb ladders become comparable to the ones of the strong forces, leading to the breakdown of the perturbative regime of non-relativistic QED. In the determination of the closed expressions for the scattering parameters in terms of the coupling constants, we take advantage of the separation of the Coulomb interaction from the strong forces, considered first in refs. [34,35,36] and eventually generalized to strong couplings of arbitrary angular momentum in ref. [37]. The importance of particle-antiparticle systems led us to the applicaton of the formalism to fermion-antifermion scattering, where the attractive Coulomb force gives rise to bound states. This case provides a laboratory for the study of pp bound [38] and unbound states [39], also referenced as protonium.
Of fundamental importance for the study of few-and many-particle systems with QED are Lattice Effective Field Theories and Lattice Quantum Chromodynamics (LQCD). The latter has matured to the point where basic properties of light mesons and baryons are being calculated at or close to the physical pion mass [40]. In particular, in the case of the lowest-lying mesons, their properties are attaining a level of accuracy where it is necessary to embed the strong interactions within the full standard model [41,42,43,44,45,46]. Despite the open computational challenges represented by the inclusion of the full QED in LQCD simulations, in the last decade quenched QED [47] together with flavour-symmetry violating terms have been included in the Lagrangian, with the aim of reproducing some features of the observed hadron spectrum [48,49,50,51,52].
Conversely, the perspective to add QED interactions in LQCD simulations for light nuclei appears still futuristic, due to the limitations in the computational resources. Nevertheless, in two-body scattering processes like π ± π ± scattering [53,54], the time is ripe for the introduction of QED interactions in the present LQCD calculations.
It is exactly in this context that, in the second part of the paper, we immerse our fermionfermion EFT into a cubic box with periodic boundary conditions (PBC). The lattice environment has a number of consequences, the most glaring of them are the breaking of rotational symmetry [55,56,57,58] and the discretization of the spectrum of the operators representing physical observables [59,60,61]. Concerning the Hamiltonian, its spectrum consists of levels that in the infinite-volume limit become part of the continuum (scattering states) and in others that are continuously transformed into the bound states. For two-and three-body systems governed by strong interactions, the shifts of the bound energy levels with respect to the counterparts at infinite volume depend on the lattice size L through negative exponentials, often multiplied by nontrivial polynomials in L. Apart the pioneering work on two-bosons subject to hard-sphere potentials in ref. [62], these effects for two-body systems have been extensively analyzed by Lüscher in refs. [63,64] ( [65]), where the energy of the lowest unbound (bound) states has been expressed in terms of the scattering parameters and the lattice size.
However, the presence of the long-range interactions induced by QED leads to significant modifications in the form of the corrections associated to the finite volume energy levels. Irrespective on whether a state is bound or unbound, in fact, the energy shifts take the form of polynomials in the reciprocal of the lattice size [40] and the exponential damping factors disappear. Moreover, the gapped nature of the momentum of the particles on the lattice allows for a perturbative treatment of the QED contributions, even at low energies [40,47,87,88]. In this regime, composite particles receive corrections of the same kind both in their mass [40] and in the energies of the two-body states that they can form [88].
As shown in ref. [88], the leading-order energy shift for the lowest S-wave bound state is proportional to the fine-structure constant and has the same sign of the counterpart in absence of QED, presented in refs. [60,66]. In the second part of this work, we demonstrate that the same relation holds for the lowest bound P-wave state, whose finite volume correction is negative as the one for the counterpart without electromagnetic interactions. Additionally, we prove that the QED energy-shifts for S-and P-wave eigenstates have the same magnitude if order 1/L 3 terms are neglected, a fact that remains valid in the absence of interactions of electromagnetic nature. At least for the = 0 and 1 two-body bound eigenstates, in fact, the sign of the correction depends directly on the parity of the wavefunction associated to the energy state, whose tails are truncated at the boundaries of the cubic box, as observed in ref. [60].
Concerning scattering states, the energy shift formula for the lowest P-wave state that we present in this paper has close similarities with the one in ref. [88], despite an overall ξ/M ≡ 4π 2 /M L 2 factor, owing to the fact that the energy of the lowest unbound state with analogous transformation properties under discrete rotations (T 1 irrep 1 of the cubic group) is different from zero. Additionally, further scattering parameters appear in the expression for the = 1 finite volume energy correction, even as coefficients of the smallest powers of 1/L.
The present article is structured into two parts and its content can be summarized as follows.
After this preamble, the theoretical framework that is the basis for both the infinite and the finite volume treatment is introduced, by starting from the Lagrangian with the strong P-wave interactions alone. Next, in the end of sec. 2.0, the T-matrix for two-body fermion-fermion scattering to all orders in the strength parameter of the potential is computed. Subsequently, in sec. 2.1, non-relativistic QED is presented in the same fashion of refs. [89,90] and the Lagrangian is reduced to the case of spinless fermions and electrostatic interactions. After displaying the amplitudes corresponding to tree-level and one-loop diagrams with one Coulomb photon exchange, the non-perturbative treatment of the Coulomb interaction is implemented. To this aim, we resort to the formalism of ref. [27], that we recapitulate in the end of sec. 2.1. As in the introductory section, the T-matrix matrix element accounting for both Coulomb and strong interactions is derived to all orders in α, thanks to the Dyson-like identities that hold among the free, the Coulomb and the full two-body Green's functions. As in ref. [27], section 2.2 closes with the expressions of the scattering length and the effective range in terms of the physical constants of our EFT Lagrangian, that are obtained from the effective-range expansion. The first part of the analysis is concluded in sec. 2.3 with the calculation of the same amplitude for the fermionantifermion scattering case. 1 Throughout, we use the abbreviation "irrep" for an irreducible representation.
Then, the two fermion-system is transposed onto a cubic lattice of size L and the distortions induced by the new environment in the laws of electrodynamics [40,91] and in the masses of possibly composite particles [40] are briefly summarized in sec. 3.0. Next, the quantization conditions, that give access to the energy spectrum in finite volume though the expression of the T-matrix elements [88], are displayed and discussed (sec. 3

Effective Field Theory for non-relativistic fermions
Our analysis of two-particle scattering and bound states in the infinite-and finite-volume context is based on pionless Effective Field Theory [9,10,13,14,92,93,94,95]. The theory, developed more than two decades ago [9], describes the strong interactions between nucleons at energy scales smaller than the pion mass, M π [7,27,92]. The action is non-relativistic and is constructed by including all the possible potential terms made of nucleon fields and their derivatives, fulfilling the symmetry requirements of the strong interactions at low energies (i.e. parity, time reversal and Galilean invariance) [96]. The importance of the various interaction terms decreases with their canonical dimension while approaching the zero energy limit. Besides, even the dominant contribution at low energies for local contact interactions between four-nucleon fields is of dimension six, thus making the theory non-renormalizable [27] in the classical sense.
Analogously to ref. [88], we begin by extending pionless EFT to spinless fermions of mass M and charge e, and we assume that the theory is valid below an upper energy Λ E * in the center-ofmass frame (CoM). More specifically, if the fermions represent hadrons, the latter energy cutoff can be chosen to coincide with the pion mass. Second, we construct the interactions in terms of four-fermion operators, selecting the ones that transform explicilty as the 2 + 1-dimensional where P is a Legendre polynomial, ±p (±q) are the three-momenta of two incoming (outcoming) particles in the CoM frame, such that |p| = |q|,V ( ) is the potential in terms of second quantized operators and the c ( ) 2j are low-energy (LECs) constants, whose importance at low-energy scales diminishes for increasing values of j . In particular, for the three lowest angular momentum couplings ( ≤ 2), the interaction potentials take the form and As shown in sec. II of ref. [88], the terms in eq. (2) (eq. (3) and (4)) proportional to even powers of the momentum (a gradient expansion in configuration space), can be encoded by a single interaction with energy-dependent coefficient C(E * ) (D(E * ) and F (E * )) for S-waves (P-and Dwaves), where E * represents the CoM energy of the colliding particles, equal to 2M + p 2 /M .
While the case of fermions coupled to zero angular momentum via a single contact interaction proportional to C(E * ) is the starting-point of the analysis in ref. [88], the fundamental P-wave interaction in eq. (3) with energy-dependent coefficient D(E * ) becomes the key tool of the present investigation. Although interactions of the same form have been already adopted in pionless EFT for nucleons (cf. eq. (4) in ref. [30]) and in EFT with dimeron fields (cf. eq. (2) in ref. [97]), the P-wave counterpart of Kong and Ravndal's analysis on fermion-fermion scattering in ref. [27] is not available in the literature. Adopting the conventions of ref. [30] for the coupling constants (cf. the Feynman rules in app. A), the Lagrangian density assumes the form where  Fig. 3 in ref. [27]. In particular, the tree-level diagram, consisting of a single four-fermion vertex, leads to an amplitude equal to −iD(E * )p · p (cf. ref. [32]) where ±p and ±p are, respectively, the momenta of the incoming and outcoming particles in the CoM frame. As a consequence, the two-body = 1 (pseudo)potential in momentum space takes the form which coincides with the tree-level diagram multiplied by the imaginary unit. Considering also the other possible diagrams in momentum space with amputated legs in fig. 1, the expression for the full scattering amplitude due to strong interactions can be written as, whereĜ E 0 ≡Ĝ (+) 0 (E) is the two-body unperturbed retarded (+) Green's function operator, withĤ 0 the two-body free Hamiltonian in relative coordinatesĤ 0 =p 2 /M , and M/2 is the reduced mass of a system of identical fermions. Inserting a complete set of plane wave eigenstates |q, −q ≡ |q in the numerator, the latter expression becomeŝ where E = p 2 /M is the energy eigenvalue at which the retarded (+) and advanced (−) Green's functions are evaluated. In configuration space the latter take the form that is diagrammatically depicted by two propagation lines. The explicit computation of the three lowest order contributions to the sum in eq. (7) yields and i q, −q|V (1)ĜE

0V
(1)ĜE 0V where ∂ i ≡ ∂/∂r i , ∂ i ≡ ∂/∂r i and ∂ i ≡ ∂/∂r i , while J 0 is a symmetric matrix whose elements are given by and Einstein's index convention is henceforth understood. Extending the computation to higher orders, it is evident that the infinite superposition of chains of bubbles translates into a geometric series in the total scattering amplitude, as in the = 0 case, and a formula analogous to eq. (2) in ref. [27] is obtained, Furthermore, performing the Fourier transform of the potential in eq. (6) into configuration space, the full scattering amplitude can be recovered independently in position space by means of partial integrations and cancellations of surface integrals at infinity, The matrix elements of J 0 can be, now, evaluated by dimensional regularization. Applying the formula in eq. (B18) of ref. [98] for d-dimensional integration, eq. (14) in arbitrary d-dimensions becomes where µ is the renormalization scale introduced by the minimal subtraction (MS) scheme. Like the S-wave counterpart, the integral proves to be finite in three dimensions and, within this limit is given by where the energy E in the CoM frame has been eventually expressed as p 2 /M . For the sake of completeness, we derive the contribution to (J 0 ) ij from the power divergence subtraction (PDS) regularization scheme, in which the power counting of the EFT is manifest [10,14]. To this aim, the eventual poles of the regularized integral for d → 2 should be taken into account. In this limit, it turns out from eq. (18) that the Euler's Gamma has a pole singularity of the kind 2/(2 − d).
As a consequence, the original dimensional regularization result in eq. (18) acquires a finite PDS contribution, transforming into This can be compared with the one in eq. (4) in ref. [27] for the S-wave interactions. Since the J 0 matrix is diagonal (eq. (20)), few efforts are needed for the computation of the fermion-fermion scattering amplitude, With reference to scattering theory [99], the T S matrix for P-wave elastic scattering with phase shift δ 1 can be written as where θ is the angle between the incoming and outcoming direction of particles in the CoM frame.
Recalling the effective-range expansion (ERE) for = 1 scattering [99], an expression for the scattering parameters in terms of the momenta of the particles, the coupling constant and the mass M can be drawn. In particular, a formula for the scattering length analogous to eq. (2.16) of ref. [14] can be recovered, Furthermore, the effective range parameter r 0 vanishes, as in the zero angular momentum case.
Plugging the PDS-regularized expression of J 0 in eq. (20) into eq. (15) and exploiting the ERE again, finally, the renormalized form of the coupling constant D(E * ) is obtained, Unlike in the = 0 case, we note that the µ-dependent version of D(E * ) = 12πa/M is quadratic in the momentum of the incoming fermions.

Coulomb corrections
We introduce the interactions of electromagnetic nature in the non Lorentz-covariant fashion of ref. [90] and [89]. The formalism of non-relativistic quantum electrodynamics (NRQED), introduced in ref. [89], is designed to reproduce the low-momentum behaviour of QED to any desired accuracy. Besides, only non-relativistic momenta are allowed in the loops and in the external legs of the diagrams. The contributions arising from relativistic momenta in the QED loops, in fact, are absorbed as renormalizations of the coupling constants of the local interactions in the nonrelativistic counterpart of QED [89]. The Lagrangian is determined by the particle content and by the symmetries of the theory, such as gauge invariance, locality, hermiticity, parity conservation, time reversal symmetry and Galilean invariance. The particles are fermionic, characterized by mass M and unit charge e, and are represented by two-component non-relativistic Pauli spinor fields Ψ. In compliance to these prescriptions, the NRQED Lagrangian density in ref. [89] assumes the form, where D = ∇ + ieA is the covariant derivative, while E = −∇φ − ∂ t A and B = ∇ × A denote the electric and magnetic fields, respectively. The terms in the first row encode the leading ones of L NRQED , containing the minimal coupling of the fermionic fields with the vector potential A and the scalar potential, φ. The interactions proportional to the constants c 1 -c 4 and d 1 in eq. (26) are next-to-leading-order terms, corresponding to corrections of order v 4 /c 4 and v 6 /c 6 , respectively [90], whereas the ellipses represent contributions containing higher order covariant derivatives, Since the Coulomb force dominates at very low energies and transverse photons couple proportionally to the fermion momenta, in the present treatment we choose to retain in the Lagrangian only the scalar field and its lowest order coupling to the fermionic fields as in ref. [27]. Moreover, we reduce the latter to spinless fields ψ, consistently with sec. 2 and with ref. [88]. As a consequence, the full Lagrangian density of the system becomes the superposition of the one in eq. (5) with the one involving the electrostatic potential and its leading order coupling to the spinless fermions, namely Alternatively, on top of the P-wave interaction in eq. (6) the Coulomb force, that in momentum space regulated by an IR cutoff λ, reads  From the Feynman rules, the amplitude for the tree-level diagram with one photon insertion in the left part of fig. 2 gives, After integrating over the free energy l 0 , the tree-level amplitude in eq. (29) can be decomposed as follows In particular, in the last rewriting the first integral on the r.h.s. turns out to be identical to the one in eq. (9) of ref. [27] except for the factor p · p = p 2 cos θ, therefore it can be immediately integrated. Conversely, the last integral in eq. (30) represents a new contribution, whose evaluation in dimensional regularization is carried out in app. B. Adding the two contributions together, the tree level amplitude with Coulomb photon insertion in eq. (30) becomes where the limit λ → 0 for the O(λ) terms is understood. From the last equation we infer that, due to the linear dependence in the momenta of the incoming particles p in the CoM frame, Pwave fermion-fermion scattering is suppressed with respect to the S-wave one in the low-p limit.
However, both the = 0 and = 1 tree level amplitudes with one Coulomb photon insertion are divergent in the λ → 0 limit.
Nevertheless, due to the fact that the latter logarithmic contribution is imaginary, the infinite term does not contribute to the O(α) corrections of the strong cross section, which is proportional up to first order in α. The corrected cross section to that order turns out to be IR finite and proportional to 1 − πη, where η = αM/2|p| for particles with equal unit charge. As observed in ref. [27], the inclusion of n Coulomb photon exchanges leads to corrections proportional to η n in the cross section. Therefore, the feasibility of a perturbative treatment for the Coulomb force is regulated by the smallness of the parameter η, i.e. by a constraint on the momenta of the incoming particles, |p| αM/2. As a consequence, if the momenta of the incoming particles are too small, the Coulomb force is expected to have a strong influence on the cross-section of the elastic process and a non-perturbative treatment becomes necessary.
Furthermore, the Feynman rules for the one-loop diagram with one photon insertion on the right part of fig. 2 yield Similarly, the contour integration with respect to the free energies k 0 and l 0 , followed by the momentum translation l → q ≡ l − k, leads to The remaining momentum integrations are performed in dimensional regularization (cf. sec. B) and give where ≡ d − 3 and γ E ≈ 0.5772 denotes the Euler-Mascheroni constant. The amplitude in eq. (34) displays a pole at d = 3 as the one in eq. (15) of ref. [27], an ultraviolet divergence that can be reabsorbed with a redefinition of the strength parameter D(E * ) via the renormalization process. However, T 1−loop SC is devoid of the logarithmic divergence in the zero-momentum limit, due to the multiplication by a factor p 2 . Consequently, in comparison with the zero angular momentum counterpart, the = 1 one-loop scattering amplitude with one-photon insertion is suppressed in the limit of zero momentum p of the incoming particles in the CoM frame. Since possesses also a pole in the d → 2 limit, the implementation of the PDS scheme results into an additional term proportional to the renormalization scale (or mass), Differently from the = 0 counterpart in eq. (14) of ref. [27], the logaritmic term in the CoM momentum of the colliding fermions does not give rise to a divergence in the zero momentum limit, due to the p 4 prefactor. Nevertheless, the dressing of the one-bubble diagram with two or more Coulomb photon insertions results in the multiplication of T 1−loop SC by one or more powers of η = αM/2|p| 2 , so that, at order higher than four in α, the amplitude becomes singular in the limit |p| → 0. It follows that the perturbative approach breaks down and the effects of Coulomb repulsion need to be treated to all orders in α.
Since our interest resides in the low-momentum sector of fermion-fermion elastic scattering, we incorporate the Coulomb ladders in the amplitude of the process to all orders in the fine structure constant. For scalar photons, this amounts to replacing the free-fermion propagators in the bubble diagrams of fig. 1 with the Coulomb propagators in fig. 3). To this aim, we follow the procedure outlined in ref. [27] and introduce the Coulomb Green's functions. The inclusion of the Coulomb potential (cf. eq. (28)) in the Hamiltonian yields the Coulomb Green's function operator, an expression that, together with eq. (8), admits a self-consistent rewriting à la Dyson [100], that can be diagrammatically represented as in fig. 3.
Moreover, the solutions of Schrödinger equation with a repulsive Coulomb potential, , can be formally expressed in terms of the free ones as see eq. (18) in ref. [27]. The above eigenstates share with the plane waves the generalized normalization property, i.e. ψ (±) p |ψ (±) q = (2π) 3 δ(q − p). If the potential is repulsive, the solution with outgoing spherical waves in the future is given by while the state with incoming spherical waves in the distant past coincides with where M (a, b; c) is a Kummer function. In particular, the squared modulus of the two given spherical waves evaluated in the origin, i.e. the probability of finding the two fermions at zero separation, is equal to known as the Sommerfeld factor [101,102]. Since the scattering eigenfunctions of the repulsive Coulomb Hamiltonian form a complete set of wavefunctions, they can be employed in an operatorial definition of the Coulomb Green's functions analogous to eq. (9), In a way fully analogous to the one with which we have defined the Coulomb Green's functions in eq. (36), we introduce the full Green's function, including both the strong and the electrostatic interactions. Therefore, we add the operatorV S ≡V (1) to the kinetic and Coulomb potential in eq. (36), so thatĜ Then, we define the incoming and outcoming wavefunctions as in ref. [27], similar to the eq. (38). Exploiting the operator relation C (E) we find the self-consistent Dyson-like identitŷ that permits to rewrite the eigenstates of the full Hamiltonian in terms of the Coulomb states, Subsequently, the scattering amplitude can be computed via the S-matrix element, given by the overlap between an incoming state with momentum p and an outcoming state p , where T (p , p) = T C (p , p) + T SC (p , p) as in eq. (4) in ref. [103] (for the complete derivation of eq. (47) we refer to chap. 5 of ref. [104]). In particular T C (p , p) = p |V C |ψ (+) p is the pure electrostatic scattering amplitude and T SC (p , p) = ψ is the strong scattering amplitude modified by Coulomb corrections. Since the eigenstates ψ p of the former are known, the scattering amplitude due only to the Coulomb interaction can be computed in closed form and admits the following partial wave expansion [27], where θ is the angle between p and p and σ = arg Γ(1 + + iη) is the Coulomb phase shift. In particular, the strong scattering amplitude T (p, p ) possesses a phase shift σ + δ . Furthermore, the Coulomb corrected version of T S can be expanded in terms of the Legendre polynomials P as where δ is the strong contribution to the total phase shift. After expressing the eigenstates of the full Hamiltonian in terms of the Coulomb eigenstates (cf. eq. (46)), we concentrate on the P-wave amplitude. Since the strong interaction couples the fermions to one unit of angular momentum and Coulomb forces are central, the only nonzero component of T SC of the expansion in eq. (49) is the one with = 1. Analogously to eq. (31) in ref. [27], we can, thus, write and we replace the ERE of the l.h.s. of the last equation with the = 1 version (cf. ref. [37]) of the generalized effective-range expansion formulated in ref. [36] for the repulsive Coulomb interaction, where a (1) 0 and r 1 are the scattering length, the effective range and the shape parameter, respectively. By comparison with the S-wave counterpart in eq. (32) of ref. [27], we can observe that, apart from the different power of the momentum of the incoming particles in front of the cot δ 1 − i term, the most significant difference is provided by the polynomial on the l.h.s. of the eq. (51), containing all the even powers of η from zero to 2 , as shown in eq. (10.10) in ref. [37].
Besides, the function H(η), that represents the effects of Coulomb force on the strong interactions at short distances, is given by where ψ(z) = Γ (z)/Γ(z) id the Digamma function. Despite the appearance, the generalized ERE is real, since the imaginary parts arising from H(η) cancel exactly with the imaginary part in the l.h.s. of eq. (51). Due to the following identity on the logarithmic derivative of the Gamma function, in fact, the imaginary part of H(η) proves to coincide with C 2 η /2η. For the sake of completeness, in the case of fermion-antifermion scattering the Coulomb potential is attractive and H(η) in the effective range expansion (cf. eq. (51)) should be replaced bȳ where η = −αM/2p is defined as a negative real parameter.

Repulsive channel
Considering the results of the previous section, all the elements for the derivation of the Coulomb-corrected strong scattering amplitude, T SC (p, p ), are available. Recalling the definition and eq. (46), the amplitude can be computed by evaluating each of the terms in the expansion, whose insertions are given by retarded Coulomb propagatorsĜ (+) C followed by = 1 four-point vertices,V S =V (1) . In particular, the lowest order contribution to the T-matrix reads where eqs. (16) and (39)- (40) have been exploited, partial integration for the two variables has been performed and the vanishing surface terms dropped. The explicit computation of the two integrals over the free Coulomb wavefunctions (cf. eqs. (39)- (40)) in the last row is carried out in app. C, and yields where C 2 η is the Sommerfeld factor, a function of η = αM/2|p|. As it can be observed, the polynomial 1 + η 2 in the l.h.s. of the generalized effective range expansion appears, see eq. (51).
As shown in app. C, the P-wave strong vertex projects out of the integral all the components of the Coulomb wavefunctions ψ (±) p with = 1 appearing in the angular momentum expansion ψ (+) wherep andr are unit vectors parallel to p and r ≡ rr respectively and F (η, |p|r) is the regular free Coulomb wavefunction. The latter functions, derived by Yost, Wheeler and Breit in ref. [105], display a regular behaviour in the vicinity of the origin, in contrast with the G (η, |p|r), linearly independent solutions of the Whittaker equation for a repulsive Coulomb potential which are irregular for r → 0. Explicitly, F (η, |p|r) has the form given in ref. [102], whereas the expansion for the incoming waves is obtained via the complex-conjugation property Next, we proceed with derivation of the next-to-leading order contribution to the scattering amplitude, where partial integration has been exploited and Einstein summation convention over repeated indices is understood. More succintly, the last equation can be recast as where in the second row, the Coulomb-corrected counterpart of the J 0 matrix defined in eq. (14) has been introduced, Analogously to the = 0 case, the higher order contributions to the T-matrix possess the same structure of eqs. (56) and (60) differ from the latter only in the powers of J C and the coupling constant D(E * ). Therefore we can again write and we can treat the terms enclosed by the round barckets as a geometric series, Recalling the tensor product between vectors and the definition of the Coulomb Green's function operators in eq. (42), it is convenient to rewrite the overall matrix as and we observe that the numerator can be considerably simplified by means of the results of app. C. In particular, eq. (281) can be applied twice, yielding Equipped with the last result together with eq. (41), we recast the components of the J C matrix as where the dependence of η on the integrated momentum s has been made explicit. As the = 0 counterpart in eq. (43) of ref. [27], the integral is ultraviolet divergent. Additionally, all the offdiagonal matrix elements of J C vanish, as the integrand is manifestly rotationally symmetric in three dimensions except for the components s i s j , that are integrated over a symmetric interval around zero, see eq. (4.3.4) in ref. [106]. In dimensional regularization, eq. (66) can be rewritten as an expression that in three dimensions, combined with the results in app. C, allows to simplify the Coulomb-corrected strong scattering amplitude as in momentum space. Returning to eq. (63) and ignoring the Feynman prescription in the denominator, we first exploit the aforementioned trick and split the integral into three parts, While the first one proves to be finite, the other two display a pole for d → 3 and the PDS regularization scheme has to be implemented. We begin with integral in the first row of eq. (69), j fin C . The numerator of the latter can be split into two parts, according to the terms of the polynomial in η inside the round brackets. Taking the limit d → 3, we observe that one of the two parts coincides with J fin 0 in eq. (45) of ref. [27], up to a proportionality constant equal to p 2 /3. In comparison with the latter, the other part of j fin C in eq. (69) is suppressed by two further powers of η, therefore it is pairwise UV-finite and the three-dimensional limit finds a justification. After these manipulations, j fin C becomes j fin Due to spherical symmetry, the integration over the angular variables in the last term can be immediately done. By performing again the substitution s → 2πη = παM/s, the second term on the r.h.s. of eq. (70) can be simplified as where a ≡ iπαM/|p|. The first of the two integrals on the r.h.s. of eq. (71) can be evaluated by means of the following identity connecting Euler's Gamma function with Riemann's Zeta function, while the second one in eq. (70) is analogous to the integral in eq. (46) of ref. [27], modulo a constant factor. Considering the last two identities, eq. (71) can be recast into where the definition of H(η) in eq. (52) and the fact that ζ(2) = π 2 /6 have been exploited. The subsequent addition of the last result to the already calculated contribution to eq. (70) yields the sought closed expression for j fin Now we focus on the term in the second row of eq. (69). By comparison with the integrand of eq. (44) of ref. [27], we expect the integral of interest to display an UV singularity. Splitting the polynomial within the round brackets on the numerator of the integrand, we recognize, in fact, the already available J div 0 in ref. [27] whose result in the PDS scheme is given in eq. (53) of the latter reference, Again, spherical symmetry permits to integrate over the angular variables of the second integral on the r.h.s. of eq. (75) and the substitution s → 2πη = παM/s allows for the exploitation of the integral relation between the Gamma-and the Riemann Zeta function in eq. (72), obtaining Unlike the first term on the r.h.s. of eq. (75), the present integral proves to be convergent in three dimensions, since ζ(2) = π 2 /6 is finite and the argument of the Gamma functions are positive integers or half-integers. Additionally, no PDS poles are found in the same expression. Therefore, the limit d → 3 can be safely taken, yielding Plugging the available result in eq. (53) of ref. [27], we can finally write a closed expression for j div,1 SC (p) in the PDS regularization scheme, Finally, we concentrate our attention on the term in the third row of eq. (69). From that equation, we infer that the only difference with respect to integrand of j div,1 C consists in the absence of the factor 1/s 2 , which enhances the divergent behaviour of the integral in the s → +∞ limit.
Therefore, we expect also this third contribution to j C to be UV divergent. After splitting the integral as in eq. (75), we obtain Now we focus on the first term on the r.h.s. of the last equation. Rotational invariance allows again for the integration over the angular variables in d dimensions. Then, change of variables s → x ≡ 2πη permits to exploit again the multiplication identity between the Riemann Zeta and the Euler's Gamma functions (cf. eq. (71)). Additionally, thanks to the fundamental properties of the Gamma function and the definiton of where, in the last step, the Gamma functions and the physical constants have been rewritten in order to highlight the dependence on the small quantity . From the last row of eq. (80), we can infer that, while the Gamma function has a simple pole for d → 3, the Riemann Zeta function analytically continued to the whole complex plane is zero in that limit, since it is evaluated at a negative even integer, i.e. ζ(−2n) = 0 n ∈ N + . Therefore, the fourth expression in eq. (80) cannot be immediately evaluated in the three-dimensional limit. Performing a Taylor expansion of the Zeta function about −2, we obtain where negligible terms in have been omitted in the intermediate step.
As it can be inferred from eq. (82), the result of the integration becomes finite in the framework of dimensional regularization, even if the corresponding integral in the first row of eq. (80) is divergent for d = 3 due to the singularity at x = 0. Since the original expression in the end of the second row of eq. (80) contains a pole at d = 2 while Γ(1) = 1 and ζ(−1) = − 1 12 in the two-dimensional limit, the PDS correction should be taken into account. Therefore, the complete application of the PDS scheme into eq. (82) Next, we switch to the evaluation of the last term on the r.h.s. of eq. (79). Proceeding exactly as in eq. (80), we find Differently from the previous case, the Riemann Zeta function is nonzero in the three-dimensional limit and the only singularity for = 0 belongs to the Gamma function in the numerator of the last row of eq. (84). Considering the expansions of all the -dependent functions about zero, the asymptotic expression for eq. (84) is recovered As it can be inferred from eq. (84), also a PDS singularity at d → 2 is present, since the Riemann Zeta function displays a simple pole at unit arguments. In particular, the Laurent expansion of the Zeta function around 1 yields Applying the PDS regularization scheme and subtracting the correction corresponding to the Thanks to the last expression and eq. (87), a closed form for the third contribution to the diagonal elements of the J C matrix is found, Finally, collecting the three results in eqs. (74), (78) and (88), the latter matrix elements are A direct comparison with the = 0 counterpart of the last expression, eqs. (47) and (53) in ref. [27], shows that the QED contributions to j C include terms of higher order in the fine-structure constant α. Moreover, owing to the elements j fin C and j div,1 C , an explicit dependence on the momenta of the incoming fermions ±p outside H(η) appears. Since j C contains quadratic terms in p, eq. (89) gives rise to a non-zero value for the effective range parameter r  (69), an expression for |p| 3 (cot δ 1 − i) can be found, Plugging the last expression into the = 1 generalized ERE formula, the term of eq. (89) proportional to H(η) cancels out with its counterpart in eq. (51), and all the momentum-independent contributions can be collected, yielding the expression for the Coulomb-corrected = 1 scattering length, which represents the measured P-wave fermion-fermion scattering length. As in the = 0 case, the ultraviolet pole is expected to be removed by counterterms which describe short-distance electromagnetic and other isospin-breaking interactions due to the differences between the quark masses [107]. The subsidiary terms will transform the coupling constant D(E * ) into a renormalization mass dependent coefficient, D(E * , µ), which allows for a redefinition of the scattering length as in eq. (55) of ref. [27], 1 The latter quantity is non-measurable and depends on the renormalization point µ, related to the physical scattering length through the relation which is the = 1 counterpart of eq. (56) in ref. [27]. Besides, grouping the quadratic terms in the momentum of the fermions arising in the l.h.s. of eq. (51), an expression for the effective range is recovered, As in the case of the inverse of the scattering length in eq. (93), r 0 possesses a simple pole at d = 3. Now the energy-dependent coefficient of our P-wave interaction D(E * ) is replaced by D 0 , the singularity can be removed by means of counterterms coming from the p 2 -dependent interactions and the one with = 1 (cf. eq. (3)) potentials give rise to a scattering amplitude T S (p, p ) whose |p| 2 +1 · (cot δ − i) factor leads to a vanishing effective range. As soon as the Coulomb interaction is included in the Lagrangian, when the potential couples the fermions to one unit of angular momentum, a purely electrostatic non-zero effective range emerges, in contrast with the = 0 case, see sec. 3.3 in ref. [27]. Therefore, we shall expect that, for higher angular momentum interactions further coefficients in the generalized expansion of |p| 2 +1 cot δ in even powers of the momentum of the fermions in the CoM frame become non-zero when the colliding particles are allowed to exchange Coulomb photons.

Attractive channel
We consider the scattering of two non-relativistic fermions with opposite charges, such as fermion-antifermion pairs. Concerning elastic scattering, the continuum eigenstates are again represented by the spherical wave solutions in eqs. (39)- (40), with η now given by −αM/|p|.
Besides, the phenomenology of the scattering process is now enriched by the presence of bound states. In addition, annihilation is possible, but this will not be considered here. The Coulomb Green's function, in fact, is enriched by discrete states, φ n, ,m (r), corresponding to bound states with principal quantum number n ≥ 1 and rotation group labels given by ( , m), where E s is equal to α 2 M/4η 2 and the bound state eigenvalues, E n , are given by Bohr's formula for a system with reduced mass equal to M/2, in natural units. As in the previous case, the Coulomb-corrected strong scattering amplitude of the elastic scattering process in configuration space takes the form where D(E * ) is the strong P-wave coupling constant in presence of attractive electrostatic interaction and the matrix J C is, now, given by which corresponds to the addition of the contributions from discrete and continuum states, and respectively. Let us start by evaluating the term J d C . With reference to the expression of the eigenfunctions belonging to the discrete spectrum, where L n k (x) are the associated Laguerre polynomials, we first evaluate the integrals containing the gradient of the latter in the expression forJ d C in eq. (99), that can be performed separately for each of the wavefunctions, since the denominator does not depend on the coordinates. The application of the gradient on the bound state wavefunctions, ∇φ n, ,m (r) where the spherical symmetry of the Dirac delta has been exploited. Of the latter equation, we consider now the first term on the right hand side. Firstly, expressing the radius vector componentwise as a spherical tensor of rank 1 (cf. eq. (5.24) and sec. 5.1 in ref. [108]), the aforementioned part of eq. (102) becomes Now, recalling the expression of the constant term of the associated Laguerre polynomials, eq. (103) can be concisely recast into where the integration over the angular variables Ω has been performed. After replacing the , and performing few manipulations, the sought expression is recovered, Concerning the second term on the r.h.s. of eq. (102), the rewriting of the gradient of a spherical harmonic into linear combination of spherical tensors (cf. eqs. (5.24) and (5.27) in ref. [108]) gives thus, allowing again for an immediate integration over the angular variables, where the eq. (104) for the evaluation of the Laguerre polynomials at the origin has been exploited.
Subsequently, the replacement (011|0mm) = 1 gives the desired expression for the second term of eq. (102), Equipped with the results in eqs. (106) and (109), the original integral can be immediately evaluated, Now, taking the tensor product of the latter expression with its complex-conjugate version, as required by eq. (99), J d C reduces to Since the diagonal form of the matrix in the spherical complex basis (cf. eq. (2.141) in ref. [108]) is preserved in the Cartesian basis and the sum over the principal quantum number can be decomposed and evaluated in terms of the Digamma function ψ(z), the contribution to the scattering matrix due to the discrete states, can be ultimately rewritten as As underlined in sec. 3.4 of ref. [27], the divergent sum of the harmonic series, ζ(1), appears in the last formula. Its presence is only due to the numerable infinity of states in the discrete spectrum, whose energy depends on the inverse square of n, while the modulus square of the gradient of the eigenfunctions evaluated at the origin yields a factor ∝ n 2 . The replacement of ζ(1) in eq. (115) by its Cauchy principal value, equal to γ E , allows to assign a finite value to j d C and, thus, circumvent the divergence.
At this stage, we switch to the continuous contribution to the auxiliary scattering matrix, J c C . As for the repulsive counterpart in sec. 2.2, the possible divergences in the three-dimensional limit require the rewriting of the relevant intergrals in arbitrary complex dimension d. Therefore, the dimensionally regularized version of the second term on the r.h.s. of eq. (98) gives where the initial integral has been split into two parts, making use of the trick in eq. (69). Due to the sign change in η, the first term on the r.h.s. of eq. (116), can be immediately evaluated, since it coincides with eq. (66). Therefore the result in eq. (89) can be directly exported, rewriting eq. (117) as j old where the H(−η) function has been replaced by its definition in terms of the Digamma function in eq. (52), in sight of the next developments. Subsequently, we evaluate the second term on the r.h.s. of eq. (116), the new part of the continuum states contribution. In order to bring s 2 to the denominator, we apply again the trick introduced in eq. (69) and split the integral into three parts, Concerning the first term on the r.h.s. of the latter equation, it vanishes in dimensional regularization, see eq. (4.3.1a) in ref. [106]. Therefore, we can switch to the subsequent term of eq. (119) and apply Feynman's trick for denominators, finding Defining again the auxiliary variable γ = −i|p|, we perform the momentum integration in eq. (120), Then, since the remaining integration over ω turns out to be finite in two dimensions and the rest of the expression does not display any PDS singularity, we can directly reintroduce ≡ 3 − d and consider the three-dimensional limit. In particular, the integral over ω in eq. (121) can be evaluated in first-order approximation in , obtaining Second, the terms depending on in the exponent can be grouped and expanded to first order in as in eq. (82), whereas the Gamma function can be expressed in Laurent series up to order 0 . Performing few manipulations and taking the → 0 limit, the original expression in eq. (120) becomes Subsequently, we compute the last term on the r.h.s. of eq. (119). As it can be inferred, the integral coincides with the one in of eq. (123), except for an overall factor of α 2 M 2 /4p 2 = η 2 .
Therefore, its evaluation is straightforward and gives Collecting both the results in eqs. (123) and (124), we obtain the sought expression for j new We now collect all the contributions in eqs. (115), (119) and (125) and write a closed form for the diagonal matrix elements of J C , where the definition of H(η) in eq. (54) has been exploited and the Cauchy principal value of ζ(1) has been taken. A direct comparison with the repulsive counterpart of the last formula in eq. (89) shows that the map between the two expression is provided by the sign reversal in front of all the terms containing odd powers of the fine-structure constant and the replacement of H(η) by H(η). This fact is consistent with the conclusions drawn from eq. (70) in ref. [27], where all the PDS-corrective terms remained unaffected by the sign change in the charge of one of the interacting fermions.
We conclude this section with the derivation of an expression for the scattering length and the effective range, by making use of the attractive counterpart of the generalized effective-range expansion in eq. (51), obtained by replacing again H(η) by H(η) with η < 0. As a consequence of the attraction of the electrostatic interaction, the Coulomb corrections in the strong scattering parameters change sign, consistent with eq. (126). Concerning the scattering length, we have where the divergence can be reabsorbed by the P-wave strong coupling constant. Analogously to eq. (72) of ref. [27], the renormalized version of the scattering length, a C (µ), can be defined in terms of the physical one, a Finally, the terms proportional to the square of the momentum of the fermions p give rise to a nonzero value for the effective range, as in eq. (94), whose divergent part, in case the energy-dependent coefficient of the = 1 interaction D(E * ) is replaced by D 0 , can be again reabsorbed by counterterms coming from p 2 -dependent = 0 interactions.

The Lattice environment
At this stage, we transpose the physical system of non-relativistic spinless fermions interacting via Coulomb photons onto a cubic lattice with L points per (spatial) dimension and spacing a. In this environment, it is customary to continue analytically the fields and the wavefunctions outside the lattice by means of periodic boundary conditions (PBCs). It follows that a free particle subject to PBCs carries a momentum p = 2πn/L, where n is a dimensionless three-vector of integers. Unlike QCD fields, the photon field in QED is truncated and modified by the boundary of the volume. In particular, when PBCs are implemented, the validity of Ampère's and Gauss's law is compromised. The problem is circumvented by introducing an uniform background charge density, a procedure that proves to be equivalent to the removal of the zero modes of the photon [88]. Once the latter are canceled, the Coulomb potential between two identical charges e becomes (cf. fig. 4) where the n ∈ Z 3 encodes the dimensionless lattice momenta. Discarding the zero modes, the latter are restricted to |p| ≥ 2π/L, whereas the viability of a perturbation treatment of QED is controlled by the parameter η = αM/2|p|, which scales as the inverse of the momentum of the interacting particles. Combining the above constraint with the definition of η, it follows that η ∼ αM L and the photon field insertions can be treated perturbatively if M L 1/α. As η grows linearly with the spatial volume, for any value of M exists a critical value of L that regulates the applicability of perturbation theory. Besides the condition η 1, we assume henceforth the limit M 1/L, since for the current Lattice QCD calculations large volumes are employed [88]. Furthermore, the finite volume QED effects are such that the energy eigenvalues of two charged fermions (e.g. hadrons) are modified in the same way by their self-interactions and by their interactions with each other, and the shifts take the form of power laws in L [40]. As a consequence, in the presence of Coulomb photons the kinematics of two-body processes receives power law modifications in the finite volume context [40,47]. In particular, if the infinite-volume effective range expansion for the P-wave scattering is rewritten in terms of the center of mass energy, then E * = 2M + T in the above expression is replaced by its finite volume counterpart 2 . Eq. (131) thus becomes The original dependence of the r.h.s. of the last equation on the powers of the finite volume kinetic energy T L = E * L − 2M L can be exactly restored by exploiting the expression of the finite-volume shift for the masses of spinless particles with unit charge in eqs. (6) and (19) of ref. [40], where the sum of the three-dimensional Riemann series regulated by the spherical cutoff Λ n is denoted with I (0) ≈ −8.913632 (cf. app. D.1). To this purpose, primed scattering parameters are and where the changes in the total energy have been incorporated in the primed scattering parameters.
Finally, also the validity region of the last expansion is modified by the lattice environment, due to the changes in the analytic structure of the scattering amplitude in the complex |p| plane. The absence of the zero mode in the Coulomb potential in eq. (130), in fact, yields a shift in the branch cut of the imaginary |p| axis from the origin to 2πM/L + O(1/M ), which fixes the inelastic threshold for the two-hadron state (cf. fig. 2 in ref. [88]).
The last version of the ERE, combined with the quantization conditions discussed below, will turn out to be the key ingredient for the derivation of the finite volume energy corrections for scattering and bound states with one unit of angular momentum.

Quantization Condition
After introducing the finite and discretized configuration space, we derive the conditions that determine the counterpart of the = 1 energy eigenvalues on the lattice. These states transform as the three-dimensional irreducible representation T 1 (in Schönflies's notation [109]) of the cubic group [56,57,58], the finite group of the 24 rotations of the cube that replaces the original SO (3) symmetry in the continuum and infinite volume context [55].
As it can be inferred from eq. (43), the eigenvalues of the full Hamiltonian of the system H 0 +V C +V S can be identified with the singularities of the two-point correlation function G SC (r, r ) and are called quantization conditions in the literature [63,64,65,88]. The Green's functions G SC (r , r) in turn can be computed from the terms in the expansion over the P-wave interaction insertions stemming from eq. (45), withV S in momentum space given in eq. (6). In particular, the three lowest order contributions in D(E * ) yield, respectively, r |Ĝ and r |Ĝ Extending the calculation to higher orders, the expression of (N + 1) th order contribution to the full two-point correlation function can be derived, thus allowing to rewrite the original Green's function in terms of a geometric series of ratio , that we identify as J C (cf. eq. (64)), where can be similarly interpreted as a source and a sink coupling the fermions to an P-wave state, respectively. As in the = 0 case, the pole in the second term of eq. (144) permits to express the infinite volume quantization conditions, where the identity matrix multiplied by the CoM energy dependent coupling constant is equal to the inverse of the matrix of the double derivatives of the Coulomb two-point Green's function evaluated at the origin, J C . Concentrating again on the retarded two-point correlation function and adopting the notation of ref. [88], the finite-volume counterpart of eq. (144) becomes . (146) Similarly, the pole in the second term on the r.h.s. of the last equation yields the finite volume quantization condition, that determines the T 1 eigenvalues. As in the = 0 case, we proceed by expanding G  fig. 2 to two-fermion vertices, evaluated at the origin in configuration space. As a consequence, the matrix elements of J L C can be interpreted as bubble diagrams with multiple Coulomb-photon insertions inside. Analytically, the two lowest order contributions to J C , that correspond to bubble diagrams, respectively with and without a Coulomb-photon insertions, read where the Dyson identity betweenĜ (±) C andĜ (±) 0 (cf. eq. (37)) has been exploited and ε has been set to zero. Replacing again the integrals over the momenta by sums over the dimensionless momenta n, m ∈ Z 3 , the lattice counterpart of eq. (148) is obtained, where the finite-volume mass M L , the speherical lattice cutoff Λ n and the dimensionless CoM momentum of the incoming particlesp have been reintroduced. With the aim of regulating the sums in eq. (149) for numerical evaluation while maintaining the mass-independent renormalization scheme (cf. sec. II B of ref. [88]), we are allowed to rewrite the finite volume quantization conditions as where J where S 2 Λ denotes the three-dimensional sphere with radius Λ. Isolating the O(α) contribution, we obtain where the isotropy of the cutoff has been exploited in the second step and O(Λ 0 ) denotes constant or vanishing terms in the Λ → +∞ limit. Concerning the O(α) term, the integral can be simplified as follows where Ξ 1 ≡ (1 − ω)q 2 − ωp 2 and Feynman parametrization for the denominators has been applied. The subsequent integration over the momentum k through eq. (B17) in ref. [98] and the exploitation of rotational symmetry in the outcoming integrand (cf. eq. (4.3.1a) in ref. [106]) gives Then, it is convenient to split the integrand of eq. (154) into two parts and to simplify the numerator, where J 1 (J 2 ) corrresponds to the first (second) integral on the l.h.s. of the last equation and it will generate the leading contributions in Λ. Beginning with the latter term, integration over the momentum l yields where Ξ 2 is an ancillary variable, .
Exploiting the fact that p/Λ 1, the integrand in the latter expression can be considerably simplified. Exploting the results the remaining integration can be performed, obtaining the desired expression for I 1 Concerning the I 2 term, its rotational symmetry and integration over the angular variables permits to split it in turn into two terms, Considering the first term on the r.h.s. of eq. (160), integration over the radial momentum yields again an arccoth(x) function, which is eventually responsible of a further logarithmic divergence in the UV region, Approximating the expression again under the assumption p/Λ 1 and performing the integration over ω (cf. eq. (157)), the expression on the r.h.s. of eq. (161) becomes i.e. it carries the second logarithmic contribution to J {Λ} C to order α in the perturbative expansion. Finally, we consider the second term on the r.h.s. of eq. (160) and introduce the auxiliary variables γ 2 ≡ −p 2 and Ξ 3 ≡ γ 2 /(1 − ω). The integration over the radial momentum |q| yields an expression that can be simplified in the large coutoff limit, p/Λ 1, obtaining Although the remaining integral is unbound, the overall expression is independent of the cutoff where the O(Λ 0 ) contributions have been discarded. Now we proceed with the calculation in dimensional regularization of J C . To this purpose, it is convenient to start from the exact expression of J C to all orders in α in arbitrary d dimensions (cf. eq. (65)), where the integrand has been expanded up to the first order in α. In particular, the α-independent contribution in eq. (166) gives where the rotational invariance of the integrand has been exploited and ε has been set to zero.
Since the integrand is a polynomial in the momentum, the first contribution on the r.h.s. of the last equation vanishes in dimensional regularization, whereas the remaining term turns out to coincide with the purely strong counterpart of J C (d; p) in eq. (18), thus is finite and the limit d → 3 can be safely taken, obtaining where γ ≡ −i|p| and the Feynman parametrization for the denominators has been adopted. The subsequent momentum integration in the latter yields where the renormalization scale µ has been introduced. Conversely, the first term in the second row of eq. (170) vanishes like the first integral on the r.h.s. of eq. (167). Introducing the small quantity = 3 − d, the integral over ω can be computed to first order in , obtaining Exploiting the result in the last formula, eq. (170) partially expanded to order becomes Expanding in turn the Gamma function and the power term µ √ π/γ in Laurent and Taylor series, respectively, and truncating the expansion to order 0 , the desired expression for T SC (p) in dimensional regularization is recovered Then, we bring further simplification to the finite volume quantization condition by taking the trace of eq. (150), thus transforming a matrix identity into a scalar one as the one for S-waves. Finally, we take the real part of the expressions in eqs. (165) and (173) and replace the fermion mass by its lattice couterpart, obtaining the regulated version of the finite volume quantization condition (cf. eq. (147)) in explicit form, Thanks to the last identity, the lattice version of the effective range expansion in eq. (139) can be rewritten in an explicit form. Besides, we drop in the next developments of the derivation the pole 1/ arising from dimensional regularization, since it does not deliver information on the energy eigenvalues of our two-body system.

Finite Volume Effective Range Expansion
First of all, we concentrate on the infinite volume version of the effective range expansion. By inserting the expression of T SC given in eq. (66) into eq. (51) and exploiting the closed form of j C in dimensional regularization to all orders in α given in eq. (89) we obtain a relation between the P-wave phase shift and the strong coupling constant D(E * ) in the presence of Coulomb photon exchanges, Since the asymptotic behaviour of the momentum integrals in the ultraviolet region is left invariant by discretization, eq. (177) can be straightforwardly adapted to the cubic lattice case, provided the infinte volume parameters are replaced by finite volume ones, which is valid to all orders in the fine-structure constant. Now, the quantization conditions derived in the previous section can be exploited by replacing the inverse of the finite-volume strongcoupling constant with the expression in eq. (175). Thus, the following equation is obtained, Similarly to the = 0 case, we notice that the finite-volume mass of the fermions is multiplied where the ellipsis stands for higher-order scattering parameters. Combining the latter expression with eq. (139) and isolating the regulated sums, we obtain the desired explicit version of the effective-range expansion to order α, in terms of the Lüscher functions, and In particular, thep-independent Lüscher sum in eq. (182) vanishes, as shown in app. D.1, whereas S 3 (p) generates the new series of double sums to be computed in the low-p limit.

Approximate Energy Eigenvalues
Since the Sommerfeld factor is not a rational function of the momentum of the colliding particles in the CoM frame, a non-perturbative counterpart in α of the eq. (181) in the lowmomentum limit would allow only numerical solutions for p 2 , which lie beyond our purpose.
Nevertheless, under the hypothesis that the expansions are perturbative in 1/L times the length scale characterizing the strength of the interaction, governed by the scattering parameters, and assuming that M L 1/α, the Coulomb photon insertions in the diagrams can be treated perturbatively. Under these conditions, the approximate expression of the ERE presented in eq. (139) can be exploited for an analytical derivation of the finite volume corrections to the energy eigenvalues.

The Lowest Unbound State
Differently from the S-wave case in sec. III D of ref. [88], the perturbative expansion of the arguments of the summations in the Lüscher functions around zero lattice momentump, corresponding to a total energy equal to 2M now looses significance, due to symmetry reasons. Given that = 1 states in the continuum and infinite volume are mapped to T 1 states in a finite cubic lattice, the latter are expected to be three-fold degenerate. While the multiplicity of the zero energy scattering state is one, its nearest neighbour with |p| = 1 is six-fold degenerate with total energy equal to 2M + 4π 2 /M L 2 , thus making it a suitable candidate for both T 1 and where the dots denote terms of order (δp 2 ) 6 and the large Λ n limit is understood. In the notation of app. D and ref. [88], eq. (185) is concisely recast into where the sums of the implied three-dimensional Riemann series are reported in app. D.1. Regarding the function S 2 (p), we proceed by isolating and expanding the double sums with |n| or |m| equal to one, The subsequent expansion of the argument of the two series in the same fashion of eq. (185) shows that all the resulting sums converge in the infinite-Λ n limit. Therefore, we are allowed to remove the spherical cutoff and the two double sums in eq. (187) merge, yielding where, differently from the ones in the appendix of ref. [88], the series χ i with i ≥ 1 include n = 0 and are defined as Furthermore, the expansion of the remaining term of the Lüscher function S 2 leads to results in analogous to the ones of the S-wave case, where the ellipsis stands for terms of order (δp 2 ) 3 . The R (1) and R (1) ij sums in the last equation coincide with the P-wave counterparts of the sums in eqs. (A6) and (A9) in ref. [88], namely and with i, j ≥ 2 and i+j ≥ 6, which are invariant under permutation of the lower indices, R ji , and convergent for i, j ≥ 1 and i+j > 2. Exploiting the symmetry property under index exchange and combining the results in eqs. (191) and (192), we find the desired result Subsequently, we treat the genuinely new Lüscher function, S 3 (p). To this purpose, we decompose the initial double series into three pieces, where a subsidiary spherical cutoff Λ m has been introduced in order to highlight the divergent nature of the three terms. In particular, the first and the second contribution to S 3 (p) on the r.h.s. of the last equation, can be recast as where the dots denote terms of order (δp 2 ) 3 , while the non-symmetric and divergent generalizations of χ 0 in eq. (283) and of R (1) ij in eq. (192) have been introduced, and respectively. Quite similarly, the third contribution to S 3 (p) in eq. (194) can be subdivided and expanded as follows where the two terms involving single sums can be, in turn, expanded in pairs, obtaining where the neglected terms are again of order (δp 2 ) 3 and the for the resulting convergent series the limit Λ n , Λ m → +∞ is understood. On the other hand, this limit can not be taken for the non-regularized counterpart of I (1) , whose divergence will cancel with the one from P (1) in eq. (196). Secondly, the expansion of the double sum in eq. (198) yields the appearance of further P (1) 2i 2j 2k terms, Now, the two partial results in eqs. (195) and (202) can be summed together, obtaining the sought expression of S 3 (p) as a power series in δp 2 . In particular, we notice that the sum of all the series appearing at each order in δp 2 has to be finite in the limit Λ n , Λ m → +∞, irrespective of the convergent or divergent behaviour of each individual sum. In particular, we observe that the latter limit can be directly taken to all orders in the small quantity, with the only exception of the sums in the δp 2 -independent contribution, that are regularized quadratically in the cutoff Λ n . As a consequence, we can simplify our δp 2 -expansion for S 3 (p) by grouping the divergent sums order by order and defining the finite coefficients and where the polynomials q S (n, m) and q X (n, m) are defined as and q X (n, m) = 3 − 3(|n| 2 + |m| 2 ) + |m| 2 |n| 2 + |m| 4 + |n| 4 .
After redefining the regularized sum in eq. (193) through the addition of the last term on the r.h.s. of eq. (181),R where the shifted higher-order scattering parameters r 2 (1) and r of ref. [88], modulo an overall damping factor ξ ≡ 4π 2 /L 2 . If the scattering parameters are small and of the same magnitude of 1/L as in the S-wave case, the importance of the auxiliary parameters just introduced can be quantitatively assessed. In particular, by assigning one unit of 'weight' for each scattering parameter in the effective-range expansion and one unit for 1/L, we find that, neglecting the fine-structure constant, the largest parameter is d 1 (order three), followed by d Due to the smallness of δp 2 , contributions multiplied by higher positive powers of the lattice momentum are increasingly suppressed. It follows that the dominant finite volume corrections are expected to be found by solving the truncated version of eq. (213) to order zero in δp 2 , Solving the last equation for δp 2 and expanding the denominator up to order twelve in the small constants, we find Denoting the outcome of the first iteration as δp 2 0 , we proceed with an improvement of the last result by plugging the latter into the term proportional to δp 2 of the expansion in eq. (213) truncated to the terms quadratic in the squared momentum shift, Solving the last equation and retaining only terms up to order twelve in 1/L and the scattering parameters, we obtain a refined version of the momentum shift, in which the new contributions appearing at each order in the auxiliary constants (cf. eq (215)) have been isolated. Restoring the dimensional constants in the momenta, the energy of the lowest T 1 scattering state can be now obtained in few steps, (1) where the differences with respect to eq. (218) involve only the linear term in the fine-structure constant and are proportional to I (0) .

The Lowest Bound State
Although bound states between two hadrons of the same charge have not been observed in nature, at unphysical values of the quark masses in Lattice QCD such states do appear [110,111,112,113]. Moreover, two-boson bound states originated by strong forces are expected to explain certain features of heavy quark compounds. In particular, the interpretation of
Additionally, loosely bound binary compounds of hadrons appearing in the vicinity of a P-wave strong decay threshold are not forbidden by the theory of hadronic molecules [115]. Possible candidates of such two-body systems are represented by the hidden charm pentaquark states P + c (4380) and P + c (4450), located slightly below theDΣ * c andD * Σ c energy thresholds at 4385.3 MeV and 4462.2 MeV, respectively. Although a wide variety of different studies on the two states have been conducted [116,117,118], a very recent one advances the molecular hypotesis [119] with orbital angular momentum equal to one in the framework of heavy quark spin symmetry (HQSS).
Therefore, it remains instructive to study the lowest two-fermion T 1 bound state, by switching to imaginary momenta p = iκ, where κ = |κ| represents the imaginary part of the momentum. To this purpose, we rewrite the FV effective range expansion in eq. (181), truncated on the l.h.s. to the sextic term in the binding momentum, First, we consider the limit of large lattice binding momentum,κ = |κ| 1, which corresponds to a tightly-bound state. Thus, an approximation for the Lüscher functions in this regime becomes necessary. In particular, we observe that the asymptotic behaviour of S 1 (iκ) and S 2 (iκ) is already available in literature and gives, and see eqs. (43,44) in ref. [88], respectively. On the other hand, the large binding energy limit of S 3 (iκ) requires a new derivation, presented in detail in app. (E), Recalling the fact that S 0 (κ)) = Υ = 0 (cf. eq. (283)) and collecting the results in eqs. (221)-(223), the above finite volume effective range expansion becomes, Second, we highlight the dependence on the fine-structure constant in the last equation by rewriting the binding momentum in a power series, Grouping all the terms independent on α, we observe that the following equality holds, since the remaining terms in eq. (226) depend linearly on α. Therefore, an expression for κ 1 can be drawn from the original eq. (226) by dropping the terms listed in eq. (227), In particular, retaining only the terms depending on the two lowest order scattering parameters as in the zero angular momentum case (cf. eqs. (45)-(46) in ref. [88]), a more approximated expression for κ 1 in terms of κ 0 , a C and r (1) 0 can be obtained, From the latter equation the first two terms generate finite-volume corrections, whereas the third one introduces QED modifications to the unperturbed binding momentum κ 0 which do not vanish in the infinite volume limit. Now, considering the binding energy of the lowest T 1 state in the linear approximation in α, and substituting the simplified expression of κ 1 in eq. (229), we find an approximate expression for the energy of the lowest T 1 bound state,  Also noteworthy is the fact that, when compared to the S-wave case in eq. (46) in ref. [88], the sign of the P-wave shift is reversed while the magnitude remains unchanged. As discussed in refs. [60,66], the significance of this behaviour can be traced back to the spatial profile of the = 0 and = 1 two-body wavefunctions associated to the considered bound eigenstates.
Qualitatively, the relationship found between the two finite volume energy corrections means that zero angular momentum states are more deeply bound when embedded in a finite volume, while the counterpart with one unit of angular momentum are less bound. In conclusion, together with the derivation of ∆E , we have simultaneously proven that the addition of a long-range force on top of strong forces in two-fermion systems within a cubic lattice produces changes of the same magnitude on S-and P-wave bound energy eigenvalues.

Outlook
Before drawing the conclusions of our work, we qualitatively outline the possible improvements to the above analysis. There are two main directions for the generalization of the present study.
These consist in the inclusion of transverse photons in the Lagrangian density of the system and in the treatment of strong interactions with higher angular momentum couplings.
Let us start from the former. Thanks to their vector nature, transverse photons can couple to the fermionic fields in several different ways, see eq. (26). Among these, here we retain only the Dipole vertex (cf. app. A), which is expected to yield the leading-order contributions to the T-matrix and the full Green's function. As we previously hinted, these photons, denoted by wavy lines, can propagate also between different bubbles. Consequently, at each order in D(E * ) a numerable infinity of new diagrams with different topologies and transverse-photon exchanges inside and outside the bubbles appears. Unfortunately, part of the amplitudes associated to these diagrams can not be written as powers of the loop (bubble) integrals, since the transverse photons propagating outside the fermion loops introduce a correlation between the bubbles (cf. fig. (4) in ref. [88] and fig. 5). It follows that an expression for the T-matrix element of the two-body scattering process, T SCT , written in terms of a geometric series of ratio proportional to the interaction strength, D(E * ), can not be found. More formally, a self-consistent rewriting of the full Green's function operatorĜ SCT in the form of a self-consistent equation à la Dyson separating the QED Green's function operator,Ĝ CT , from the strong interaction operator V (1) (cf. eq. (45)) does not exist. These facts prevent the exact determination of T SCT to all orders in the fine-structure constant.
Nevertheless, the gapped nature of the momentum in the lattice environment allows for a perturbative treatment of the whole non-relativistic QED. Therefore, approximate expressions for T SCT that incorporate the effects of the transverse photons up to the desired order in α can be derived.
One of these approaches consists in writing the infinite-volume T-matrix element T SCT exactly as T SC in eq. (69) with D(E * ) at the denominator replaced by a dressed strong P-wave coupling constant D T (E * ), that includes the effects of transverse photons up to first order in the finestructure constant. Analytically, this energy-dependent constant can be derived by evaluating the contributions of all the possible bubble diagrams with one tranverse-photon exchange in fig. 5.
Moreover, considering the fact that diagrams with one transverse photon across n-bubbles are suppressed by a factor ( |p|) n , the numerable inifnity of contributions on the r.h.s. of fig. 5 can be reduced to a finite set. The amplitudes corresponding to these diagrams can be evaluated via dimensional regularization as the ones containing radiation pions in refs. [14,24,120] or via the cutoff approach and, in finite volume, they can be constructed by replacing the relevant integrals with three-dimensional sums, eventually regularized by a spherical cutoff. In the present linear approximation in the fine-structure constant, the latter operation is, in fact, justified. Due to their dependence on the scattering parameters, the new terms arising from the transverse-photon interactions are not expected to bring the dominant O(α) contibution to the finite-volume corrections of the lowest T 1 scattering state. Additionally, such terms may have no effect in the final expression for the corrected version of ∆E As we previously noted, the second main generalization of our work consists of the adoption of strong interactions coupled to more units of angular momentum. The most significant of these extensions is represented by the D-wave case, where the strong part of the EFT Lagrangian density becomes where F (E * ) is a new suitable energy-dependent coupling constant (cf. sec. 2). Due to the presence of higher-order differential operators acting on the fermionic fields, the computation of the strong scattering amplitude, T S , via a geometric series on the loop integrals, as in eq. (15), involves a rank-four tensor as a ratio. The elements of the D-wave counterpart of J 0 correspond to double mixed derivatives of the free Green's function, except for the diagonal terms, in which an additional contribution proportional to G 0 (0, 0) is expected to be present. Also the full Green's function G SC (0, 0) is likely to undergo similar changes, which lead to the onset of more rapidly UV-divergent integrals in the D-wave counterpart of J SC . Besides, some novelties are expected to arise from the quantization condition stemming from the full two-point correlation function As a consequence, the one-to-one correspondence between the transformation properties of the selected multiplet of states under the operations of SO (3) and O is no longer valid. It is, thus, possible that a derivation of the finite-volume corrections that makes use of the effective range expansion for D-waves in ref. [37], has to follow two separate paths for the E and the T 2 states. New challenges arise also in the

Conclusion
In this paper, we have extended the investigation of QED for low-energy scattering to the case of spinless fermions strongly coupled to one unit of angular momentum in the framework of pionless EFT. Among the forces of electromagnetic nature, the electrostatic interaction represents the dominant contribution to T-matrix elements in the low-energy sector and the Coulomb ladders have to be resummed to all orders in the fine-structure constant.
A pivotal role in the procedure has been covered by the non-perturbative formalism based on the full Coulomb propagator developed comprehensively in ref. [27] in the context of pp scattering.
Differently from the transverse photons, the Coulomb ones do not propagate between the fermionic bubbles in the diagrams, a crucial feature that permits to rewrite the full two-body Green function operator in terms of the Coulomb one and operator representing the P-wave strong interaction. It is exactly the same property that led us to write the T-matrix element of the scattering processes and the full Green's function in closed form, thus paving the way to the derivation of the quantization conditions. Moreover, though the analysis of the attractive Coulomb case, the observations on the scattering parameters pointed out in sec. 3.4 of ref. [27] have found here another confirmation.
Second, the infinite-volume analysis of fermion-fermion scattering in secs. 2-2.3 allowed us to attain our main goal, the derivation of finite-volume energy corrections for two-body P-wave bound and scattering states, by providing an extension of the analysis on S-wave states in ref. [88].
Motivated by the growing interest for lattice EFTs and, above all, LQCD, we have transposed our system of charged particles in a cubic box with periodic boudary conditions.
Having regards to the prescriptions from the literature [87,91], we have removed the zero momentum modes from the relevant three-dimensional sums and considered the QED corrections to the mass of the spinless particles [40,47], in sight of the application of our results to realistic baryon-baryon systems on the lattice. Furthermore, the characteristic size of our cubic box has been chosen to fulfill the constraint M L 1/α, which is required for the viability of the perturbative approach of QED on the lattice. Under this hypothesis, the non-relativistic relation between the finite-volume energy of two composite fermions in the T 1 representation of the cubic group and its P-wave scattering parameters receives QED corrections obtainable in closed form.
Although more cumbersome than the S-wave counterpart, the expression we have presented in secs. 3.3.1 for the energy shift of the lowest unbound state resembles the features of the one in sec. III D 1 of ref. [88], except for the appearance of higher-order scattering parameters. Besides, the finite-volume corrections for the P-wave bound state prove to have the opposite sign and the same magnitude of the ones for the S-wave state in sec. III D 3 of ref. [88], up to contributions of order 1/L 3 . This fact confirms the long-standing observations on the A 1 and T 1 lattice energy eigenvalues in the analysis of refs. [60,66], drawn in the context of a two-body system governed by finite-range interactions in the non-relativistic regime.
In the latter work, the interplay of parity and angular momentum quantum numbers in the wavefunctions was found to be responsible of the relation between the leading-order S-and P-wave energy shifts. Only the generalization of our analysis will tell if the existing relationships between the finite-volume shifts in tab. I of ref. [60] for two-body states with higher angular momentum remain at least approximately valid in presence of QED.

Acknowledgements
First of all, we express our gratitude to Akaki Rusetski, Andria Agadjanov and Wael Elkam-

Appendix B Integrals in Dimensional Regularization
In this section the evaluation of the second term on the r.h.s. of eq. (30) for T tree SC (p, p ) is explicitly shown, that is the derivation of the new contribution to eq. (31). In d dimensions, the last integral in eq. (30) becomes Since the loop integral is not going to be performed in the complex plane, the parameter can be set to zero. Applying the Feynman parametrization for the denominators, we obtain Then, we rewrite the polynomial in the denominator as l 2 − 2p · l(1 − ω) + ωλ 2 and consider the application of the dimensional regularization formula in eq. (B.17) in ref. [98] for the carrying-out of the momentum integral, As the r.h.s. of the last equation is non-singular in three dimensions, the limit d → 3 can be safely taken. Furthermore, we define the auxiliary quantities Performing the change of variables ω → ω ≡ (ω − β)/γ in the integrand of eq. (237), we obtain The last expression can be integrated, After evaluating the results of the integration over ω in terms of the original variables, p and λ, a closed form for K is found, where the conventions log(−i) = −iπ/2 and √ −1 = −i are understood. Plugging eq. (241) into eq. (237), the desired result is obtained, In the second part of this appendix, we evaluate the momentum integrals appearing in the 1-loop diagram with one Coulomb photon exchange inside the loop (cf. right part of the fig. 2 and eq. (34)). To this aim, we first introduce the notation γ 2 ≡ −p 2 for the physical momenta, we set ε to zero and rewrite eq. (33) in arbitrary dimensions d as where the fictitious photon mass has been set to zero, since no infrared divergences occur in the integration (cf. Appendix A of ref. [27]). Secondly, we decide to carry out the integration over k first and merge the relevant denominators again by means of Feynman's trick, obtaining where ∆ 2 ≡ γ 2 + q 2 (1 − ω) and Einstein's convention for repeated indices is understood. Now, we can proceed by evaluating the two integrals over k generated by q i k j and k i k j separately. The former integral turns out to be an application of eq. (B.18) in ref. [98] and gives .
The second term of eq. (244), instead, yields two contributions, one of the two being proportional to the r.h.s. of eq. (245). The application of eq. (B.18) in ref. [98] indeed leads to By comparison of the last equation with eq. (244), we observe that the q integrals involving q i q j can be performed together. Therefore, we group the two terms together with the necessary constants and, after few manipulations, we define Then, with the help of the auxiliary parameter Ξ 2 4 ≡ γ 2 /[ω(1 − ω)] we apply again the Feynman parametrization for the two denominators and rewrite the last equation as In this form, the application of eq. (B.18) in ref. [98] with q = 0 suffices for the carrying out of the momentum integral in L 1 , and eq. (248) becomes where the Gamma functions have been simplified. While the integration over φ is straightforward and gives 2/d, the remaining one can be performed in an analogous fashion as eq. (A.6) of ref. [27].
By considering ≡ 3 − d, in fact, the integrand can be expanded in a power series in , giving Introducing also the renormalization scale µ, the part of the scattering amplitude of interest can be recast as It displays a pole singularity in the limit d → 3 and a simple PDS pole at d = 2. By exploiting the Laurent series expansion of the Gamma function for small arguments and truncating it at NLO, L 1 finally becomes where the regular parts of the amplitude have been evaluated in the three-dimensional limit. Now, due to the presence of a d = 2 singularity, the PDS correction to L 1 is non-zero. Introducing again the renormalization scale µ of the MS scheme and noticing that the integrand in eq. (252) in the d → 2 limit coincides with 1, the correction turns out to be Subtracting the latter equation to eq. (253), the PDS corrected part of the amplitude can be obtained, Next, we concentrate on the evaluation of the second term of eq. (244). Considering the original constant factors and the integral over the q, we define and the constant Ξ 2 5 ≡ γ 2 /[ω(1 − ω)], so that we can exploit again the Feynman parametrization for the denominators, obtaining The integral over the q can be now carried out as a straightforward application of eq. (B.16) in ref. [98], yielding In the last rewriting, the integral over the φ can be immediately performed, while the Gammas can again be simplified and reduced, so that eq. (258) transforms into where the conventional renormalization scale factor (µ/2) 3−d has been introduced as in eq (252).
The remaining integral has been already met in eq. (254) and it can be evaluated exactly in d = 2 or expanded in powers of 3 − d in the three-dimensional case. Plugging eq. (251) into eq. (258), it turns out that the integral is again divergent in the limit d → 3 and includes also a threefold PDS pole in d → 2. Taking the former limit, the Gamma function can be expanded in power series as before and the amplitude can be re-expressed as where the ultraviolet divergence of the original integral is made explicit. Then, the application of the PDS scheme to the d → 2 pole yields the following correction, Finally, we subtract the PDS contribution just determined to eq. (260), obtaining Now we collect the two results in eqs. (255) and (262) and write the one-loop scattering amplitude with one photon exchange in the power divergence subtraction scheme, Here, we focus our attention on the computation of the leading order matrix element of the Coulomb-corrected strong fermion-fermion scattering amplitude in eq. (56). The process can be reduced to the evaluation of only one of the two integrals presented in the first row of the last equation, by virtue of the complex-conjugation property satisfied by the repulsive Coulomb In particular, we choose to concentrate on the first term on the l.h.s. of eq. (56), Recalling the parity rule of the spherical harmonics (cf. sec. A V I in ref. [125]), the repulsive Coulomb eigenstate turns out to be given by Now, we start by observing that the three-dimensional Dirac delta function peaked at the origin can be rewritten in spherical coordinates as due to rotational invariance. The last equation can be proven by exploiting the integral representation of δ(r) [122], then expressing the exponential in terms of spherical harmonics and spherical Bessel functions of the first kind and finally performing the angular and radial integrations by means of the identity in sec. 11.2 of ref. [124]. Besides, the rewriting in eq. (266) paves the way for the integration of the angular variables and on the radial distance on which the Dirac delta effectively acts separately. Equipped with the last equation, we can split the expression in eq. (264) into two parts, where and In the last equation, we note that the application of the gradient effectively reduces to the application of the derivative with respect to the radial variable r , therefore the resulting vector is parallel to r . It is then convenient to exploit the expression of the latter vector in terms of the spherical harmonics given in eq. (5.24) and sec. 5.1 of ref. [108], in order to perform the integration on the radial and the angular variables associated to r in eq. (268) separately. Making use of the complex conjugation (cf. eq. (4.31) in ref. [123]) and the orthonormality (cf. chap. VI of ref. [122]) properties of spherical harmonics, the integral over the angular variables can be carried out rapidly, obtaining where Ω denotes collectively the angular variables associated to r . Recalling eq. (5.24) and sec. 5.1 in ref. [108], the remaining spherical harmonic on the second row of eq. (272) together with the round bracket with the Kronecker deltas can be identified as the unit-vector parallel to p , up to a multiplication factor that cancels out with 2 π/3 on the left of the last integration sign. Then, exploiting eq. (271) for = 1, eq. (272) can be rewritten as (1 + i|p |r)M (2 + iη, 4, −2i|p |r ) − i|p |r 2 + iη 2 M (3 + iη, 5, −2i|p |r ) p .
The explicit evaluation of the limit yields immediately to the disappearance of the terms depending linearly on the radial coordinate r , since the Kummer functions are equal to unity for zero values of the third argument, Besides, by exploiting the fundamental property Γ(z + 1) = zΓ(z) of the Gamma functions, the constants outside the limit in eq. (273) can be rewritten in terms of the Sommerfeld factor (cf. eq. (41)), |Γ(2 + iη)| e πη/2 = 1 + η 2 e πη/2 |Γ(1 + iη)| = 1 + η 2 C η , thus recovering the polynomial on the r.h.s of the generalized effective range expansion (cf. eq. (51)). Equipped with the two last results we can, finally, obtain the desired expression for the r.h.s. of eq. (272), Now, we can proceed with the application of the gradient to the spherical harmonics (cf. eq. (269)).
From ref. [108], the result of the latter derivative can be rewritten as a linear combination of shperical harmonics as in eqs. (5.24) and (5.27) in ref. [108]. In particular, by the introduction of 1 = √ 4πY 0 * 0 (θ, ϕ) in the relevant integral of eq. (269), we can observe that the surface integrals over the Y µ +1 (θ, ϕ) yield no contribution, since cannot assume negative values. It follows that our term of interest becomes where the Clebsch-Gordan coefficients (j 1 j 2 J|m 1 m 2 M ) ≡ JM, j 1 j 2 |j 1 m 1 , j 2 m 2 vanish whenever m = µ + µ . The evaluation of the latter in the second row of eq. (277) leads to Combining the last result together with the one in eq. (276), we obtain the expression of the complete integral in eq. (264), By exploiting the complex-conjugation property of the regular Coulomb repulsive wavefunction , also the second integral in eq. (56) can be evaluated, by noting that the latter property implies only the disappearance of a (−1) factor in eq. (265), that for = 1 is compensated by the overall minus sign in front of the integral, Therefore, the full leading order T SC matrix element in eq. (55) reads where θ is the scattering angle.
addressing the Ϙ 1 calculation. Due to the rapid oscillation of the sum of the series for similar values of the cutoff constant Λ n , the original series in eq. (290) has been recast as where is a small real constant and the exponential factor proves to quench the oscillations of G (0) for neighbouring values of Λ n . Considering the interval 0.1 ≤ ≤ 1, the sum of the series proves to decrease monotonically towards = 0 n the Λ n → +∞ limit and the behaviour is linear with , with small quadratic corrections. The subsequent quadratic interpolation, in fact, returns a value of G (0) compatible with the exact one in literature (cf. eq. (290) and ref. [40]), It follows that the chosen approach (cf. eq. (290)) is successful in the evaluation of the sum of the series and can be promoted to more involved cases. Moreover, exploring a larger interval of ε towards larger values, further deviations form linearity are likely to appear in the fit. In particular, the following class of fitting functions, is expected to provide a satisfactory description of the behaviour of of G (0) with ε both in the vicinity of zero and in the infinite ε limit.
Correlated to G (0) is the I (0) series, which appears in both in the QED leading order corrections to the scattering parameters (cf. eqs. (134)-(139)) and in the large binding momentum limit of the Lüscher functions S 2 (iκ) and S 3 (iκ), see eqs. (222) and (223). The sum of this series is already known in literature [40,64,88] and is given by where the first equality is shown in tab. 1 and eq. (2.61) of ref. [64]. A precise evaluation of I (0) can be attained by isolating the cutoff-dependent part of the series via the Poisson summation formula. In particular, the addition and subtraction of a 1/(|n| 2 + 1) term in the original series yields where the first term on the r.h.s proves to converge in the Λ → +∞ limit as fast as the J ≡ J (0) series in ref. [88] and the linear divergence is confined into the second summation. Exploiting the Poisson formula in eq. (285), the non-regularized series can be evaluated as follows, where the zero n term has been added to the sum and, then, the spherical cutoff has been moved from the sum to the integral over the lattice momenta. Separating the zero modes from the others in the result of eq. (294), we obtain where the third and the fourth term on the r.h.s. turn out to be finite in the infinite cutoff limit and the last integral can be performed in the complex plane, by collecting the residues according to Jordan's Lemma. Performing the remaining integration, in fact, the last equation becomes an expression that can be directly plugged into eq. (294), obtaining the desired result, As it can be observed, the cutoff dependent term in the original series has been removed and, at the same time, two rapidly convergent sums, Ϡ 0 = 0.0125 and γ 1 (cf. eq. (302)) replaced the divergent one, thus reducing significantly the computational efforts. The procedure is completely analogous to the one adopted in Appendix B 1 of ref. [61] and can be applied to other cutoff-regulated single sums.
At this stage, we focus on the single sums appearing as coefficients in the δp 2 expansion of the Lüscher functions S 1 (p), S 2 (p) and S 3 (p) aroundp 2 = 1. Adopting a notation similar to the one used for eq. (290) and (293), the sums of relevant three-dimensional Riemann series yield where the first three coincide respectively with I 1 , J 1 and K 1 in Appendix C of ref. [61]. In particular, all the series not regulated by a cutoff can be computed directly, without the need to resort to the techniques outlined above. On the other hand, the sum of I (1) can be obtained rapidly from the existing result for I (0) . In fact, the addition and subtraction of a 1/|n| 2 term gives where in the first term on the r.h.s. the limit Λ n → +∞ has been taken. Moreover, replacing I (0) with its expression given in eq. (293), the last formula can be rewritten in a compact fashion as Both the classes of series in eqs. (300) and (301) do not display convergence issues and can be directly evaluated. Note that ϝ 1 coincides with the sum listed as χ 3 in eq. (A1) of ref. [88].

D.2 Double Sums
Differently from their single counterparts, double sums appear only in the purely Coulombic contributions in the = 1 ERE and arise from thep 2 → 1 limit of the Lüscher functions S 2 (p) and The presence of the factor 1/|m − n| 2 in the sums χ 1 − χ 4 of eq. (304) ensures convergence without the need for the introduction of spherical cutoffs and regulators. Besides, the sum χ 0 is analytical and denoted as χ 1 in ref. [88], whereas the series χ 2 and χ 3 coincide with the χ 2 and χ 5 , respectively, in the latter work, except for the inclusion of the zero mode in the sum over m.
Another group of series belonging to the same category is provided by the sums which are not present in the expansions of the S 2 (p) and Now we switch to the second category of double sums, the one consisting of two threedimensional sums performed on Z 3 . First, we consider the series stemming from thep ∼ 1 approximations of S 2 (p), see eq. (193). The latter sums, in fact, are the counterpart of divergent double integrals contributing to the amplitudes of the relevant two-particle scattering processes.
Due to the large increase of the configuration space, for the numerical calculation of the sum of such series it is advisable to parallelize the operations via the development of GPU codes (e.g. in Cuda C++). The computational efforts can be significantly reduced by subdividing the original double sum into an arbitrarily large finite number of single sums, characterized by a three dimensional vector of integers. Then, assigning each of the outcoming single sums to a different subunit of a graphic card, the sum of the original double series is derived by gathering the results obtained simultaneously by each operating unit.
Once in this form, the sum of R (1) can be obtained by exploiting the existing result for the cutoffregularized sum R (0) in ref. [88], together with the single sums in eqs. (303), (304) and (305).
The only additional computational effort is given by the double sum explicitly shown on the r.h.s of eq. (307), which proves to converge rapidly, differently from R (0) . For the last sum, in fact, an approach analogous to the one shown in eq. (290) or to the tail-singularity separation (TSS) in ref. [126] is recommendable. In the latter method, summarized in detail in ref. [88], a threedimensional Riemann sum is subdivided into an IR part, dominated by the the singularities of the summand and an UV part, expressed in the form of a three-dimensional integral and describing the behaviour of the argument of the original sum towards the infinity. As shown in the Appendix A of ref. [126] for the sums θ As and θ Bs , the TSS approach holds also for double sums regulated asymmetrically.
Second, we switch to the series of the kind R As noticed in ref. [88], the only contribution of such sums in the expression of the finite volume energy corrections for the lowest energy T 1 eigenstate (cf. eq. (218) and (219)) is provided by 42 , whose explicit evaluation gives −93.692 (10). Finally, we present the double sums arising from thep ∼ 1 approximation of S 3 (p) in eq. (209), starting with the cutoff-regularized double sum Ϙ 1 in eq. (204). Since the argument of the sum is quadratically divergent with the spherical cutoff, it is convenient to adopt a stabilization technique for the evaluation of the sum. To this purpose, we chose to apply the approach in eq. (290) to the cutoff-regularized sum over n, Since the ε-dependent Ϙ 1 sums display a non-linear behaviour in the interval 0 ≤ ε ≤ 1.1, we choose to interpolate the data with the fitting function in eq. (292). As shown in fig. 7, f (ε) describes the behaviour of the sum of the stabilized series as a function of ε satisfactorily, therefore the sum of the series becomes Conversely, the double sum Ϙ 2 (cf. eq. (205)) appearing at order δp 2 in the power series expansion of S 3 (p) can be calculated efficiently even without stabilization approaches, despite its sign-changing numerator. Its numerical evaluation yields 1 − |m| 2 |n| 2 + n · m(|n| 2 + |m| 2 − 2) (|m| 2 − 1) 2 (|n| 2 − 1) 2 |m − n| 2 = −315.981 (74) .

Appendix E Lüscher functions
In this appendix, we concentrate on the derivation of the large imaginary momentump = iκ where Λ m is a spherical cutoff accounting for the divergence of the following sums over m, Due to the presence of a cutoff in the inner sum, the translation in momentum space operated in the S 2 (iκ) case (cf. Appendix A3 in ref. [88]) is no longer allowed in the individual sums (cf. eqs. (316)-(318)). Nevertheless, since the purpose is the extraction of the finite andκ-dependent contributions from each of the three double sums in eq. (315), terms depending on nonzero powers of Λ m and Λ n can be neglected without loss of information for the FVECs. Therefore, we assume henceforth the limits Λ n , Λ m → +∞ and extract the finite parts from the sums depending on the binding momentum. Now, we consider the first of the three double sums in the second row of eq. (315). After the translation in the momenta m → m − n ≡ p, we rewrite the argument of the inner sum in integral form, Then, we apply the Poisson summation formula to the unconstrained sum over n and isolate the zero modes from the nonzero ones, e −i2πq·n |q + p| 2 +κ 2 . (320) As a consequence the second term on the r.h.s. of eq. (330) can be finally rewritten as where the last two terms on the r.h.s. are exponentially suppressed and they have to be neglected for consistency. Second, we switch to the single sum on the r.h.s. of eq. (330). Introducing the integral sign and exploiting again the Poisson summation formula, the last contribution to S 3 (iκ) where the zero modes have been again isolated from the non-zero ones. The first integral in eq. (334) can be carried out after few manipulations, Then, the integration over the angular variables associated to q in eq. (335) gives where q ≡ |q|. As in the case of eq. (324), the integrand is an even function of q, thus the integration region can be extended to the whole real axis. Besides, the integrand tends uniformly to zero in the limit l → ±∞ and is analytical all over the complex plane l ∈ C, except for two double poles at q ± = ±iκ, located along the imaginary axis. It follows that the integration region can be extended to an arbitrary large circular region about the origin of the complex plane encompassing the two singularities. Moreover, the integrand can be split into two functions of complex variable q h ± (q) = q i|n| e ±i2πq|n| (q − q + ) 2 (q − q − ) 2 , so that h + (q) (h − (q)) can be integrated in a semicircumference with arbitrarily large radius about the origin in the upper (lower) part of the complex plane picking up the q + (q − ) singularity, according to Jordan's Lemma. Again, the residues about the two double poles turn out to coincide and to depend on |n| through negative exponentials, Now, observing again that e −2π|n|κ ≤ e −πκ e −π|n| , the sum in eq. (338) can be bound from above where Ϡ 2 is a small constant equal to 0.485647 (1). It follows that the contribution of the nonzero modes associated to the single sum on the r.h.s. of eq. (330) is exponentially suppressed and can be neglected in the large binding momentum limit. Collecting all the results in eqs. (328), (329), (333) and (335), the large binding momentum limit of the double sum S 3 (iκ) is found, where the ellipses include the cutoff-dependent divergent terms and functions ofκ which are suppressed by negative exponentials.