Variational approach to $N$-body interactions in finite volume

We explore variational approach to the finite-volume $N$-body problem. The general formalism for N non-relativistic spinless particles interacting with periodic pair-wise potentials yields N-body secular equations. The solutions depend on the infinite-volume N-body wave functions. Given that the infinite-volume N-body dynamics may be solved by the standard Faddeev approach, the variational N-body formalism can provide a convenient numerical framework for finding discrete energy spectra in periodic lattice structures.


I. INTRODUCTION
Three-particle channels are important in spectroscopy of excited mesons. For example, G parity forbids a two pion decay of the a 1 (1260) meson so that the width of this axial vector resonance is determined by its coupling to three pions. The three pion decay channel is allowed for the exotic, J P C = 1 −+ , π 1 meson, while the exotic h 2 with J P C = 2 +− could decay to b 1 [ωπ]π and ρπ. Because of their possible hybrid nature and relation to confinement [1] these mesons are a subject of intense theoretical and experimental investigations [2]. In the baryon sector, the situation is even more complex because two and three-particle systems often mix. Phenomenologically it is observed that the ππN channels have significant influence and cannot be neglected when determining resonance parameters of nucleon excitations. The number of possible quantum numbers (π∆ ( * ) , πN * , σN , ρN ,. . . ) for a given spin and parity J P is large so that one has to truncate. Yet, phenomenology can be a guideline to select relevant channels, e.g., by inspecting the involved centrifugal barriers. For example, it is reasonable to assume that the Roper resonance N (1440)1/2 + is dominated by the σN channel (all three particles in relative S-wave) and maybe the π∆ channel, while the ρN channel could be less important at smaller energies where only pions from the low-energy tail of the ρ can be on-shell.
Pioneering calculations in lattice QCD of the spectra of light excited mesons [3] and baryons [4,5] deliver a semi-quantitative picture in which pion masses are large and finite-volume effects are usually neglected. Energy eigenvalues have been calculated for the a 1 (1260) [6] using qq operators but also a ρπ meson-meson operator, though taken at zero momentum, π(0) ρ(0). Recently, * pguo@csub.edu the HadronSpectrum collaboration calculated isospin I = 2 πρ scattering although the ρ meson is stable at the pertinent pion mass [7].
Few-body systems above threshold represent the next milestone for the ab-initio understanding of the strong interaction through lattice QCD calculations. Therefore, the infinite-volume extrapolation of three-body systems has attracted much interest recently , including an extension to coupled 2 to 3-body channels [18] and an study of connection between low temperature condensation and scattering in lattice φ 4 theory [31]. For the three-particle system in the finite volume, the energy eigenvalues are expected to change not only quantitatively but also qualitatively above threshold, compared to the two-body case. One can easily understand this by considering the non-interacting energies. For two particles of equal mass m in a cube of side lengths L with periodic boundary conditions, they are given by E = 2E n where E n = m 2 + p n 2 and p n = (2π/L) n, n ∈ Z 3 . In contrast, noninteracting three-body energies are given by E = E n + E n ′ + E n+n ′ where n, n ′ ∈ Z 3 . This pattern is entirely different from the two-body case which is expected to reflect in the interacting spectrum, as well. Apart from conceptual challenges, there is also the problem of underdetermination, i.e. there being a plethora of possible channels in contrast to scarce data, and, for the selection of relevant channels, similar considerations as in the infinite volume will have to be made.
Four-body systems above threshold have so many possible combinations of quantum numbers that one can only hope to address the simplest cases. A prime example is the energy region above the four-pion threshold in ππ scattering. In analyses of experiment, usually the fourpion channel below 1 GeV is neglected given that the experimentally measured partial wave amplitudes are almost elastic in this energy region (see Refs. [32,33] for a calculation of the inelasticity). Similarly, the existing lattice studies of ππ partial waves exclude the data above the four pion intermediate states for the infinite-volume extrapolation (see, e.g., Refs. [34][35][36][37]).
So far, in the three-body calculations only two mesonmeson or meson-baryon operators are used, and not yet three of type πππ or ππN . An exception is the N -body threshold calculation of positively charged pions by the NPLQCD collaboration [38,39]. Energy eigenvalues for the Roper [40] resonance above the ππN threshold have been calculated recently including non-local πN and σN operators. This calculation was performed close to the physical pion mass such that the extracted energy levels lie above the πN but also above the ππN threshold. An unusual pattern was observed that could originate from the afore-mentioned three-body dynamics in coupled πN , σN , . . . channel, where the σ has to be understood as an interacting two-body subsystem and not a stable particle.
In a rough and incomplete classification, there are three major approaches to solving the three-body problem in the finite volume, in the momentum space representation: An relativistic, all-orders perturbation theory pursued by Briceño, Hansen, Sharpe [12-14, 18, 21], a non-relativistic dimer formalism by Hammer, Pang, Rusetsky et al. [10,[15][16][17]20] and a method based on three-body unitarity to identify on-shell configurations and, therefore, power-law finite volume-effects by Döring and Mai [8,19,22]. For the latter two approaches, the partial diagonalization of the amplitude according to cubic symmetry was discussed in Ref. [22].
All these approaches aim at fully mapping out the few-body dynamics but there are also attempts to obtain essential information from these system without the need to take explicitly all degrees of freedom into account [25,26]. The first-ever prediction of excited threebody energy eigenvalues of a physical system (π + π + π + ), from two-body scattering information and lattice threshold eigenvalues [38,39] was achieved recently [8].

