Two-particle scattering from finite-volume quantization conditions using the plane wave basis

We propose an alternative approach to L\"uscher's formula for extracting two-body scattering phase shifts from finite volume spectra with no reliance on the partial wave expansion. We use an effective-field-theory-based Hamiltonian method in the plane wave basis and decompose the corresponding matrix elements of operators into irreducible representations of the relevant point groups. The proposed approach allows one to benefit from the knowledge of the long-range interaction and avoids complications from partial wave mixing in a finite volume. We consider spin-singlet channels in the two-nucleon system and pion-pion scattering in the $\rho$-meson channel in the rest and moving frames to illustrate the method for non-relativistic and relativistic systems, respectively. For the two-nucleon system, the long-range interaction due to the one-pion exchange is found to make the single-channel L\"uscher formula unreliable at the physical pion mass. For S-wave dominated states, the single-channel L\"uscher method suffers from significant finite-volume artifacts for a $L=3$~fm box, but it works well for boxes with $L>5$~fm. However, for P-wave dominated states, significant partial wave mixing effects prevent the application of the single-channel L\"uscher formula regardless of the box size (except for the near-threshold region). Using a toy model to generate synthetic data for finite-volume energies, we show that our effective-field-theory-based approach in the plane wave basis is capable of a reliable extraction of the phase shifts. For pion-pion scattering, we employ a phenomenological model to fit lattice QCD results at the physical pion mass. The extracted P-wave phase shifts are found to be in a good agreement with the experimental results.