A. Variational approach
Most of the approaches to the few-body problem, including the examples discussed above, rely on momentum-representation of the reaction amplitude or correlation function in the finite volume. There is, however, an alternative approach based on Faddeev equations and two-and three-body wave functions in configuration space [27-30, 41, 42]. The finite volume wave function is related to the infinite volume wave function by a linear superposition over infinite sets of periodic cubic boxes, the quantization conditions are subsequently obtained from matching conditions [27-30, 41, 42]. One of the advantages of this approach is that the connection between long-distance correlation over boxes and shortdistance interaction within each cubic box is made explicit by construction of finite volume wave function.
In the present work, we set up and explore a foundation to a potentially convenient numerical approach to the N -body interaction in finite volume. The formalism presented in this work is based on the variational method [43,44] combined with the Faddeev approach [45][46][47][48]. A brief summary of the variational principle and Faddeev approach are presented in Appendix A and B as a short reference for a reader who is not familiar with above mentioned methods.
Based on the traditional variational principle, the secular equation may be obtained by considering: δΛ = 0 with Λ = Φ|E −Ĥ (L) |Φ , where both the N -body HamiltonianĤ (L) =T + i<jV (L) (ij) and the trial wave function Φ display periodicity in the cubic lattice,T is the kinetic energy operator andV (L) (ij) stands for the periodic pair-wise interaction between the i-th and j-th particle. Instead of the traditional approach, in this work we write the total wave function as a sum of multiple terms [45][46][47][48] (ij) |Φ . Therefore, a single Schrödinger equation, (E −Ĥ (L) )|Φ = 0, is turned into N (N − 1)/2 coupled equations. The advantage of Faddeev approach is to allow one to incorporate the dominant subsystem structures in an adequate way and to use two-body scattering amplitudes as input for the N -body dynamical equations. By splitting the complete N -body wave function, ultimately N (N − 1)/2 secular equations may be obtained by considering (1) As will be shown in Section III, the secular equations obtained from Eq. (1) resemble two-body secular equations.
In solid state and condensed matter physics, the Linear Combination of Atomic Orbital (LCAO) method for the calculations of the electronic structure of periodic systems [49] provides a elegant way to construct a wave function that satisfies periodic boundary conditions. The trial wave function that describes electrons traveling in a periodic crystal is constructed by linear superposition of all atomic orbital solutions of an isolated atom centered at each cell of the crystal. In this way, the periodic boundary conditions of the wave function is automatically guaranteed, and the energy spectra are given by secular equations from the variational principle. Similarly as also suggested in Refs. [27][28][29]42], the trial finite-volume N -body wave function may be constructed from the N -body infinite-volume wave function, where Ψ is the solution of the corresponding infinite-volume Schrödinger equation, (E −Ĥ)|Ψ = 0, and the set {x} stands for the complete set of particle positions. Hence, the trial finitevolume N -body wave function satisfies periodic boundary conditions by construction. In principle, the infinitevolume wave function may be solved by standard methods, such as Faddeev's approach [45][46][47][48].
The immediate gain of the variational approach for finite volume systems is evident. (1) The construction of the finite-volume wave function from the infinite-volume wave function presents a clear connection between shortrange N -body dynamics within each image of the cubic box and long-range correlations in the entire periodic lattice structure. The N -body dynamics within each cubic box is determined by N -body Faddeev equations. The long-range correlation effect accumulated from all cubic boxes is implemented by the linear superposition of all wave functions centered at each image of cubic boxes.
(2) The quantization conditions (secular equations) due to the periodic structure of the lattice are imposed by a variational approach, and eventually yield the discrete energy spectra of the system in the finite volume. (3) The variational formalism is mathematically transparent and may be suitable for the numerical evaluation of the N -body finite volume problem. However, it comes at the price of sacrificing the explicit analytical expressions of quantization conditions as present in the Lüscher formula [50] for the two-body problem and at the high cost of computation of the N -body phase space integration.
As we present here the first attempt to calculate the N -body interaction in finite volume with the proposed methods, and also for the sake of simplification of the discussion, in this work, we only consider a simple model with N non-relativistic spinless particles interacting through pair-wise short-range potentials. Given the infinite-volume wave function, Ψ, that may be solved by the Faddeev approach and is used as input to the finitevolume problem, the quantization conditions for the finite volume N -body interaction are obtained and presented in Section II. In principle, three-body forces and coupled-channel effects may be included in the formalism as well. However, such type of effects are not considered in the present work. In addition, bound states below threshold may be described by the analytic continuation of scattering amplitudes. In the present work, the focus is only given on the presentation of the N -body problem with pair-wise interactions. The discussion of threebody force, coupled-channel effect, and bound states below threshold will be given in the future publications.
An effective two-body formalism may be accomplished by integrating out the rest of the degrees of freedoms of the N -body system. The resemblance of the N -body quantization condition to the two-body quantization condition will be discussed in Section II E. At last, in the Appendix we also present some main results for a pairwise short-range δ-function potential. The renormalization issue due to the singular nature of the δ-function interaction in 3D is also discussed.
The paper is organized as follows. In Section II we present the derivation and main results of the variational approach to the N -body interaction in finite volume. The discussion and summary are given in Section III. The dynamics of N non-relativistic particles interacting via pair-wise interactions in the finite volume with periodic boundary conditions is determined by the Schrödinger equation,  (2) where we take all masses to be equal to m, x i denotes particle position and {x} ≡ {x 1 , x 2 , · · · , x N }. The potential between the i-th and j-th particles is described by where V (ij) is the potential between the i-th and the j-th particle in the same box, and L is the size of the threedimensional cube. Because of the periodicity of the finitevolume potential, the finite-volume N -particle wave function, Φ({x}), must also satisfy the periodic boundary condition, where {n x } = {n x1 , n x2 , · · · , n xN } and n xi ∈ Z 3 . Following the Faddeev approach [45][46][47][48], that is briefly summarized in Appendix B, the N -body finite-volume wave function may be expressed as the sum of N (N −1)/2 terms, with wave function Φ (ij) required to satisfy the equation 2m . Hence, the Schrödinger equation, Eq. (2) is converted to N (N − 1)/2 coupled equations. The N -body finite volume wave function is normally expanded in terms of a set of periodic basis functions, where [J] refers to a complete set of quantum numbers for the N -particle basis wave functions, and c [J] stands for the expansion coefficients that may be determined by the variational principle. Similar to Eq. (5), the basis function Φ [J] is also given by the sum of N (N − 1)/2 terms: [J] . We remark that the choice of basis functions may be arbitrary, depending on the symmetry of the specific physical system and the convenience of numerical computation. However, the different choices of basis functions should lead to consistent results. The reasonable choice of basis functions should preserve the symmetry of physical system, such as periodic boundary conditions in finite volume, and numerical results should be stable and converge.
In the finite volume, to fulfill the periodic boundary condition, the basis functions, Φ [J] , are constructed from infinite volume solutions of the Schrödinger equation. The complete details are given in following Section II B. In the following we use the two-body problem as a specific example to explain our choice of basis functions Φ [J] in finite volume. We denote the infinite volume wave functions by Ψ [J] . In the case of two-particles in the center-of-mass frame, the solutions of the Schrödinger equation in the infinite volume may be determined by the partial wave expansion of the free incoming wave: = JM are partial wave quantum numbers, r stands for the relative position of the two particles, and q = √ mE is the relative momentum of the particles. Thus, the asymptotic solution of two-body scattering in infinite volume is given by where t J (q) stands for the two-body scattering amplitude. Hence, Ψ [J] can be used as basis functions in the infinite volume to construct the corresponding finitevolume basis functions, Φ [J] in order to fulfill the requirement of periodic boundary conditions, Such finite volume basis function reflect the periodic lattice structure. It will be shown in Section II D, that the two-body basis functions, Φ [J] , will be related to Lüscher's zeta function [50]. As a general remark, the basis functions should be constructed respecting the cubic symmetry of the problem; yet, we are not explicitly addressing this topic in here. Such basis functions will be given by linear superpositions of infinite-volume basis functions; see, for example, Ref. [22] in which basis functions for "shells" are derived as linear combinations of cubic harmonics which by themselves are superpositions of spherical harmonics. With the finite volume wave function expanded in terms of the basis function, Eq. (7), the variational prin-ciple, ∂Λ (ij) /∂c * [J ′ ] = 0, with Λ (ij) given by yields a set of coupled secular equations, and summing the above N (N − 1)/2 equations leads a secular equation of familiar form (see also Eq. (A5)), are satisfied simultaneously. To summarize, given a set of basis functions, Φ (ij) [J] that satisfy periodic boundary conditions, the variational solution for the spectrum of the N -body finite volume problem is reduced to a solution of the determinant conditions given by Eq. (13).