The Lüscher method assumes that the box size L is much larger than the interaction range R, such that exponentially suppressed corrections ∼ e −L/R are negligible. In currently feasible lattice simulations, the box size can usually not be chosen very large given the available computational resources. On the other hand, in certain cases such as e.g. the two-nucleon system, the long-range interaction is known and appears to play a very prominent role. Another example of a system featuring long-range dynamics is theD * D/DD * system, in which the one-pion exchange (OPE) potential plays an important role for the formation of exotic hadrons X(3872) and Z c (3900) [18][19][20]. The accidental mass relation mD * ≈ mD + M π makes the exchanged pion in the OPE potential almost on-shell, which results in the interaction range much larger than 1/M π . In such situations, taking into account the long-range dynamics may allow one to reduce the cost of lattice QCD simulations by considering smaller volumes. It, however, requires the finite-volume Lüscher quantization conditions to be generalized by taking into account exponentially suppressed effects due to the long-range interaction.
In single-channel cases, each energy level in a finite volume (FV) can be related to the scattering phase shift at this energy using interaction-independent Lüscher's formula. On the other hand, one has to deal with coupled-channel effects when taking into account inelastic channels or partial wave mixing. In the multichannel case, Lüscher's formula becomes a determinant equation det[f (T, E)] = 0, which relates the T -matrix and the FV energy levels. Though this result is still model-independent, there is no one-to-one correspondences between phase shifts and FV energy levels. It does not allow one to extract the full T -matrix, i.e. phase shifts for every channel and the corresponding mixing angles. In practice, the T -matrix can be parameterized within a particular theoretical framework, and several FV energies can be used as input to determine the parameters in the T -matrix. This may, however, introduce some model dependence. Another complication is related to rotational symmetry breaking in finite boxes, which results in the energy levels generally receiving contributions from multiple partial waves.
In the literature, much effort has been devoted to cure or estimate the two drawbacks of the Lüscher formula, namely exponentially suppressed corrections and mixing effects in multichannel scattering. In refs. [21,22], the authors estimated exponentially suppressed corrections for ππ and NN scatterings in the FV. It was shown that for realistic pion masses, simulations in a box of the size L > 5 fm would result in exponential corrections to the phase shifts less than one degree. In ref. [23], finite volume corrections to the X(3872) were analyzed in an effective field theory (EFT) with perturbative pions. It was shown that FV effects are significant for box lengths as large as 20 fm. Unitarized chiral perturbation theory in a finite volume was employed in refs. [14,[24][25][26], in which some exponentially suppressed effects were included by adopting a fully relativistic propagator. Such exponentially suppressed differences with the Lüscher formula for ππ scattering were calculated explicitly in refs. [27,28]. In order to deal with coupled-channel effects in a FV, an effective Hamiltonian approach was developed in refs. [29][30][31][32]. Within different approaches, partial wave mixing effects were estimated in refs. [7,14,33,34]. The sensitivity to the second lowest partial wave in Lüscher's quantization conditions was discussed in ref. [35].
Starting from Lüscher's seminal work, FV effects were usually investigated using the partial wave expansion technique. In the infinite volume, partial wave expansion allows one to simplify the scattering problem thanks to the rotational symmetry. However, in the finite volume, continuous rotational symmetry is broken and the angular momentum is no longer a good quantum number. To some extent, introducing the partial wave expansion in a finite volume complicates the problem. For small boxes and/or systems with coupled channels as well as in the presence of long-range interactions, the above mentioned drawbacks of the Lüscher approach get amplified. A natural alternative method is to employ vectors of discrete momenta. In refs. [7,34], the authors used the discrete momentum basis to investigate partial wave mixing effects. Recently, a Fourier basis {cos(2πn i r i /L), sin(2πn i r i /L)} in coordinate space was adopted to diagonalize the Hamiltonian in FV [36].
In this work, we employ a similar approach as that of ref. [36] but in momentum space. We consider the Lippmann-Schwinger equation (LSE) in the non-relativistic case or a reduced Bethe-Salpeter equation (BSE) in the relativistic case using the plane wave basis with discrete momenta. Next, we reduce the matrix LSE and BSE into the irreducible representations (irreps) of the cubic group. With the reduced matrix equation, we compute the energy levels in the FV corresponding to specific irreps. In such a framework, effects from mixing with higher partial waves are embedded naturally. We also extend this approach to systems with nonvanishing total momentum. We consider NN and ππ scattering as examples of the non-relativistic and relativistic systems, respectively, to illustrate the calculations and compare our results with those obtained using the single-channel Lüscher approach. In the case of the NN system, we also study the role played by the long-range interaction due to the one-pion exchange.
This work is organized as follows. In section 2, we review the discretization conditions of three-momenta and discuss the relevant symmetry groups for two particles in a cubic box with periodic boundary conditions. In section 3, we construct the representation spaces of the corresponding point groups with the plane wave basis and perform their reduction into irreps using the projection operator technique. In section 4, we introduce a general approach to obtain the FV energy levels and to fit lattice QCD energy spectra in the FV. Next, in section 5, we consider spin-singlet NN scattering as an explicit example of a non-relativistic system to demonstrate our method. In this section, we also discuss in detail partial wave mixing effects in the FV. In section 6, we consider ππ scattering in the ρ-meson channel as an example of a relativistic system. Using a phenomenological model to parameterize the ππ scattering amplitude, we calculate the FV energy levels and fit the parameters of the model to the lattice QCD spectra. Finally, in section 7, we summarize the main results of our study and discuss possible generalizations. Some basic information about the point groups used in this work and Lüscher's quantization conditions is given in appendices A and B.

Two particles in a finite volume
Discretization of three-momenta of two particles in a cubic box with periodic boundary conditions is considered e.g. in refs. [4,9]. In the following, we briefly review the relevant findings for the sake of completeness. Throughout this paper, we focus on systems of two spinless particles of the same mass. It is straightforward to generalize our results to particles with arbitrary spin and different masses.

Non-interacting case
In the box frame (BF), the momenta p 1 and p 2 of two non-interacting particles are discrete, where P and E are the total momentum and energy of the two-particle system, respectively. For the non-interacting case, each particle is on-shell with E i = » p 2 i + m 2 i . In the center-of-mass frame (CMF), the related momenta and energies read where the quantities in the CMF are labeled by superscript " * ". In eq. (2.3), the on-shell relation of each particle is not used. Thus, eq. (2.3) is also valid for interacting systems. For non-interacting systems, the momentum and energy in the CMF are related to those in the BF by the Lorentz transformation where γ = E E * is the Lorentz factor. Equations (2.4) and (2.5) can be combined and simplified to For a given vector u, u and u ⊥ are defined via u = (u·P )P P 2 and u ⊥ = u − u . In this paper we are interested in systems of two particles with equal masses (and thus A = 1). With eqs. (2.2) and (2.6), one obtains the discretization relations in the CMF In figure 1, we show a set of discrete momenta with the corresponding meshes in different frames. The first mesh in each row corresponds to the momenta in the BF. The second and third ones illustrate the steps needed to obtain the momenta in the CMF according to eq. (2.8), namely the shift of the origin by d/2 and the deformation in the direction of d. In figure 1, we take d = (0, 0, 1) and d = (1, 1, 0) as examples. The meshes show the symmetries of the system in the FV. In the infinite volume, the momenta are continuous and satisfy the rotational symmetry  Figure 1. Meshes of the set of discrete momenta for a two-particle system with equal masses in the BF (left column) and CMF (right column). The second and third meshes in each row illustrate the shift of the origin by d/2 and the deformation in the direction of d to obtain the mesh in the CMF.  figure 1 is not in the center of its adjacent points anymore after shifting by 1 2 Ad, and the system is not invariant with respect to space inversions [9]. In appendix A, we list the group elements, conjugacy classes, and character tables for the relevant point groups.
In the literature, two conventions to denote irreps of D 4h , D 2h and D 3d are used. For example,  Table 2. Irreps of the relevant point groups using the two conventions as discussed in the text and their corresponding relations. The convention I is adopted in ref. [4,37] and in this paper. The convention II is employed in refs. [6,38].
the D 4h group can be formed either as The elements of D 4 are all proper rotations. The C 4v group contains two improper rotations and is the symmetry group for systems with particles of unequal masses. In the first convention, the irreps of D 4h are denoted by the irreps of the D 4 group with parity [4,37]. In the second convention, the irreps of D 4h are denoted by the irreps of the C 4v group with parity [6,38]. We list the irreps of the two conventions and their corresponding relations in table 2. In this work, we adopt the first convention.
For non-relativistic systems, the above discussion can be simplified since the Lorentz factor becomes γ = 1 and the Lorentz transformation reduces to the Galilean one. Thus, the second step in figure 1 is unnecessary. One consequence is that the meshes for non-relativistic systems might be more symmetric than the ones with the same d in a relativistic theory. For example, in a non-relativistic theory, the mesh for the d = (0, 0, 2) system is the same as that of the d = (0, 0, 0) system. The shift of origin by d/2 = (0, 0, 1) in figure 1 makes no difference for momentum discretization. Thus, the discrete symmetry group of a non-relativistic system with d = (0, 0, 2) is O h . However, the relativistic system with the same d only possesses the D 4h symmetry. In this work, for calculational convenience, we adopt an unified approach for both non-relativistic and relativistic systems. E.g., for a given system with d = (0, 0, 2), we classify the energy levels according to the irreps of D 4h . Thus, for non-relativistic systems, we will end up with several degenerate states in different irreps reflecting the additional degeneracy of the more symmetric group O h .

Interacting case
To compute FV energy spectra in the interacting case, we need to perform boost transformations for loop momenta with particles being off shell, i.e. E In the literature, different schemes have been proposed to define the quantization conditions for the cases with P = 0. The one of ref. [14] starts with the prescription which is the same as in the non-interacting case and ensures that E * 1 + E * 2 = E * , and E * 1 = E * 2 = E * 2 for m 1 = m 2 . By substituting eq. (2.9) into eqs. (2.4) and (2.5), one obtains the same quantization conditions as those in eqs. (2.6) and (2.8).
Recently, an alternative scheme was proposed in ref. [32]. Their relation of p * 1 with p 1 reads where ω i = » m 2 1 + p 2 i . The Lorentz factor γ is defined as In the first scheme, the Lorentz factor γ is defined rigorously, and it only depends on the total energies of the two-particle system in the BF and CMF. In the second scheme, the Lorentz factor is taken as an approximation using the expression for the non-interacting case. Such an approximation makes the transformation independent of the energy, which simplifies the determination of poles of the T -matrix. It was argued that the difference between the two schemes is exponentially suppressed [32].
In the infinite volume, the LSE and BSE are usually written in the CMF. We need to replace the integration over the loop momentum in the LSE or BSE with a summation over the discrete momentum values. For a general case with P = 0, the transformation reads where J is the Jacobi determinant for the momentum transformation from the CMF to LF. The explicit expressions for J are different for the two schemes in eqs. (2.6) and (2.10), The ambiguities in choosing the transformations in the interacting case are related to the need to define the interactions in the off-shell kinematics. However, for non-relativistic systems, the two schemes reduce to the same one, and the Jacobi determinant in eq. (2.12) becomes J = 1.
3 Construction of irreducible representations using the plane wave basis

Representation space
The rotation operator in the Hilbert space isD(g(n, θ)) = e −in·Ĵθ , whereĴ is the angular momentum operator and g(n, θ) is the rotation around n axis with angle θ. For spinless systems in the infinite volume featuring the SO(3) symmetry, the space spanned by |lm , m = −l, ..., l, forms an irrep. The matrix elements of the rotation operator in the irrep give rise to the Wigner D-function defined via l m |D(g)|lm ≡ D l m m (g)δ ll , (3.1)  where l and m are the quantum numbers of angular momentum and its third component, respectively. In the finite volume, the Hilbert space with fixed l also forms a representation space of O, D 4 , D 2 , and D 3 groups, but it becomes reducible. The procedure of reducing the Hilbert space with a fixed l to the irrep ones of the corresponding point group is well established, see e.g. ref. [39].
In the infinite volume, the Hilbert space spanned by the plane wave basis |p also forms a representation space of the O(3) group, In the finite volume, a natural choice for the representation space of the cubic O h group is {n 1 , n 2 , n 3 } ≡ {|n 1 , n 2 , n 3 with permutations of n 1 , n 2 , n 3 and changing signs}, where n i is a multiple of 2π L for each component of the momentum. We use the notation |n i , n j , n k to denote the plane wave basis in a finite cubic box. It is obvious that n i are integer numbers for the O h group when d = (0, 0, 0). The dimensions of {n 1 , n 2 , n 3 }, which we refer to as patterns throughout this paper, depend on the number of zeros in three components and the number of equality relations among them. In table 3, we list seven different patterns and their dimensions. The group representations in the {n 1 , n 2 , n 3 } space are reducible. Fortunately, the reduction only depends on the patterns of this space. For the O h group, we only need to deal with the reduction for seven patterns as shown in table 3. Notice that the A i and B i irreps are one-dimensional ones, while the E and T i representations are two-and three-dimensional irreps, respectively.
For moving systems in a box with d = (a, a, a), the corresponding point group is D 3d and we can choose a similar representation space as in eq. (3.3). In this case, n i are unlikely to be integer numbers according to the quantization condition in eq. (2.8). In order to obtain the irreps, we need to deal with at most seven different patterns shown in figure 3. If a is odd, the patterns with zero components will not appear.
For moving systems with d = (0, 0, a) or d = (a, a, 0), it is more convenient to employ the representation spaces {n 1 , n 2 ; n 3 } ≡ {|n 1 , n 2 , n 3 with permutations of n 1 and n 2 and changing signs}.
(3.4)   .4) we can obtain the representation matrices for each pattern using eq. (3.2). Next, we are going to reduce the representations into the direct sums of irreps.

Reduction of the representations
The procedure of reducing the representations with the help of projection operators is well established [40]. It has been used to reduce representations in the partial wave basis [39] and in the Fourier basis [36]. For a related discussion see ref. [41]. We will employ the same technique to reduce representations in the plane wave basis. The projection operator is defined aŝ whereD(g i ) is the (im)proper rotation operator in the Hilbert space corresponding to the group element g i as discussed in section 3.1, Γ a denotes the corresponding irrep, while N (Γ a ) and n G refer to the dimensionality of the irrep and the order of the group, respectively. Further, R Γa is the unitary matrix irrep of Γ a . For the point groups relevant for this study, the unitary irreps can be taken from the literature [6,38,39]. We will specify the procedure of constructing R Γa later. Here, we list some properties of the projection operator without proof [40], When acting on a general state |ψ defined as the projection operator projects the state onto a single component of a fixed irrep, For the case at hand, we choose |ψ as the basis of the representation space |n 1 , n 2 , n 3 , fix α and vary |α to obtain the basis states |Γ a , α . With the transformation from |n 1 , n 2 , n 2 to |Γ a , α , we reduce the original representation to a direct sum of irreps. The unitary matrix irreps R Γa in eq. (3.5) can be constructed with the character projection operator, which is defined aŝ Acting with the character projection operator on |ψ , we obtainP One can vary |ψ to obtain a sufficient number of vectors, which span the irrep space after orthogonalization. In particular, in order to obtain the R Γa , one can start with the faithful representation. For point groups, |ψ andD(g i ) are replaced with the vector (x, y, z) and (im)proper rotation matrices in the Euclidean space, respectively. For point groups, the character table is unique but the unitary irrep matrices are not. One can choose a convention to make R Γa real symmetric matrices [6,39].
In summary, the general procedure to reduce a representation to irreps is as follows: 1. Identify the symmetric group and its elements and character table,

Construct the unitary irrep matrices with character projection operation,
3. Reduce the representation to a direct sum of irreps.

Finite volume energy levels: determination and fitting
In general, FV energy levels can be obtained by finding the poles of the T -matrix in the box, calculated as described in sections 5 and 6. Using the technique described in section 3, one ends up with a determinant equation for a specific irrep Γ The specific expressions for the energy-dependent matrix M will be presented for the non-relativistic and relativistic systems in the next two sections. In general, one can adopt the root-finding algorithm to search for the roots of eq. (4.1). To this end we define Ω [33] where λ Γ,i denotes the eigenvalue of the matrix M Γ . The root of eq. (4.1) corresponds to at least one λ Γ,i = 0. Thus, the zeros of Ω Γ (E) are the roots of eq. (4.1). The nonzero µ can be chosen to optimize the root-finding procedure. The quantity Ω Γ (E; µ) is bounded between −1 and 1.
Introducing Ω Γ (E; µ) thus allows one to avoid the numerical complication when finding the roots with the determinant in eq. (4.1) becoming very large for large-dimensional matrices.
Meanwhile, one could use eq. (4.1) to determine the parameters of the interaction within a given theoretical framework by matching to the FV energy levels from lattice QCD. To this end, we have two options, the spectrum method and the determinant residual method [33]. In the spectrum method, one compares the FV energy levels from the calculation with the lattice QCD spectrum and minimizes the residuals. In the fitting procedure, the FV energy levels need to be obtained by the root-finding algorithm repeatedly, which is time-consuming, see also ref. [42] for a discussion of further challenges. Alternatively, in the determinant residual method, the determinant det(M Γ ) that becomes zero when the lattice QCD energy level is the solution of eq. (4.1), is regarded as a residual. As already pointed out above, it is more practical to employ the quantity Ω Γ (E; µ) defined in eq. (4.2) as a residual. Then, the χ 2 -function to be minimized is defined via is the error of Ω Γ propagated from the error of E Γ,i . In this exploratory study, we neglect possible correlations among the different FV energy levels.
Here and in what follows, we employ the determinant residual approach.