B. Construction of finite volume wave function and reduction of the secular equations
To proceed, a proper choice of basis functions, with periodic symmetry, Φ (ij) [J] , has to be made. There exist numerous examples of periodic basis functions, e.g. devised for calculations of electronic structure in a periodic lattice and other condensed matter or solid state systems. Specifically, in the LCAO method [49], the periodic variational basis functions for describing electronic states tunneling in a crystal are constructed as a linear superposition of all atomic orbital solutions for an isolated atom located at each unit cell of the crystal. In application to hadron scattering, instead of using bound state solutions, as done previously in [27][28][29]42], the basis functions may be constructed from the linear superposition of infinite-volume scattering wave functions centered at each image of the cubic box, where {n x } = {n x1 , n x2 , · · · , n xN }, N = n∈Z 3 is a normalization factor, and Ψ (ij) [J] satisfies the Schrödinger equation in the infinite volume with potentials V (ij) , The total N -body wave function in the infinite volume, Ψ [J] , is the scattering solution of the Schrödinger equation, (E −T −V )|Ψ [J] = 0, and may be expressed as where Ψ [J] stands for the free incoming wave and satisfies the free Schrödinger equation: For example, the two-body partial wave free incoming wave are given by Ψ while the corresponding wave function in the three-body case in given in [27]. The periodic symmetry of the finitevolume wave function Φ is satisfied automatically by the construction in Eq. (14). The infinite-volume N -body wave function, Ψ, may be solved by a standard Faddeev approach. Locally, at each image of the cubic box, assuming the size of the box is large enough, the short distance N -body dynamics is thus described by the solution of the infinite-volume Schrödinger equation, Ψ. The long distance correlations at scales larger than the range of the potential is taken in to account through the linear superposition of the infinite-volume wave functions centered at each image of the periodic cubic box, see Eq. (14).
Using Eqs. (14), (15) and taking into account periodicity of the finite-volume potential, the secular equations given by Eq. (11) reduce to where and {n x } (i,j) stands for the set {n x } excluding two elements: n xi and n xj . Furthermore, the N (N − 1)/2 determinant conditions given by Eq. (13) become Given the scattering solution of the N -body problem in infinite volume, Ψ, as input, the finite volume wave functions, Φ and X (ij) , can be constructed by using Eq. (14) and Eq. (18)  , momenta dependence remains in finite volume wave functions as well. Therefore, using the infinite volume scattering solutions, Ψ [J] , as inputs, the finite volume quantization conditions, Eq. (19), yield discrete energy spectra as the consequence of periodic cubic lattice structure. In other words, the discrete energy spectra in finite volume is the result of long distance correlation effects of particles in a periodic lattice structure. Meanwhile, the specific patterns of discrete spectra rely on the short distance interaction that is described by scattering amplitudes or "amplitudes carrier" wave functions, Ψ [J] , in infinite volume. Ultimately, the quantization conditions play the role of imposing constraints on energy spectra due to periodic boundary conditions. The scattering information in infinite volume and periodic boundary conditions are combined together by finite volume wave functions, Φ [J] and X (ij) [J] . The secular equations, Eq. (17), and corresponding quantization conditions, Eq. (19), can be further reduced by removing the center-of-mass motion. However, the choice of relative coordinates for the N -body system normally can be made quite arbitrary. In Section II C, the specific choice [48] is described by forming the succession of subsystems of N particles in such a way that the subsystems are obtained by successive joining of particle-2, particle-3, · · · , particle-N to particle-1, i.e. (12), ((12)3), (( (12)3)4), · · · .

C. Removal of center of mass motion
Since in this work, we are mainly interested in scattering solutions, the set of incoming particle's momenta is also introduced to label the N -body wave functions: {p} = {p 1 , p 2 , · · · , p N }, where p i stands for the momentum of the i-th incoming particle. The total energy of N particles is given by E = N i=1 p 2 i /2m. The center-ofmass motion and internal motion of the N -body system can be separated out by changing variables of coordinates system. One particular choice of the relative coordinates and momenta [48] may be given respectively by where n < N . The quantities ρ (12),n and q (12),n may be interpreted as the relative coordinate and momentum between the (n + 1)-th particle and the center of mass coordinate or momentum of the cluster of particles (1, 2, · · · , n), respectively. The index (12) is used to label the special choice we made on relative coordinates and momenta, so that for n = 1, the relative coordinate and momenta are defined between particle-1 and particle-2, The complete sets of relative coordinates and momenta are given by {ρ (12) } = {ρ (12),1 , ρ (12),2 , · · · , ρ (12),N −1 } and {q (12) } = {q (12),1 , q (12),2 , · · · , q (12),N −1 } respectively. The center of mass coordinate and momentum of the N -body system are The total energy of the N -particle system is given by E = N −1 n=1 q 2 (12),n /m + P 2 /2mN , and the kinetic energy operator of the N -body system iŝ In terms of the set {ρ (12) }, the infinitevolume N -body wave function is given by where the center of mass motion is represented by a plane wave, e iP·R , and the wave function ψ describes the internal motions of the N -body system. The same applies to the finite-volume wave function. The N -body wave function representing the relative motion in finite volume is thus given by where the n-th element of the set {ρ (12) } is ρ (12),n = ρ (12),n + 2n n + 1 1 n n k=1 n x k − n xn+1 L. (25) In order to remove one redundant element of the set {n x }, a subset {n (12) } is introduced, (12),1 , n (12),2 , · · · , n (12) where the k-th element is given by n (12),k = n x k − n x k+1 . With the introduction of the set {n (12) } in this particular way, all the elements of set {n (12) } still belong in Z 3 . We also find ρ (12),n = ρ (12),n + 2 n(n + 1) n k=1 kn (12),k L, and Hence, after removal of the center-of-mass motion, the finite-volume relative wave function, φ, is related to ψ by and φ satisfies periodic boundary condition, φ({ρ (12) }, {q (12) (12),k )L φ({ρ (12) }, {q (12) }). (30) Clearly, the set {ρ (12) } cannot be the only choice of independent relative coordinates. By exchanging the labels of particle, such as 1 ↔ i and 2 ↔ j in set {ρ (12) }, another independent coordinate set, {ρ (ij) }, may be obtained. For example, relabeling 2 ↔ 3 in {ρ (12) }, the elements of set {ρ (13) } are given by The different sets {ρ (ij) } and {ρ (i ′ j ′ ) } are linked by a linear transformation, for example, Because of the way of construction, the transfer matrix Γ between two different sets is an orthogonal matrix, Γ T Γ = I. Hence, both the sum, With introduction of these sets of relative coordinates, the finite volume N -body wave function, φ [J] , may be written as, (12) }, {q (12) }), (35) and we introduced rescaled total energy and kinetic energy operators, [J] is given by where {n (ij) } = {n (ij),1 , n (ij),2 , · · · , n (ij),N −1 }, and the set {n (ij) } may be related to {n (12) } by relabeling 1 ↔ i and 2 ↔ j. The set {ρ (ij) } is defined bÿ ρ (ij),n = ρ (ij),n + 2 n(n + 1) n k=1 kn (ij),k L, (n (ij),k ∈ Z 3 ). (37) The wave functions, ψ (ij) [J] , are the solutions of infinitevolume Schrödinger equations, (12) }, {q (12) }), (38) and the total infinite-volume N -body wave function is given by, [J] refers to the incoming free wave.