Application I: Spin-singlet two-nucleon scattering
We now consider nucleon-nucleon scattering in spin-singlet channels as an example of a nonrelativistic system. This is a particularly interesting case for several reasons. First, the formalism described in sections 2 and 3 for spinless particles of equal masses is well suited for NN scattering in spin-zero channels assuming the exact isospin symmetry. Secondly, because of the appearance of the strong long-range interaction due to the OPE, the NN system is expected to feature severe FV artifacts due to partial wave mixing, which is one of the central focuses of our study. Last but not least, the quantization conditions for momenta in section 2 take a particularly simple form for non-relativistic systems with energy-independent potential and allow one to compute the FV energies by simply solving the eigenvalue problem.
In the past decades, nuclear forces have been extensively studied in the framework of chiral EFT. See refs. [43][44][45] for review articles. State-of-the-art NN interactions at fifth order (N 4 LO) in chiral EFT have been shown to provide an excellent description of the available neutron-proton and proton-proton scattering data up to pion production threshold [46][47][48]. For the purpose of this proof-of-principle study, we restrict ourselves to next-to-next-to-leading order (NNLO) and employ the nonlocally regularized potentials of ref. [49]. The NNLO potential reads cont are the leading-order (LO) and next-to-leading order (NLO) contact interactions, respectively. V (0) 1π and V (2) 1π denote the LO one-pion exchange interaction and its correction at NLO. The NLO correction to the OPE is included by taking a larger value of the (effective) nucleon axial vector coupling g A = 1.29 to account for the Goldberger-Treiman discrepancy [50,51]. The  [53]. Solid lines correspond to the NNLO chiral EFT results [49,52].
two-pion-exchange interactions start contributing at NLO. For spin-singlet NN channels, we can replace the spin operators in the potentials in Ref. [49,52] as where q = p − p. For spin-singlet channels, it simplifies to At NLO, the contact interactions read where k = (p + p )/2, while C S , C 1 and C 2 are low energy constants (LECs). Here, the general expression for the contact interactions including spin operators and consisting of 9 independent terms has been reduced using eq. (5.2). The LECs C S , C 1 and C 2 are thus linear combinations of the original LECs appearing in the chiral NN potential. We take the specific values of the LECs in Table 4 of ref. [49] and employ the same numerical values for all the other parameters. The expressions for the two-pion-exchange interactions can be found e.g. in ref. [49]. Following the standard procedure, we introduce the regulator by the replacement, where the cutoff is taken as Λ = 0.55 GeV.
To calculate phase shifts in the infinite volume, we solve the LSE in the CMF where E is the energy of two nucleons and m N is the nucleon mass. Clearly, for infinite-volume calculations, it is more convenient to use the standard partial wave basis instead of solving the equation in real space. In figure 2, we compare the resulting neutron-proton phase shifts at NNLO with the ones from the Nijmegen partial wave analysis (PWA) [53]. This comparison clearly demonstrates that the NNLO potential already captures the most important features of the nuclear force. Below, we will regard the chiral EFT potential introduced above as the underlying NN interaction to generate the corresponding FV energy spectra. We then compare our approach to reconstruct phase shifts from the FV energies with the single-channel Lüscher formula and explore FV effects due to partial wave mixing.

Lippmann-Schwinger equation in a finite volume
To solve the LSE in a finite volume, we employ the plane wave basis with quantized momenta instead of expanding it in partial waves. In the non-relativistic limit, eq. (2.8) reduces to We define the matrices T and G according to where q n = 2πn/L. Similarly, we define V n ,n . The integration has been replaced with the summation according to eq. (2.12). We introduce a truncation of the discretized momentum q n < q max . For a non-relativistic theory, the Jacobi determinant is J = 1 and the LSE turns to a matrix equation The poles of the T -matrix are obtained by solving equation If the potential is energy-dependent, the root-finding algorithm should be adopted. For energyindependent potentials, eq. (5.11) becomes an eigenvalue problem, where I is the identity matrix. We see that the equation becomes a Hamiltonian equation. We, however, want to stress that this Hamiltonian equation is different from the one emerging in the approach of refs. [29][30][31][32]. In that case, partial wave expansion was made prior to discretizing the magnitudes of the momenta. In our calculation, we discretize the momenta as three-vectors with no reliance on the partial wave expansion. For a related approach see ref. [54].
Using the projection operator technique described in section 3, we reduce the matrices in eq. (5.10) to block diagonal ones. Therefore, equation (5.12) turns into a set of matrix equations according to different irreps. The FV energy levels for a given irrep Γ are then obtained by solving the eigenvalue problem for this irrep where H Γ belongs to the block of H corresponding to the irrep Γ.

The single-channel Lüscher method
We consider boxes of three different lengths: L = 3, 5, 8 fm. In our calculations, we truncate the momentum at leading to q max ≈ 4.1, 2.5, 1.5 GeV for L = 3, 5, 8 fm, respectively. These momenta are much larger than the cutoff parameter in the employed NN potential so that the solution of the LSE is well converged. We have verified that the error due to the truncation is negligibly small.

Benchmark calculation using contact interactions
To confirm the robustness of our approach, we first calculate the FV energy levels in pionless EFT. Specifically, we include the LO and NLO contact interactions, which contribute to the 1 S 0 and 1 P 1 channels. The coefficients in front of the contact interactions have been chosen the same as those in chiral EFT to the NNLO in ref. [49]. The regulators and cutoff parameters are introduced in eq. (5.6). The corresponding phase shifts, calculated by solving the LSE in the infinite volume and shown by solid lines in figure 3, feature a qualitatively similar behavior to the empirical ones in NN scattering. We then compute the FV energy levels by solving the LSE in the discretized plane wave basis as described in sections 3, 4. We regard the calculated energy spectra as synthetic data and apply the single-channel Lüscher formula to extract the phase shifts at these energies. The Lüscher quantization conditions are outlined in appendix B. The full quantization conditions are determinant equations involving all partial waves. In this study, we will refer to the approximation involving only a leading partial wave, which sets one-to-one relations between FV energy levels and phase shifts, as the single-channel Lüscher approach. In figure 3, the extracted 1 S 0 and 1 P 1 phase shifts in the boxes with L = 3, 5, 8 fm are compared with the infinite-volume results. For the positive parity states, we only list the A + 1 states, since the other irreps have no 1 S 0 components. With increased box sizes, the unit momentum (2π/L) becomes smaller leading to denser energy spectra. In the infinite-volume limit, these discrete energy levels turn to the continuum spectrum. Meanwhile, one can see many energy levels corresponding to vanishing phase shifts (i.e. to the non-interacting case). In the FV, partial wave states with the same parity mix with each other due the violation of rotational symmetry. For example, the 1 S 0 states mix with the 1 D 2 , 1 G 4 , etc., ones. Thus, in general, energy levels with positive parity receive contributions from the interaction in the 1 S 0 , 1 D 2 and even higher partial waves. Contact interactions at NLO, regularized with an angle-independent regulator, do, however, not contribute to D-and higher partial waves. Thus, the states in figure 3 with vanishing phase shifts correspond to the non-interacting 1 D 2 , 1 F 3 and higher partial wave states. Apart from these points, the phase shifts calculated using the single-channel Lüscher formula are very close to the ones obtained in the infinite volume calculation. As expected, the single-channel Lüscher approach is found to work accurately for this case without partial wave mixing. Notice that since for the case at hand the range of the interaction is given by the inverse cutoff, the exponentially suppressed corrections to the Lüscher formula are negligible even for the smallest considered box size.

Chiral EFT at NNLO
We now turn to the realistic case and repeat the analysis outlined in the previous section for the chiral EFT NN interaction at NNLO. The results of our FV calculations are compared with those in the infinite volume in figure 4. Notice that for positive parity states, we only show the results using the A + 1 irrep. We use the single-channel Lüscher formula suitable for the 1 S 0 partial wave in the case of the positive parity states and for the 1 P 1 partial wave in the case of the negative parity states. Obviously, for the realistic NN interaction at the physical pion mass, the single-channel Lüscher formula leads to significant deviations from the infinite volume calculations.
For positive-parity states shown in the upper row of figure 4, the deviation of single-channel Lüscher's results from the exact ones is significant for the smallest considered box with L = 3 fm. When the box size is increased to 5 fm, the Lüscher single-channel formula performs reasonably well. Apart from the S-wave dominated states, one can identify several D-wave dominated states in figure 4 leading to smaller-in-magnitude phase shifts. Though we employ the Lüscher formula for the S-wave, its leading part w 00 is the same as that for the D-wave, see appendix B. Therefore, we obtain in figure 4 an approximation of the D-wave phase shift from D-wave-dominated states using the S-wave Lüscher formula.
The states with negative parity in the lower row of figure 4 are more interesting. Even for the largest considered box with L = 8 fm, one observes large deviations from the infinite-volume results when using the single-channel Lüscher formula. One identifies two regions where the points are concentrated, namely well above the solid line (region-I) and below the solid line (region-II). One may expect the points in the region-I to emerge from states dominated by the F-wave (and possibly by even higher partial wave) interactions. The points in the region-II are found to gradually  deviate from the exact result with increasing energies. Only in the very near-threshold region, the phase shifts for several energy levels are close to those from the infinite volume calculation. Notice that further increasing the box size 1 is found to yield no improvement for the negative-parity phase shifts.

The one-pion-exchange potential
The results of the previous two sections suggest that the failure of the single-channel Lüscher formula for the NNLO chiral EFT potential is caused by partial wave mixing effects. For interactions with a finite non-zero range R, the near-threshold behavior of the phase shifts δ l (p on ) with E = p 2 on /m N is given by δ l (p on ) ∼ p 2l on . This asymptotic expression, however, only applies to on-shell momenta well below the lowest t-channel singularity, i.e. for p on < R −1 . For the NN interaction, this restriction translates into E lab ∼ 2M 2 π /m N ∼ 10 MeV. At such low energies, the single-channel Lüscher formula is expected to be valid since the contributions from higher partial waves are kinematically suppressed. This expectation is in line with the findings of the previous section for the largest considered box size. For p on R −1 , higher partial waves are still suppressed due to the centrifugal barrier, but the convergence of the partial wave expansion becomes slow. For example, to compute the NN differential cross section at E lab = 300 MeV at the 1% accuracy level, it is necessary to take into account partial waves up to the total orbital momentum of j max = 16 [55]. For very long-range interactions such as e.g. the magnetic moment interaction, the scattering amplitude is not converged even for j max ∼ 1000 [56]. One, therefore, may expect partial wave mixing effects in the FV calculations to be dominated by the longest-range OPE interaction. To get further insights into this issue and to validate this conjecture, we switch off all shorter-range interactions in the chiral EFT potential and consider the case of pure OPE.
In the upper row of figures 5 and 6, we show the corresponding positive-and negative-parity phase shifts extracted from the FV energies using the single-channel Lüscher formula, along with the results in the infinite volume. The deviations in the obtained phase shifts are qualitatively similar to those found in the calculations based on the full NNLO chiral potential. For S-wave dominated states, the single-channel Lüscher formula gives rise to a significant deviation for the L = 3 fm box, but it works reasonably well for L ≥ 5 fm. This is in line with the findings of ref. [22]. For P-wave dominated states, we again observe large deviations which persist even for the box size of L = 8 fm. Meanwhile, the deviations are found to increase with energies.
We stress again that increasing the box size does not allow one to improve the results of the single-channel Lüscher formula at higher energies. For example, focusing on the energy levels around 0.15 GeV in the first row of figure 6, one sees no obvious improvement of the results using the single-channel Lüscher method with increasing L. If one, however, considers the trajectory of a single state such as e.g. the ground state of the A − 2 irrep in the d 2 = 1 system as a function of the box size, one does observe a significant improvement when it approaches the threshold region for a sufficiently large L.
To unambiguously identify partial wave mixing effects as the source of failure of the singlechannel Lüscher formula, we consider the S-wave and P-wave projected OPE interaction. Specifi-  Figure 6. Phase shifts in the negative-parity channels extracted from the FV energy spectra using the P-wave Lüscher formula (various symbols) in comparison with the infinite-volume results for the OPE potential (upper row), the P-wave projected OPE potential (middle row) and the P-and F-wave projected OPE potential. Solid and dashed lines show the 1 P 1 and 1 F 3 phase shifts calculated in the infinite volume. For remaining notations see figure 3.
cally, we start with the partial wave expansion of the potential where z = cos θ with θ being the angle between the initial and final momenta and P l (z) denotes the Legendre polynomial, and define the potentials V S−wave (p, p ) = (4π) −1 V 0 (p, p )P 0 (z) and V P−wave (p, p ) = 3(4π) −1 V 1 (p, p )P 1 (z). Clearly, the resulting potentials do not generate any partial wave mixing effects when used to compute the FV energy spectra, so that the single-channel Lüscher approach is expected to become applicable. This is indeed fully consistent with our results, see the second rows in figures 5 and 6. For S-wave dominated states, the small deviation disappears after switching off the interaction in higher partial waves. For states with negative parity, the change is more pronounced. After switching off the interaction in higher partial waves, the P-wave phase shift is reproduced accurately with the single-channel Lüscher formula except for the smallest box size (presumably due to exponentially suppressed corrections). Finally, the lower row in figure 6 shows the effect of including the F-wave interaction, i.e. we consider V P,F−wave (p, p ) = 3(4π) −1 V 1 (p, p )P 1 (z) + 7(4π) −1 V 3 (p, p )P 3 (z). While we have not explicitly investigated the impact of even higher partial waves, a comparison of the upper and lower rows in figure 6 suggests that their mixing effects may also be significant at higher energies.