D. Relation to two-body Lüscher formula
In the simplest case, N = 2, there is only one relative coordinate and relative momentum so that the particle index can be dropped and in what follows they are denoted by r and q, respectively. Notice that χ (ij) [J] is now reduced to ψ [J] for the two-body system. Therefore, the determinant condition in Eq. (39) is given by where [J] = JM stand for quantum numbers of a specific partial wave, and Assuming the potential is of finite range, V (r) = 0 for R < r < L, in the region outside of the potential the infinite-volume wave function has the form, where t J is the partial-wave two-body scattering amplitude, and t J is normalized by the unitarity relation, Im [1/t J ] = −1. Inside the potential region, the wave function is given by the solution of the Schrödinger equation, and is denoted by ψ (in) [J] (r, q) from now on. The matrix element of the determinant condition, Eq. (45), is given by Using the fact that and the relation [J],[j] , we find Therefore, the determinant condition for the two-body problem, Eq. (45), yields Lüscher's formula, [J],[j] (q) = 0.
Hence, the variational approach to the finite -volume fewbody system is consistent with the Lüscher approach.

E. Effective two-body formalism
In this section, we would like to show that the N -body quantization conditions may be recast in a similar form as the two-body quantization condition in Eq. (45).
Before proceeding to N -body interaction, first of all, let us rewrite the two-body quantization condition of Eq. (45) by using the periodicity of the finite-volume wave function. We find where χ When all the interactions V (i ′ j ′ ) except the interaction between the i-th and j-th particle are turned off, clearly, the quantization condition in the (ij) channel in Eq. (57) is thus reduced to the two-body quantization condition given in Eq. (54).

III. DISCUSSION AND CONCLUSION
A. Resemblance to isobar model In the past, isobar models [51][52][53][54][55][56][57][58][59][60][61][62][63][64] have been a useful tool to describe few-body interactions, in which the fewbody interaction is treated by taking into account all possible recombinations of two-body subsystems. The twobody subsystems are considered as the dominant contribution compared to three-body force, and the few-body interaction correction to the two-body subsystem is generated by rescattering between all possible pairs. In order to show the similarity of this approach to the isobar formulation [51][52][53][54][55][56][57][58][59][60][61][62][63][64], we consider a special case, i.e., a three-body system with two light spinless particles and one infinitely heavy spinless particle stationed at the origin. The heavy particle is labeled as third particle. The system may be described by where V (L) represents the interactions between the heavy particle and one of the light particles, and U (L) stands for the interaction between the two light particles. Again, we use the superscript (L) to identify the periodic potential in finite volume. We also assume that U (L) is a weak interaction, so that the interaction between the two light particles is treated as a perturbation, which serves the purpose of discussion in this work. Therefore, the corresponding infinite-volume wave function must have the form of ψ(r 1 , r 2 ; q 1 , q 2 ) = ψ(r 1 ; q 1 )ψ(r 2 ; q 2 ) + δψ(r 1 , r 2 ; q 1 , q 2 ), (59) where the first term is the solution of the system with zero interaction between two lights particles. The second term, δψ, can be considered as perturbative contribution when the weak U -potential is turned on. The two-body infinite-volume wave function, ψ(r i ; q i ), is given by the solution of with E = q 2 1 /2m + q 2 2 /2m. Following the argument provided in previous sections, the finite-volume wave function is constructed from the infinite-volume wave function. The two ingredients φ and χ (1/2) of the secular equations thus also have the forms, φ(r 1 , r 2 ; q 1 , q 2 ) = φ(r 1 ; q 1 )φ(r 2 ; q 2 ) + δφ(r 1 , r 2 ; q 1 , q 2 ), (61) and (62) where φ(r i ; q i ) = ni∈Z 3 ψ(r i + n i L; q i ). The construction of δφ and δχ (i) can also be obtained based on δψ, accordingly, but the specific expressions are not crucial for our brief discussion. The two secular equations may also be treated as a perturbation; for example, for channel (13), we obtain where δU stands for the perturbative contribution to the secular equation from the weak U -potential. The threebody quantum numbers set, [J], is constructed from twobody quantum numbers, [L i ], by [J] = [L 1 ] [L 2 ]. In the limit of δU → 0, two secular equations yield two independent two-body quantization conditions, and both q i 's are quantized independently according to the corresponding Lüscher formula. The three-body correction with weak U -potential may be obtained by perturbation. In the current approach, the physical picture is presented in a way that is similar to the three-body rescattering effect corrected isobar model, [51][52][53][54][55][56][57][58][59][60][61][62][63]. That is to say that the three-body system considered in this subsection can be treated as two two-body isobar subsystems: (13) and (23), which yield two independent two-body quantization conditions to q 1,2 . The threebody rescattering correction to isobar subsystems produces an energy shift on quantized two-body energies: q 2 1 /2m + q 2 2 /2m + δE.