Phase shifts from FV energies using EFT
In section 5.2, we have shown in detail that partial wave mixing effects can not be neglected when calculating FV energies of two interacting nucleons at the physical pion masses. Such large partial wave mixing effects are caused by the long-range (pion-exchange) interaction, and they are responsible for the failure of the single-channel Lüscher approach to extract phase shifts from FV energy spectra. This will pose a significant challenge for future lattice QCD calculations in the NN sector close to the physical point. In the lattice QCD community, the issue is addressed by including several partial waves in the Lüscher's quantization conditions, see e.g. refs. [57,58]. When more than one partial wave is included, there is no longer a one-to-one mapping between energy levels and phase shifts, and one has to choose a theoretical framework to parameterize the T -matrix. As an alternative to the Lüscher method, one can benefit from the known model-independent OPE interaction using an EFT-inspired approach. Specifically, we propose to determine the short-range part of the NN interaction, parametrized in a systematic way by means of contact interactions, via a direct matching to lattice QCD FV energy levels. The method is, to some extent, similar to the low-energy theorems used in refs. [59,60] to restore the energy dependence of the NN scattering amplitude at unphysical pion masses.
To illustrate the method, we first define a toy model comprising the OPE and the heavy-mesonexchange potentials to generate synthetic data for the FV energies. Specifically, we consider where V 1π is the OPE potential considered in the previous sections (up to an S-wave contact interaction). For the various parameters, we choose the numerical values of M π = 139 MeV, F π = 92.4 MeV and g A = 1.26. For the heavy-meson-exchange interaction, we introduce both the isospin-triplet and isospin-singlet potentials with the same meson mass m h = 0.5 GeV. Further, c h1 and c h2 denote the corresponding dimensionless coupling constants. In order to regularize the UV divergences, we introduce a nonlocal Gaussian cutoff 17) and choose Λ = 0.45 GeV. The couplings c h1 and c h2 are adjusted in such a way that the toy-model interaction mimics the behavior of the NN 1 S 0 and 1 P 1 phase shifts as shown in figure 7. Using the toy model introduced above, we compute the corresponding FV energy levels, which are regarded as synthetic lattice data, see the right-most symbols in figures 8 and 9 for a given d 2 -value, where the results are only shown for the box with L = 5 fm. Having generated the synthetic data as described before, we are now in the position to describe our approach for extracting the corresponding phase shifts. To this aim, we exploit the knowledge of the long-range interaction in the underlying model to construct the EFT interaction We use the same regulator as in eq. (5.17) and employ the same value for the cutoff parameter Λ = 0.45 GeV. In order to fit the FV energy levels with positive parity, we introduce the interaction at NLO, whereC1 S 0 and C1 S 0 refer to the corresponding LECs. For the states with negative parity, we include the contact interactions up to NNLO, where C1 P 1 and D1 P 1 are the LECs. Since the contact interactions introduced above only contribute to the S-and P-wave channels, the FV partial wave mixing effects at the considered EFT order only arise from the OPE interaction.
To fix the LECsC1 S 0 , C1 S 0 , C1 P 1 and D1 P 1 , we employ the determinant residual method. For the considered toy-model example, we neglect the uncertainties of the synthetic data. First, we perform single-parameter fits by including only the dominant contact interaction in the corresponding parity channel, i.e. at LO (NLO) for positive-(negative-) parity states. In the second step, we also take into account the corresponding subdominant contact terms and perform two-parameter fits to the FV energies. As for the synthetic data, we only include the ground state energy of each irrep as input. Meanwhile, we ignore the energy levels of the d = (0, 0, 2) systems because they are identical to those of the d = (0, 0, 0) system in the non-relativistic case. For positive-parity channels, this leaves us with three and four energy levels for the boxes with L = 3 and 5 fm, respectively, up to the maximal energy of 0.23 GeV. For the negative-parity case, we have one and eight 2 energy levels for the boxes with L = 3 and 5 fm, respectively. Thus, for the smallest considered box size, our calculations for the negative-parity states are restricted to NLO in the EFT expansion due to the lack of input information.
In figures 8 and 9, we show the quality of the reproduction of the FV energy levels. The results from a single-parameter fit using the dominant contact interactions are shown by the dotted lines and the leftmost symbols for the considered d 2 value, while the dashed lines with symbols in the middle correspond to two-parameter fits including the contributions of the subdominant contact terms. As expected, one observes a clear improvement in the description of the energy levels by including the subdominant short-range interactions. The largest deviations from the synthetic data are observed for higher-energy states, namely for the third A + 1 state and the second B − 3 state in the d 2 = 2 box (which have not been used in the fits). This pattern is expected and reflects a slower convergence of the EFT at higher energies.
Having determined the values of the LECsC1 S 0 , C1 S 0 , C1 P 1 and D1 P 1 from the FV spectra as discussed below, we are now in the position to compute the phase shifts. This is achieved by solving the partial wave LSE in the infinite volume using the standard methods. In figure 10, we compare the phase shifts resulting from matching the EFT in the finite volume with the ones from the underlying model. We also show in this figure the phase shifts extracted using the singlechannel Lüscher formula. 3 For the 1 S 0 channel, the one-parameter fit at LO only qualitatively captures the behavior of the underlying phase shift. This can be attributed to the strongly finetuned nature of the interaction in this channel as reflected by the very large absolute value of the scattering length and to the well-known important role played by the range corrections. The fit results at NLO are strongly improved, and the underlying phase shifts are correctly reproduced even using the smallest considered box size of L = 3 fm. For the 1 P 1 channel, we only have a single data point for the smallest box at our disposal. Therefore, we can only perform a single-parameter fit for the L = 3 fm box. Using the single channel Lüscher formula, the resulting phase shift deviates strongly from the underlying one, thus pointing towards a significant contribution of higher partial waves to this energy level. Nevertheless, in the approach we propose, the one-parameter fit to this energy level already allows one to capture the qualitative behavior of the phase shift in this channel. In the larger box with L = 5 fm, we have more synthetic data at our disposal. The lower interacting energies as compared with the smaller box size allow for a more accurate determination of the LEC C1 P 1 in the NLO fit, which results in an improved description of the phase shifts at low energies as compared with the one using the smaller box of L = 3 fm. Including the subdominant contact interactions ∝ C1 S 0 , D1 D 2 , the description of both the 1 S 0 and 1 P 1 phase shifts improves considerably in a wide energy range. It should be emphasized that our input includes several energy levels dominated by higher partial wave components such as the ground state in the A − 2 irrep with d 2 = 3 and the ground state in the B − 3 irrep with d 2 = 2. Contrary to the single-channel Lüscher approach, which in these cases leads to large deviations, our method is insensitive to such FV partial wave mixing artifacts.

Application II: P-wave pion-pion scattering
We now turn to our second application and consider ππ scattering as an example of a relativistic system. In this exploratory example, we employ a simple phenomenological model for ππ interaction instead of using EFT. For a formulation of chiral perturbation theory with resonances see e.g. refs. [61,62].