B. Summary
In this work, we propose a variational approach to the finite-volume N -body problem. In order to fulfill the periodic boundary conditions, the trial wave functions are constructed by linear superposition of all the solutions of the infinite-volume wave functions centered at each image of the periodic cubic boxes, given that the infinite-volume wave functions may be obtained by standard methods. In this approach short-range N -body dynamics is local to each box and long-range correlations correspond to particles travel through the entire periodic structure of the lattice. In other words, the short-range dynamics are determined by infinite-volume wave functions, and finite-volume wave functions control the long-range correlations that eventually yield the discrete energy spectra because of the periodic structure of lattice. No explicit analytic expressions between discrete lattice eigenvalues and the scattering amplitude, such as present in the twobody Lüscher formula, can be given by the variational approach for N > 2. Instead, the discrete energy spectra and N -body scattering amplitudes are linked in a rather complicated way. Nevertheless the method is has potential advantages for systems with N >> 2. By combining the variational approach with the Faddeev approach, the N (N − 1)/2 quantization conditions are obtained ultimately. In the end, these quantization conditions can be expressed in a way that resemble rescattering in the isobar approach, The overall rescattering corrections can also be written as a collective effect by integrating out the of degrees of freedom other then as selected two-body subsystem.

IV. ACKNOWLEDGMENTS
We acknowledge support from the Department of Physics and Engineering, California State University, Bakersfield, CA. P.G. would like to thank Vladimir Gasparian and David Gross for numerous fruitful discussions. For a complex system, most calculations are based on approximate methods, the variational principle is one of most commonly used approaches. In this section, we briefly outline the main idea of the variational principle as the approximate solution to a general quantum system [43,44], which satisfies the Schrödinger equation, The trial wave function Ψ may be expanded in terms of a set of basis functions which satisfy certain boundary conditions or symmetries of the system, The solution of the Schrödinger equation is thus given by the variational principle: Thus, the variational principle yields secular equations In this section, for completeness, we give a brief summary of the Faddeev equations for the interaction of N particles in general. The non-relativistic N -body dynamics is described by the Schrödinger equation, Assuming only pair-wise interactions among particles, V = N (i<j)=1V (ij) , the scattering solution of N -body wave function is normally written as the sum of multiple terms [45][46][47][48], where Ψ (0) stands for the free initial incoming wave of the N -particle state: (E −T )|Ψ (0) = 0, and Ψ (ij) satisfies the equation In this way, the N -body Schrödinger equation is turned into N (N − 1)/2 coupled equations, where the Green's function operatorĜ (ij) is given bŷ TheĜ (ij) is related to the two-body scattering amplitude,t ij , byĜ whereĜ (0) = (E −T + iǫ) −1 stands for free two-body Green's function. The total N -body scattering amplitudeT is introduced byT and Eq. (B4), we find that thê T (ij) satisfy the coupled equationŝ TheT (ij) and |Ψ (ij) are related in a simple form by One of the advantages of the Faddeev approach is that it demonstrates the general relations between the subsystem amplitude and the N -body amplitude in a natural way, see Eq. (B8), and the two-body amplitude is used as input for N -body dynamics. In addition, the unitarity relation of the N -body scattering amplitude is also automatically guaranteed by the Faddeev equations.
Appendix C: Solutions of N -body problem with δ-function potentials To give readers a concrete example of our proposed approach to N -body finite volume problem, in this section, we also include a specific example of short-range interactions between two particles with a δ-function potential.