Reduced Bethe-Salpeter equation in the finite volume
The relativistic Bethe-Salpeter equation is a four-dimensional integral equation. One has to reduce it into three-dimensional one to perform the plane wave expansion. In the literature, there are various approaches to achieve it [5,27,63]. In this work, we choose the Blankenbecler-Sugar equation as discussed e.g. in ref. [63] as our starting point. One can analytically perform the integration over the 0-th components of the relativistic two-body propagator as follows, where ω i (q) ≡ » m 2 i + q 2 and m 1 = m 2 = M π . The reduced BSE for the ππ scattering reads Notice that one can choose different reduced BSEs, but the procedure of performing the plane wave expansion is the same. We expand the reduced BSE in the plane wave basis with discrete momenta to obtain the matrix equation where the discretized propagator G is defined as Further, J is the Jacobi determinant arising from the transformation between the BF and CMF as shown in eq. (2.12), which is given explicitly in eq. (2.13) for two different schemes. The FV energy levels are obtained by solving the determinant equation The procedure for dealing with relativistic systems is more demanding than that for nonrelativistic ones since the effective interaction V in the three-dimensionally reduced eq. (6.3) is, in general, energy-dependent. Meanwhile, the Jacobi determinant can introduce additional energydependence. Such energy-dependence prevents one from reducing eq. (6.6) to the eigenvalue problem like in the non-relativistic case. In a special case when the interaction can be assumed to be energy-independent and the scheme-II is employed to relate the momenta in the BF and CMF as discussed in section 2.2, leading to the Jacobi determinant given in the second line of eq. (2.13), equation (6.6) reduces to the eigenvalue problem. This case is considered e.g. in refs. [29][30][31][32].
In this paper we do not make use of the above-mentioned simplifications and keep the energy dependence of the interaction and the Jacobi determinant. We decompose eq. (6.6) into decoupled matrix equations corresponding to different irreps. The FV energy levels are obtained by finding roots of the determinant equations for a given irrep Γ: 6.2 FV energy levels using a phenomenological model for P-wave ππ interaction We employ a phenomenological interaction model for P-wave ππ scattering rather than chiral perturbation theory. Specifically, we parametrize the interaction via an energy-dependent potential where f , M 0 and G V are free parameters. We introduce the CDD (Castillejo, Dalitz, Dyson) pole in the interaction [64], which naturally appears in the description of the vector form factor of the ππ system [65] based on the Lagrangians from refs. [66,67]. In this formalism, we neglect possible coupled-channel effects from the KK system. The interaction in eq. (6.8) corresponds to a short-range contact potential. The angular-dependent term p · p ensures that the interaction can only contribute to the P-wave, which switches off possible partial wave mixing effects in the finite volume. Thus, for such a short-range potential contributing to the P-wave alone, the single-channel Lüscher formula is perfectly applicable. The potential in the partial wave basis reads Notice that if we would make an on-shell approximation, we would obtain the interaction of ref. [28] (up to the 4π-factor), which was used to investigate finite volume effects in the same channel within the chiral unitary approach. In our normalization, the phase shift can be extracted from the partial wave T -matrix as Recently, lattice simulations of ππ scattering in the ρ-channel at the physical value of the pion mass were performed [68]. In our calculation, we choose the same box size of L = 4.3872 fm and pion mass of M π = 132 MeV as used in that work. Since the root-finding algorithm is very time-consuming, we choose a smaller cutoff q max = 1.5 GeV as compared to the one employed in ref. [28]. The three parameters f , M 0 and G V are determined by fitting the experimental ππ scattering phase shift [69,70]  Using the constructed model, we compute the FV energy levels of the ππ system in two different schemes to discretize the momenta as discussed in subsection 2.2. We then extract the phase shifts corresponding to these FV energy levels using the single-channel Lüscher formula as shown in figure 12. As expected, for the considered short-range interaction without partial waving mixing, the single-channel Lüscher formula leads to accurate results, which holds true for both considered discretization schemes. The outlier points with phase shifts about 180 degrees correspond to the non-interacting higher partial wave components. Further, in figure 13 we compare our calculated FV energy levels with those from lattice QCD simulations [68,71]. The differences of the FV energy levels from two discretization schemes appear to be tiny. Meanwhile, one can see the clear correspondence between our calculation and lattice QCD results for the low-lying energy levels in each irrep. For states below the inelastic threshold 4M π , our results are close to the lattice QCD ones.

P-wave ππ scattering from lattice QCD data for FV energies
We are now in the position to apply our approach for extracting the phase shifts from the FV spectra by fixing the model parameters f , G V and M 0 directly from the lattice QCD energy levels. To this  Figure 13. Comparison of the FV energy levels from lattice QCD and our calculations using the two different schemes for discretizing the momenta in the moving frames as explained in section 2.2. The notation for various symbols is the same as those in figure 12. The uncertainties of the lattice QCD data are taken from ref. [71]. end, we adopt the determinant residual method [33] and use the six energy levels below the 4M π threshold from the lattice QCD simulations of ref. [68] shown in figure 13.
We define the χ 2 -function as in eq. (4.3). To calculate the uncertainties σ[Ω Γ (E Γ,i )], one should propagate the errors of the energy levels from lattice QCD to the residual function Ω Γ (E Γ,i ), which is a tedious and time-consuming task. In our exploratory calculation we neglect possible correlations among the lattice data, which will have to be taken into account for more precise  Figure 14. P-wave ππ phase shifts extracted by matching the phenomenological model to the FV energy levels below the first inelastic channel computed in lattice QCD [68]. Experimental data are taken from refs. [69,70]. The red line shows the central result of the fit with the shaded area representing the fit uncertainty. The six lattice QCD energy levels with uncertainties are shown by various open symbols using the same notation as in figure 12. determinations of the phase shifts in future studies, and make an approximation where E upper Γ,i and E lower Γ,i are the upper and lower limits of the energy levels from the lattice data. We use the discretization scheme-I to fit the lattice QCD energy levels. We minimize the χ 2 -function using the package MINUIT [72] and obtain Having extracted the parameters of the model from the lattice QCD FV energy levels, we calculate the ππ P-wave phase shifts by solving the reduced BSE for the T -matrix in the infinite volume. In figure 14, the resulting phase shifts are shown, along with the uncertainty stemming from the errors in the extracted values of the model parameters in eq. (6.13). The pole mass and width extracted from our fitting results read, −195 MeV, Γ = 247 +114 −107 MeV. (6.14) As a comparison, the mass and width obtained within the inverse amplitude method in ref. [68] are 786 (20) MeV and 180(6) MeV, respectively. Though the uncertainties of the extracted phase shifts, mass and widths appear to be large, owing to the employed simplistic model and the restriction of the used lattice QCD data to energies below the first inelastic channel, these results confirm that our method can also be successfully applied to relativistic systems.