Two-body interaction with δ-function potential
We consider two-particle scattering with a pair-wise δ-function potential, assuming that all particles are spinless and have equal mass. The bare strength of the δfunction potential between the two particles is described by V 0 . With the same convention as in Section II D, the infinite-volume wave function is given by the Lippmann-Schwinger equation, where G (0) stands for the free two-body Green's function, The two-body Lippmann-Schwinger equation, Eq. (C1), exhibits an solution, where only S-wave contributes to the two-body scattering amplitude, It has been known that singular potentials, such as the δ-function potential, require in two or higher dimensions regularization and renormalization [65]. Adopting the renormalization scheme proposed in Ref. [65], a renormalized strength of the δ-function potential, V R , is introduced to absorb the divergent part of h (+) 0 (qr)| r→0 = 1 − i qr , and renormalized and bare strengths are related by In terms of the renormalized quantity, the S-wave twobody scattering amplitude now reads and the unitarity relation of the two-body scattering amplitude is guaranteed by Im[t −1 0 ] = −1.
In the finite volume, based on the discussion given in Section II D, we find that The non-trivial solutions of the δ-function potential are given only by S-wave, which are determined by N -body interaction with δ-function potential

body secular equations
Given all the ingredients of the finite-volume wave functions from Eq. (C23) up to Eq. (C27), the discrete energy spectra are thus given by N (N − 1)/2 secular equations, As ρ (ij),1 → 0, both χ and φ (ij) [J] appear to have the exact same ultravioletly divergent behavior because of the spherical Hankel function, (ij),n ρ (ij),1 ) ∼ 1 ρ (ij),1 , thus, the subtraction of two terms is completely free from ultraviolet divergence. The ultraviolet divergence appears in φ * as well. Since V 0 is a bare parameter, in order to remove all the ultraviolet divergence, a renormalization of the energy is required. Notice the fact that after the shift on the energy: E → E + δE, the divergent contribution of secular equation can be completely canceled out by the counter term, δE. Thus the renormalized quantity φ * [J ′ ] mV 0 may be given by where T (ij) [J] q; q (12) , q 3 = 4πt 0 ( σ 2 − q 2 ) σ 2 − q 2 × dr k e −iq·r k ψ (0) [J] (0, r k ; q (12) , q 3 ). (C35) In terms of the T -amplitude, the finite volume threebody wave function is taken from Eq. (C23), φ (ij) [J] (r (ij) , r k ; q (12) , q; q (12) , q 3 ) .

(C42)
Normally, the infinite sum of Bessel functions in the above equation has poor convergence, and for numerical purposes a better expression may be given by Eq. (50) in a partial-wave expanded form. It may be also convenient to use an alternative form in Eq.(D7) without partial wave expansion, and also splitting the integration by an arbitrary parameter, η, we can rewrite the lattice sum of the Green's function in Eq. (50) as (D2) For the first term in Eq.(D2), using the identity 1 2π √ π e −r 2 t 2 = 1 2t 3 dq ′ (2π) 3 e − q ′ 2 4t 2 e iq ′ ·r , and also applying Poisson summation, we find For the second term in Eq. (D2), except for n = 0, the convergence of the integration is well-defined for a finite value of η, thus, we would like to isolate the n = 0 piece, with the help of the identity − 2i √ πq ∞ η dte −r 2 t 2 + q 2 4t 2 = in 0 (qr) + i e −iqr erf(− iq 2η + rη) − e iqr erf(− iq 2η − rη) 2qr .
(D5) Therefore, we find Putting everything together, we thus obtain a fast convergent expression of the lattice sum of the Green's function without partial wave expansion, where