Summary and outlook
In this work, we propose an alternative approach to Lüscher's method for extracting two-body scattering information from FV energy levels. We employ the plane wave basis instead of relying on the commonly used partial wave expansion. Using the projection operator technique, we reduce the discrete plane wave basis into a direct sum of irreducible representations of the corresponding discrete groups. The FV energy levels are computed by solving the LSE or BSE for nonrelativistic or relativistic systems, respectively, for a given irrep and finding the poles of the resulting amplitude. The formalism is applied to both static and moving two-particle systems in finite periodic boxes. Since we do not use the partial wave expansion, all partial wave mixing effects due to rotational symmetry breaking in a cubic box are naturally embedded. We have used the above method to study spin-singlet NN scattering below pion production threshold as an example of a nonrelativistic system. Specifically, our goal was to study the impact of the long-range interaction due to the one-pion exchange on the severity of partial wave mixing effects when using the Lüscher formula to relate scattering phase shifts and FV energy levels. Throughout this study, we have restricted ourselves to the physical value of the pion mass. For S-wave dominated states, we found the single-channel Lüscher formula to lead to significant deviations for the smallest considered box size of L = 3 fm. For larger boxes with L 5 fm, the 1 S 0 phase shift is, however, well reproduced for energies below E lab 150 MeV, indicating that the contributions of the 1 D 2 and higher partial waves to the corresponding FV energy levels are negligible. On the other hand, we found severe partial wave mixing effects for P-wave dominated states, which originate from the longest-range interaction due to the OPE and make the single-channel Lüscher method inapplicable regardless of the considered box size (except for the near-threshold kinematics). While present-day lattice QCD calculations of the NN system are limited to unphysically large pion masses, this issue will pose a challenge once the lattice results will get closer to the physical point.
A natural and efficient solution to this problem can be achieved using the EFT framework, as it allows one to take into account the implications of long-range interactions in a systematic and model-independent way. To illustrate the method, we have considered a toy model of the NN interaction comprising the long-range OPE accompanied by the heavy-meson exchange potentials, which reproduces the qualitative behavior of the 1 S 0 and 1 P 1 neutron-proton phase shifts. In the absence of lattice QCD results for the NN system near the physical point, we used this model to generate synthetic data for FV energy levels. We then constructed the corresponding EFT featuring the same long-range interaction due to the OPE and involving contact terms with an increasing number of derivatives. Having determined the corresponding LECs by fitting them to the synthetic lattice QCD data, we have succeeded to reproduce the phase shifts of the underlying model using even the smallest box of L = 3 fm. All FV calculations are carried out using the discretized plane wave basis. The accuracy of the resulting phase shifts is systematically improvable by going to higher orders in the EFT expansion, and the method is completely insensitive to partial wave mixing effects that plague the applications of the Lüscher formula. Our results suggest that a precise determination of spin-singlet NN phase shifts from lattice QCD simulations at the physical point should be possible with box sizes no larger than L ∼ 4 . . . 5 fm. The lower bound on the box size in our approach is set by the requirement to have a sufficient number of resolved energy levels within the EFT applicability domain rather than by the interaction range as it is the case for the Lüscher approach.
As a second application, we considered P-wave ππ scattering. Partial wave mixing effects are not expected to be an issue for the ππ system in the FV, and our main motivation here was to demonstrate the applicability of our method in a relativistic setting. We employed a phenomenological model to parametrize the ππ interaction instead of relying on EFT. The three adjustable parameters of the model were tuned to reproduce the experimental behavior of the P-wave phase shift. We then computed the FV energy levels using the plane wave basis and compared them with those from lattice QCD simulations for the same pion mass and box size, finding a good agreement for the low-lying states in every irrep. We also applied our method to extract the phase shift by tuning the parameters of the model directly to the lattice QCD energy levels below the four-pion threshold. Though the uncertainties in the resulting determination are large, our central results agree well with the experimental ππ scattering phase shifts.
The considered examples provide a proof-of-principle that our approach can be successfully employed to extract the infinite-volume scattering information from FV energy levels accessible in lattice QCD. The essential ingredients of our method are (i) EFT or, possibly, phenomenological approaches to parametrize the interaction between the considered particles in a systematic way, (ii) the plane-wave basis in momentum space to avoid complications caused by partial wave mixing in the FV and (iii) matching to the FV energy spectra computable in lattice QCD. Compared with the conventional Lüscher approach, our method offers the advantage of being insensitive to partial wave mixing artifacts in the FV energy levels, and it also allows one to considerably reduce the box size by taking into account exponentially suppressed corrections to the Lüscher formula in the cases where the long-range interaction is known. Compared to alternative schemes discussed in literature such as the FV unitarized chiral perturbation theory for Goldstone boson scattering [14,24] and the effective Hamiltonian approach of refs. [29][30][31][32], our method allows for a clean decomposition of two-particle states into irreps of the FV symmetry groups and does not rely on the partial wave expansion known to converge slowly in the presence of long-range interactions. We have also addressed the complication due to a possible energy dependence of the effective interactions. In that case, one has to resort to the root-finding algorithm to obtain the poles of the T -matrix, which is computationally time consuming. We have explored an alternative technique to speed up the calculations by using the determinant residual method when tuning the parameters of the interaction to the FV energies.
The method contains two different ingredients, the EFT as the dynamics framework and plane wave expansion as the tool for numerical solution of FV energy levels. One can choose to combine one of them with other frameworks. For example, the plane wave expansion method in this work can be used to obtain the FV energy levels in the interaction of refs. [30,31]. The EFT framework in this work can be used to parameterize the multi-channel Lüscher quantization conditions. At the cost of more numerical effect, one might expect similar results to those in this work if the partial wave expansion is truncated such that the expansion of the amplitudes converges below a given energy. In the partial wave expansion, one should use the infinite-volume EFT framework to calculate NN phase shifts in different partial waves and then write down the multi-channel Lüscher equation containing all these partial waves. Using plane wave expansion, one can write down the FV quantization conditions directly without considering the scattering in infinite volume, which saves considerable effort.
As already pointed out above, we expect our method to show largest benefits for lattice QCD studies of systems featuring long-range interactions such as e.g. nucleon-nucleon scattering with the pion mass near or even smaller than the physical value. In the past few years, chiral EFT for the two-nucleon system has been developed into a precision tool, see e.g. [47,51,73,74] for recent precision studies. In ref. [51], a full fledged partial wave analysis of the available neutron-proton and proton-proton scattering data up to the pion production threshold has, for the first time, been performed within the framework of chiral EFT, thereby achieving a statistically perfect description of the experimental data. However, away from the physical quark masses, no experimental data are available, and one will have to rely on lattice QCD simulations, see [75][76][77][78][79][80] for examples of lattice QCD studies in the NN sector and ref. [81] for a recent review article. For heavy pion masses [77,82], lattice QCD results have already been matched to pionless EFT [83], see also ref. [81,84]. Our study provides a first step towards establishing an efficient and robust interface between lattice QCD and chiral EFT that will become necessary once the simulations at lower pion masses will become available. While we have restricted ourselves to spinless systems of two particles with equal masses in the present work, a generalization to nonzero spin and to particles with different masses in elongated boxes is straightforward. Work along these lines is in progress.

A Discrete groups
In this appendix, we briefly describe the elements, conjugacy classes and characters of the O h , D 4h , D 2h and D 3d groups needed in our analysis. For more details see textbooks about point groups, such as e.g. ref. [40]. In table 5

B Lüscher quantization conditions
For the spinless systems, the Lüscher quantization conditions read [6], ln,l n is an infinite-dimensional matrix. We truncate the angular momenta to l ≤ 2.
To make the expression concise, we adopt the short-hand notation, where Z P lm (1, q 2 ) are the conventional Lüscher zeta functions, which can be evaluated numerically [9].