A lattice model of heavy-light three-body system

We present a study of a $1+1$ dimensional heavy-light three-body system in finite volume. The heavy-light system is simulated by a coupled-channel $\phi^4$ type lattice model, and both ground state and excited states of multiparticle energy spectra are measured on various lattices. The lattice simulation data analysis is performed based on variational approach.


I. INTRODUCTION
Much of strong interaction phenomenology manifests itself in few-body systems. Due to the many degrees of freedom their quantitative description is more complicated than in the two-body case. Lattice QCD calculations, e.g., of multi-pion systems [1][2][3], provide an ab-initio understanding of few-body dynamics. However, the dynamics is encoded in a discrete spectrum of energy eigenvalues corresponding to the cubic volume with usually periodic boundary conditions in which these calculations are performed. In the past few years, much progress toward mapping such spectra to infinite-volume amplitudes has been made . Some approaches provide connections to the infinite volume without fully resolving the few-body dynamics [28][29][30][31]. In first applications of infinite-volume mappings, the two-and three-body lattice QCD spectra by the NPLQCD collaboration [1,2] was analyzed in Ref. [32], and the recent data of Ref. [3] was analyzed/predicted in Ref. [33] and Ref. [34], respectively.
Most of these developments are along the line of building connections between reaction amplitudes in infinite volume and long-range correlations due to the periodic structure of a finite box. To a certain extent, these developments may be regarded as extensions of the Lüscher formula [35][36][37][38][39][40][41][42][43][44][45][46]. In the two-particle sector, Lüscher's formula demonstrates a clear separation of two physical scales: (i) short-range physics in a single box that is parametrized by a scattering amplitude, and (ii) longdistance correlations in a periodic lattice structure that is described by zeta functions [35]. Although formulating quantization conditions by using reaction amplitudes as input to producing discrete energy spectra or vice versa presents a more conventional foundation of dealing with multiparticle dynamics, one is confronted to deal simultaneously with questions regarding infinite and finite-volume physics. As pointed out in Ref. [25], the quantization condition of multiparticle system in finite volume may be constructed directly from Faddeev * pguo@csub.edu † doring@gwu.edu equation-type coupled-channel integral equations. The discretized finite-volume wave functions may be treated as coefficients of secular equations, hence, the quantization condition given by the determinant of secular equations is free of infinite volume reaction amplitudes, and is presented in terms of finite-volume Green's function and particle interactions. By avoiding the most difficult part of multiparticle dynamics, the quantization condition may be more efficient for practical data analysis of lattice calculation results. The aim of this work is to demonstrate the feasibility of the variational approach, and explore a robust form of a quantization condition which depends only on the lattice structure and interaction potentials. The interaction potentials may be parameterized and treated as inputs to fit the lattice simulation results. Once the parameters of interaction potentials are extracted, in principle, all the physics information is complete, and the dynamical information, such as scattering amplitudes, may be computed in a separate step.
In the present work, using a non-relativistic heavylight three-body system with short-range pair-wise interactions as a pedagogical example, we show in details that the three-body problem in finite volume may be turned into a coupled two-body like homogeneous Lippmann-Schwinger equations. The quantization conditions may be obtained in terms of only three-body finite-volume free Green's function and interaction potentials, and are free of the specific form of finite-volume variational basis functions. In the second part of this work, a coupled-channel φ 4 theory lattice model is used to simulate the heavylight three-body system proposed in this work. The multiple levels of multiparticle spectra are extracted from lattice simulation. In order to compensate finite lattice spacing, finite volume and relativistic dynamics in lattice simulation, the non-relativistic formalism derived in the first part of the presentation is thus reformulated to a relativistic lattice version for the practical fit of lattice simulation results. At last, the quantization conditions are applied to extract the parameters of the heavy-light three-body system: the mass of particles and the coupling strength of short-range interactions.
The paper is organized as follows. Using a heavy-light system as a specific example, the variational approach is explained in detail in Section II. A lattice model of the heavy-light system, using a Monte Carlo updating algorithm of the lattice model, and the construction of multi-particle operators and multi-particle spectra in lattice simulation are described in Section III. The strategy of data analysis, reformulating the finite-volume heavylight three-body system to a relativistic lattice version and numerical results are presented in Section IV. The summary and outlook are given in Section V.

II. A NON-RELATIVISTIC HEAVY-LIGHT THREE-BODY SYSTEM IN FINITE VOLUME
Although the non-relativistic formalism may not be the most suitable framework for describing relativistic lattice simulations, it still provides a clean and simple presentation of physics in a finite box, and is able to capture all the key features of finite-volume dynamics without extra complication of relativistic effect. Therefore, in this Section, all the presentations are given in a non-relativistic framework.
We start our discussion with the simple example of a two-body finite volume system in 3D. The wave function of two identical bosons in the center-of-mass frame satisfies the homogeneous Lippmann-Schwinger equation, where q stands for the relative momentum of the two identical particles and periodic the Green's function satisfies The wave function in finite volume has to be periodic as well, such as, Φ(r + nL; q) = Φ(r; q). According to the variational approach [25], the solution of the homogeneous Eq. (1) is given by the linear superposition of independent solutions in infinite volume, where Using the asymptotic form of the infinite-volume wave function, is the partial-wave scattering amplitude, and also performing the partial wave expansion of G L [J], [J ] , the homogeneous Eq. (3) leads to the Lüscher formula [35], det D(q) = 0, where The equivalent derivation in differential form of the variational approach by constructing a finite-volume wave function in terms of an infinite-volume wave function can be found in Ref. [23]. Using the infinite-volume wave function solutions as basis functions shows some advantages in the two-body sector: (i) the quantization condition is given in terms of physical on-shell scattering amplitudes; (ii) the relation between finite-volume solutions and infinite-volume solutions is clearly demonstrated. However, in case of multiparticle interaction, imposing physical constraints, such as the asymptotic form of the infinite-volume wave functions, and unitarity of multiparticle scattering amplitudes, means an additional difficulty on top of the task of obtaining the quantization condition itself. All the difficulties in multiparticle interaction in finite volume in fact starts with an ambitious goal at the very beginning: expressing quantization conditions in terms of infinite scattering amplitudes or wave functions, and handling both infinite-volume and finite-volume dynamics simultaneously. The multiparticle interaction in the infinite volume alone is already difficult to solve.
In this work, we explore the possibility of separating finite-volume and infinite volume dynamics; the two scenarios are connected by interaction potentials instead of scattering amplitudes. In this way, the quantization condition may be expressed in terms of periodic lattice structure and potentials alone, and is free of multiparticle dynamics in the infinite volume. The multiparticle dynamics in infinite volume may be computed separately once the information of the potentials is extracted from fits to lattice results. The proposed approach is far less ambitious, however, it provides a way of avoiding dealing with scattering amplitudes, and a much more practical formalism for just the purpose of lattice data fitting.
Before we move on to the three-body problem, we would like to demonstrate the key idea of our proposal by using the finite-volume two-body problem as a simple example. Starting the from homogeneous equation, Eq.(1), again, according to variational approach [23,25], the continuous equation, Eq.(1), may be turned into a homogenous matrix equation if the wave function is expanded in terms of certain basis functions, and then the quantization condition may be obtained by the determinant condition of this matrix equation. In principle, the discrete energy spectra do not depend on the specific choice of basis functions. In Ref. [23], the basis functions were constructed explicitly, but the specific form of finite-volume wave functions is not the focus of the present work. Alternatively, the quantization condition may be obtained from the homogeneous Lippmann-Schwinger equation directly, Eq. (1) or Eq. (3). One way of numerically solving the LS equation is by discretizing it in coordinate space, in a single box, where r i and w j denote the chosen coordinates and associated integration weights in a single box. Hence, the non-trivial solutions exist, provided that the determinant of the homogeneous equation vanishes, det D(q) = 0, where The determinant condition, Eq. (9), is equivalent to the Lüscher formula and yields the discrete energy spectra. The quantization condition, Eq. (9), may also be transferred into momentum space representation by using the Fourier transformation relation in finite volume, and the inverse Fourier transform, We find the momentum space representation of Eq. (1), where V (p) = L 3 dre −ip·r V (r). Therefore, the quantization condition may also be given by The momentum-space representation of the two-body quantization condition in Eq.(13) is in fact consistent with the two-body quantization conditions in [40,41] obtained by a different method. The approach described above may be generalized to multiparticle system as well [25]. In the following, we consider a simple 1 + 1 dimensional heavy-light threebody system with one static heavy particle interacting with two light bosons. The heavy-light system resembles an atomic system with a periodic boundary condition, and only pair-wise δ-function short-range interactions between the heavy and the light particle, and between the two light particles are considered in this work, which may be tested by a simple φ 4 -type lattice model. We remark that the short-range interactions may be described by δfunction potential and its derivatives, where δ-function term may be considered as leading order contribution, and its derivatives are considered sub-leading order effects. As for the purpose of demonstrating the feasibility of our approach in present work, only leading order contribution is considered. It turns out that the lattice model data can be fairly well described even without subleading order contributions, see Sec. IV. The more general form of interactions, including three-body potential, will be discussed in our future publications.
A. A non-relativistic heavy-light three-body system in 1 + 1 dimensions In finite volume with periodic boundary conditions the dynamics of a three-body system with two light scalar particles at x 1 and x 2 and one static heavy boson at the origin is described by the Schrödinger equation, where σ 2 = 2mE,T = ∂ 2 The potentials U L (x) = U 0 n∈Z δ(x + nL) are the short-range interactions in the heavy-light two-body sub-systems and V L (r) = V 0 n∈Z δ(r + nL) is the interaction between the light particles. L stands for the size of the onedimensional periodic box. The wave function in finite volume satisfies periodic boundary conditions, The integral representation of Eq. (14) is given by a homogeneous Lippmann-Schwinger equation, (16) or more explicitly where the periodic three-body Green's function G L 0 satisfies δ(x 1 +n 1 L)δ(x 2 +n 2 L) .
(18) The analytic expression of G L 0 is given by Typically, the representation of the lattice Green's function in Eq. (19) exhibits poor convergence, and Ewald's method has been widely used to improve it [47]. For the completeness of our presentation, rapidly converging representations of G L 0 are given in Appendix A. Using the fact that the two light particles are identical, Φ(x 1 , x 2 ) = Φ(x 2 , x 1 ), Eq. (16) can be reduced to coupled equations, By introducing a column vector, T , these coupled equations may be expressed in a simple form, where G L is a 2 × 2 matrix function, The quantization condition may be obtained by applying the variational method [23,25]. Assuming that the wave function ξ(x) may be expanded in terms of some known orthonormal basis functions Eq. (22) leads to a secular equation of familiar form, is satisfied.

Quantization condition in coordinate representation
The major goal of this work is to obtain a quantization conditions which is easy to be used to produce the discrete energy spectra and fit data of lattice model simulations. As far as the basis functions used in the variational method converge, the quantization conditions obtained by the variational method should not depend on the choice of basis ultimately. However, in reality, projecting out the matrix elements in the variational method [23,25] by multiple dimensional integration is still an computationally intense task. Fortunately, by discretizing the homogeneous Eq. (22), where ξ i = ξ(x i ), and x j and w j stand for the nodes and weights of coordinate discretization in the range of , the ξ i can be treated as basis functions for the variational method. Therefore, a quantization condition without any dependence of specific choice of basis functions is given by One last remark about the quantization condition given by Eq. (27) is that due to the singularity of , have to be regularized. A modified subtraction quadrature method [48] may be used for the regularization scheme of singularities, the details are presented in Appendix B.

Quantization condition in momentum representation
Introducing Fourier transformation in finite volume, L n (n ∈ Z), and also with the help of the identity, an equivalent representation of Eq. (22) in momentum space may be found, where Therefore, in momentum space, ξ(p) with discrete momenta, p = 2π L n (n ∈ Z), may be treated as variational basis functions, so, the quantization condition is accordingly given by where (p, p ) ∈ 2π L n, n ∈ Z.

B. Consistency check at extreme limits
We consider two extreme limits of the heavy-light system as a simple check.

Limit of U0 = 0
In the limit U 0 = 0, the two light particles decouple completely from the interaction of the static heavy boson, and the quantization condition Eq. (31), hence, is reduced to a familiar form, stands for the relative momentum of two light particles, and p = 2π L n with n ∈ Z represents the total momentum of two light particles. This is exactly what we expect for two interacting particles in a periodic box [20][21][22]46].

Limit of V0 = 0
In the limit V 0 = 0, the two light particles do not interact with each other, hence the coupled homogeneous Lippmann-Schwinger equation is reduced to a single one, where To check consistency of Eq. (33), we note that for the case of zero interaction between the two light particles, the wave function of the two light particles is simply given by the product of two single-particle wave functions, where is the solution of the finitevolume single-particle wave function in the presence of the static heavy-particle interaction, and t(q) = − U0 2q+iU0 stands for light-heavy two-body scattering amplitude. The allowed momenta of the light particle are determined by the two-body quantization condition: cot piL 2 = 2pi U0 [20,46]. Therefore, the expression of ξ 2 is given by Using the identity and also the two-body quantization condition, we find which is indeed consistent with Eq. (33).

III. A 1 + 1 DIMENSIONAL HEAVY-LIGHT SYSTEM OF LATTICE MODEL
In this section, we present a simple 1 + 1 dimensional coupled-channel φ 4 -type lattice model to simulate a heavy-light three-body system which is described in Section II. The classical action of a heavy-light bosonic system of the lattice model in two-dimensional Euclidean space is given by where x = (x 0 , x 1 ) are temporal and spatial coordinates in two-dimensional Euclidean space, respectively. The light scalar particles are represented by the φ(x 0 , x 1 ) field, and the σ(x 0 , 0) field denotes the static heavy particle at the origin. The interaction between the light particles is simply described by a φ 4 model, and the interaction between light and heavy particles is described by the coupling between the σ and φ fields at the origin: No kinetic term is needed for a static heavy σ particle.
In infinite space with an open boundary condition, the heavy-light system given in Eq. (38) is equivalent to a non-relativistic one-dimensional multiparticle interacting system with N -light-scalar plus one static heavy boson. Interactions between both light-light and heavy-light particles are described by short-range pair-wise δ-function potentials, where x 1,i refers to the spatial position of the i-th light scalar particle, and m stands for the mass of the light bosons. The coupling strengths of δ-function potential among light-light and heavy-light particles are denoted by V 0 and U 0 , respectively. When the periodic boundary condition is considered, the lattice model designed in Eq. (38) can be used to simulate a finite-volume heavylight multiparticle system that is the simple model we study in this work. A.
The lattice model action of heavy-light system The lattice action that describes heavy-light particles living in a discrete 1 + 1 rectangle lattice is obtained by replacing the continuous derivative with discrete difference in Eq. (38): where x = (x 0 , x 1 ) now refers to the discrete coordinates of the Euclidean T × L rectangle lattice sites.

B.
Monte Carlo updating algorithm for heavy-light system A combination of Hybrid Monte Carlo updating algorithm [49,50] for the φ field and the standard Metropolis-Hastings updating algorithm [51,52] for the σ field is adopted in our simulation. The φ and σ fields are updated alternately. For updating the φ field with the Hybrid Monte Carlo algorithm, an auxiliary Hamiltonian and a fictitious conjugate momentum of the φ field, the π field, is introduced, The auxiliary Hamiltonian in Eq. (41) defines the classical evolution of both π and φ fields over a fictitious time within an interval [0, τ ]: The trajectory of (φ, π) over the time interval [0, τ ] is determined by the solutions of the equations of motion in Eq. (42). The pair of (φ, π) fields and σ field are updated alternately for each sweep over entire lattice: • Updating the pair (φ, π) with the standard Hybrid Monte Carlo algorithm: (i) the trajectory begins with a random distribution of fields (φ, π) at initial fictitious time. The initial conjugate momenta, π(0), are generated according to the Gaussian probability distribution: (ii) (φ, π) are evolved over the trajectory up to a time τ according to equations of motion in Eq. (42). The equations of motion are solved by the leapfrog method [50]. The (φ, π) fields evolve along a trajectory with a fixed length τ = 8 over 100 discrete steps.
In our simulations, the lattice model parameters are chosen as: κ = 0.1275, λ φφ = 0.02 and λ σφ = 0.007. The temporal extent of the lattice is fixed at T = 100, and the spatial extent of lattice, L, varies from 10 up to 55 with an increment of 5. Two million measurements are generated for each lattice size. Plot of single light-particle spectra in the presence of the static heavy-boson potential vs. free singleparticle energy levels (red dashed curve): E f ree 1b,n (L) = cosh −1 cosh m + 1 − cos 2π L n with m = 0.163 and n = 0, 1, 2.
Lowest few orders diagrammatic representation of heavy-light two-body interactions from point of view of perturbation theory. The solid black circles stand for the static heavy particle and the light particle is represented by the solid lines. The contact interaction between light and heavy particles is denoted by a dotted line. The intersection of two solid lines represents the contact interaction between the two light particles. U and V represent the bare coupling strengths of heavy-light and light-light particles respectively.

C. Operator construction and particle spectra
In this section, some details regarding the construction of multi-particle operators and results for the heavy-light particle spectra are presented.
1. The spectra of one light scalar in the presence of a static heavy particle potential In the presence of the static heavy-boson potential, the spectra of the light scalar particle interacting with the heavy boson are extracted from the exponential decay of the correlation functions Plot of light two-particle spectra in the presence of a static heavy-boson potential (black) vs. energy levels of E 1b,n 1 + E 1b,n 2 (red), where E 1b,n are the single-particle energy spectra presented in Fig.1. The red dotted curves represent spectra of two free light particles: L ni with m = 0.163.
where the light particle propagator, φ n (x 0 ), is defined by Multiple E 1b,n are extracted for each lattice size, see Fig. 1. Because of the interaction with the static heavyboson potential, the energy spectra E 1b,n can no longer be described by a simple free-particle dispersion relation with the light particle mass as a single free parameter, i.e., they are no longer of the form E f ree 1b,n (L) = cosh −1 cosh m + 1 − cos 2π L n [53], as indicated in Fig. 1. Instead, the E 1b,n depend on two parameters: the renormalized particle mass, and the renormalized heavy-light interaction strength. The diagrammatic representations of the heavy-light two-body interactions are shown in Fig. 2. 2. The spectra of two light scalars in the presence of a static heavy particle potential The matrix element of the correlation function of the two light particles interacting with a static heavy boson is defined by where d ∈ Z is related to the center of mass momentum of the two light particles by P = 2π L d. The disconnected contribution needs to be subtracted in the CM frame (d = 0). Typically, three or four two-light-particle operators are used in our simulation for d = 0, 1, 2, The spectral decompositions of the correlation function matrices are given by where v 2b,i (0)|0 , and n labels the n-th energy eigenstate E (d) 2b,n . In order to extract excited energy states, a generalized eigenvalue method [54,55] is used, wherex 0 is a small reference time that is set to zero in this work. A simple form of λ (d) 2b,n (x 0 , 0) = e −E (d) 2b,n x0 is used in the data fitting for x 0 ∈ [0, 8] as no contamination from high excited states is observed. The energy spectra of two light particles for various lattice sizes and d are presented in Fig. 3. Examples of both one light particle and two light particles correlation functions, C(x 0 , 0), and effective masses, ln [C(x 0 , 0)/C(x 0 + 1, 0)] are given in Fig. 4 and Fig. 5, respectively.
Due to the presence of the static boson potential, the spectra E C(x 0 +1) , for one light particle (black) and two particles (blue) at L = 50 and d = 0, and corresponding fitting curves (red band). particle mass, m, heavy-light interaction strength, U 0 , and light-light interaction strength, V 0 . The diagrammatic representation of connected heavy-light three-body interactions is shown in Fig. 6. The sum of two onelight-particle energy spectra from the previous subsection III C 1, E 1b,n1 + E 1b,n2 , are also presented in Fig. 3, see red error bars, which represent the results of two light particles spectra in the absence of the interaction between the two light particles. The difference between full two light particles spectra, E 2b,n , and sum of two onelight-particle energy spectra, E 1b,n1 + E 1b,n2 , indicates the energy shift due to the non-zero interaction between the two light particles. One interesting observation is that in the limit V 0 = 0, E (d=0) 2b = E 1b (p 1 ) + E 1b (−p 1 ) (see second red energy level on left panel in Fig. 3) and E (d=2) 2b = E 1b (p 1 ) + E 1b (p 1 ) (the first red energy level on right panel in Fig. 3) are degenerate. However, the degeneracy is lifted because of non-zero interactions between the two light particles, see the full two light particle spectra in second level (black error bars) on left panel and first level on right panel in Fig. 3. We also remark that although we still use P = 2π L d with d = 0, 1, 2 to label two-light particles spectra, we do have to keep in mind that because of the presence of the static heavy boson, the total momentum of the two light particles, p 1 + p 2 , is not a conserved quantity, and p 1 + p 2 = P = 2π L d.

IV. DATA ANALYSIS
As demonstrated in Section II, although the finitevolume wave function given in Eq. (16) has no unique solution, the energy spectra of the heavy-light multiplebody system in finite volume are uniquely determined by the multiparticle interaction and the periodic lattice structure. The quantization conditions which produce discrete energy spectra in Eq. (27) and Eq. (31) do not depend on the specific form of the finite-volume wave function, nor any particular choice of basis functions in the variational approach [23,25]. As demonstrated in Section II, quantization conditions may be constructed in such the way that energy spectra are given by parameterization of potentials. Although in QCD multiparticle systems in general exhibit complicated dynamics, for some simple systems at low energies, the dynamics may be determined by a few parameters, [12,13], such as the particle masses and the interaction strength. Therefore, the main tasks of this work are to determine these fundamental parameters by the using quantization conditions to fit energy spectra of lattice model simulations. In this particular lattice model, there are three free parameters: the light particle mass, m, and the coupling strengths of the δ-potential, U 0 between the heavy and the light particles, and V 0 , between the light particles. The lightparticle mass, m, and coupling strength, U 0 , may be extracted from spectra of one light particle in the presence of the static heavy boson. The coupling strength, V 0 , may then be determined by studying the spectra of two light scalars in the presence of the static heavy boson.
At this point we face two major obstacles: (1) The lattice model presented in Section III represents a relativistic model, and the quantization conditions given by Eq. (27) and Eq. (31) are based on a non-relativistic framework. Hence, though that framework presents a clean and simple illustration of how the quantization conditions of three-body systems in a finite volume arise, the non-relativistic nature of the energy-momentum dispersion relation may not be capable to describe the lattice spectra; (2) since the simulation of the lattice model is done in discrete rather than continuous spacetime, there is the effect of finite lattice. Here, we take the discrete space into account but simply set the spacing to one. In the two-body sector, these two challenges may be remedied simply by adopting lattice dispersion relation [24], where momentum of the light-particle, p, may be deter-mined by the quantization condition In the three-body sector, one could be tempted to proceed similarly by using a dispersion relation such as [24] E (d) Unfortunately, the currently proposed approach is limited to determine σ 2 = p 2 1 + p 2 2 as apparent in the nonrelativistic quantization conditions, Eq. (27), or Eq. (31). For these reasons, it may make more sense to find a relativistic formalism which is able to incorporate the presence of discrete space on the lattice. The relativistic framework may be achieved by replacing the nonrelativistic Green's function in the Lippmann-Schwinger equation by a relativistic one. This simple prescription may be justified by the reduction of the relativistic Bethe-Salpeter equation to a relativistic Schrödinger equation under the assumption of an "instantaneous kernel function". Of course, this reduction serves only to include relativistic kinematics and is sufficient for the purpose of mapping out finite-volume effects. The details of the reduction procedure are presented in Appendix C.

A. Quantization condition in a discrete finite box
In order to take account of both the relativistic dynamics and the effect caused by discretized Euclidean spacetime, we reformulate the continuous-space finite-volume formalism of the heavy-light system presented in Section II to a relativistic finite-volume formalism defined in discrete space. Therefore, the particle coordinates are now defined only on discrete integer points in a finite box of size L, x ∈ [0, 1, · · · , L−1], where we have assumed a lattice spacing a = 1. The Fourier transform of a periodic function, f (x + nL) = f (x) is thus defined by where allowed lattice momenta are p = 2π L n, n ∈ [0, 1, · · · , L − 1], and f (p) is also a periodic function that satisfies f (p + 2πn) = f (p), see Ref. [56]. The inverse Fourier transform into the discrete finite lattice is given by Hence, the continuous dispersion relation, E 2 p = p 2 + m 2 , is replaced by (2 sinh word, we may use the following corresponding relations between continuous and discrete lattice: E p ↔ 2 sinh Ep 2 , p ↔ 2 sin p 2 , m ↔ 2 sinh m 2 , and the on-shell energymomentum relation in the discrete lattice is now given by 2 sinh E p 2 = 2 cosh m − 2 cos p .

Two-body quantization condition
For the two-body heavy-light system, the relativistic Lippmann-Schwinger equation in discrete space may be defined by The relativistic Green's function is given by  2 cosh m − 2 cos p. In the continuum limit, it then reduces to a relativistic-type Green's function, For the δ-function potential, the quantization condition for the two-body system is simply given by The difference between G D 0 (0; E) defined in Eq. (56) and its non-relativistic counterpart, G L 0 (0, p) = cot pL 2 2p , combined with the lattice energy-momentum dispersion relation, cosh E = cosh m + 1 − cos p, is shown in Fig.7.
We test the relativistic framework by using the twobody quantization condition, Eq. (58), to fit the singlelight-particle spectrum by using the light particle mass, m, and the heavy-light short-range potential coupling strength, U 0 , as fitting parameters. The results are show in Fig. 8, and the fitting parameters are: m = 0.163 ± 0.001 and U 0 = 0.07 ± 0.03.

Three-body quantization condition
For the three-body system, we need to consider a lattice version of the Lippmann-Schwinger equation defined on a discrete lattice, where (x 1 , x 2 ) ∈ [0, L − 1], and the Green's function's are defined by where 2 sinh Ei 2 = √ 2 cosh m − 2 cos p i . In the continuum limit, a relativistic three-body Green's function is obtained where E i = p 2 i + m 2 . a. Coordinate space representation of the lattice version of the quantization condition: again, by introducing a column vector, simple homogeneous equation which is defined in a discrete space is found, where G D has a similar expression as its non-relativistic counterpart in the continuum limit in Eq. (23), The coordinate-space representation of the quantization condition defined in a discrete finite box is therefore given by b. Momentum-space representation of the lattice version of the quantization condition: in momentum space, using the discrete-space Fourier transformation in Eqs. (52)(53), and also using the identity, we find where p = 2π L n, n ∈ [0, L − 1]. The kernel function, G D , is defined by where (p, p ) ∈ 2π L n, n ∈ [0, L − 1], and .
Hence, the momentum-space representation of the quantization condition in discrete space is given by det δ α,β δ p,p − G D α,β (p, p ; σ) = 0, In the non-relativistic and continuum limits, Eq. (64) and Eq. (69) are reduced to Eq. (27) and Eq. (31). In the limit U 0 = 0, in which the interaction between the light and the heavy static particles vanishes the threebody quantization condition is, thus, reduced to a simple form, (70) The non-relativistic counterpart of the above equation is given by Eq. (32), hence p = 2π L (0, · · · , L − 1) may be related to the total momentum of the light-particle system. The comparison of 1 L p = 2π L n n ∈[0,L−1] G D 0 (p , −p ; E) in the CM frame of the two light particles and its nonrelativistic counterpart, 1 2σ , combined with lattice dispersion relation, cosh E 2 = cosh m + 1 − cos σ √ 2 , is shown in Fig. 9. non-relativistic Green's function (dashed red), using the lattice energy-momentum dispersion relation, cosh E 2 = cosh m + 1 − cos σ √ 2 . The lattice size and mass of the light particle are: L = 20 and m = 0.163.
Using the light particle mass of m = 0.163 ± 0.001, and heavy-light short-range potential coupling strength, U 0 = 0.07 ± 0.03 from the fit to the two-body spectra as input, the remaining parameter, i.e., the coupling strength of the short-range potential between the light particles, V 0 , is determined by fitting the two-lightparticle spectra using the three-body quantization condition, either Eq. (64) or Eq. (69). The fitting results are show in Fig. 10. The coupling strength is determined as V 0 = 0.43 ± 0.03 through the fit. At last, we would like to remark that due to the finite lattice spacing a, the light particle mass and coupling strengths extracted from the lattice simulation are in principle a-dependent as well. The physical value of these quantities in continuum limit may be extrapolated by using multiple lattice spacing simulation results combined with the renormalization group method that yields the explicit lattice spacing dependence of renormalized physical quantities based on perturbation theory. In principle, the multi-particle dynamics in finite volume and in infinite volume are related though same set of interaction parameters. Once all the interaction parameters are extracted from lattice results, the physical scattering amplitudes in infinite volume may be determined and computed though standard procedures, such as Faddeev's approach. The example of non-relativistic Faddeev's approach is listed in Appendix D.

V. SUMMARY
In summary, a 1 + 1 dimensional heavy-light threebody system is simulated by a coupled-channel φ 4 -like lattice model. An improved variational approach for the finite-volume multiparticle problem is proposed in this work. The advantage is that the quantization conditions are given in terms of the periodic lattice structure and interaction potentials instead of explicit finite-volume scattering amplitudes. This opens up the possibility of finding a more practical formalism which can be easily used for data fitting. The interaction potentials are varied to match the energy spectrum; once the parameters of the potentials are extracted by data fitting, the dynamics of the entire system is determined. The multiparticle dynamics in infinite volume hence may be computed separately. The periodic Green's function, G L 0 , defined in Eq. (18), has the analytic expression where x = (x 1 , x 2 ) and n = (n 1 , n 2 ). Equivalently, it can also be expressed as Unfortunately, both expressions in Eq. (A1) and Eq. (A2) suffer poor convergence, especially for σ values on the real axis.
In this section, we present two equivalent fastconverging expressions of G L 0 . (i) First of all we use the identity where η is an arbitrary parameter. For real σ values, the integration path has to be defined carefully in the complex plane to warrant the convergence of integration [47], see Fig. 11. Hence, we can rewrite G L 0 in Eq. (A1) as In this expression, the first three terms are all wellbehaved for real σ values. For the last term in Eq. (A4), using the identity and also applying Poisson summation, we find Putting everything together, we find The integration in the last term in Eq. (A7) is normally highly suppressed for L > 10 and η > 1. Therefore, with a modest choice of η ∼ 3, the last piece in Eq. (A7) can be safely ignored in the numerical evaluation.
(ii) Secondly, using the relation, the second fast convergent expression of G L 0 can be obtained from Eq.(A2), subtraction quadrature method [48] for the regularization scheme of singularities. First of all, we can subtract the same term on both side of Eq. (22), . (B2) Using relations, we can discretize Eq.(B1) to Hence, all the singular terms cancel out as x j = x i on the left hand side above equation, discretized Eq.(B5) now is well-defined. Reorganizing Eq.(B5), we find where

Reduction of two-body system
Let's consider the general form of Bethe-Salpeter equation [57] for one heavy and one light scalar particles bound state system, where p 1 = (p 1,0 , p 1 ) and p 2 = (p 2,0 , p 2 ) are the four momenta of two particles, m and M represent mass of light and heavy particles respectively. q = M p1−mp2 m+M denotes the relative four momentum of two particles. ψ(q) is the Bethe-Salpeter wave function, and I(k) represent ladder approximation interaction kernel function. Using approximation of "instantaneous interaction kernel", I(k) = I(k), in which interaction kernel has no dependence on the energy components of four momenta, and also introducing the Schrödinger equation wave function by ψ(q) = dq0 2π ψ BS (q), we obtain, (C2) Next, using identity (C3) where E 1 = p 2 1 + m 2 , E 2 = p 2 2 + M 2 and E is total energy of two-particle system. In the heavy static limit, E 2 → M → ∞ and p 1 → q, we also shift total energy by M to subtract the heavy mass, E → E − M , and introducing potential by V = 1 2M I, then we find (C4) In non-relativistic limit, Eq.(C4) hence yield familiar Lippmann-Schwinger equation, Reduction of three-body system Now, let's consider one heavy and two light scalar particles bound state system, the Bethe-Salpeter equation is given by where p i = (p i,0 , p i ) are the four momenta of three particles, particle-1 and -2 are labeled as light particles and particle-3 denotes heavy particle. The relative momenta are introduced by q 1 = (M +m)p1−m(p2+p3) 2m+M and q 2 = (M +m)p2−m(p1+p3) 2m+M . Again, assuming "instantaneous interaction kernel", I(k 1 , k 2 ) = I(k 1 , k 2 ), and introducing Schrödinger equation wave function, ψ(q 1 , q 2 ) = dq1,0 2π dq2,0 2π ψ BS (q 1 , q 2 ), hence, we get The integration over propagators can be carried out by where E 1,2 = p 2 1,2 + m 2 and E 3 = p 2 3 + M 2 . When only pair-wise interactions are considered, the three-body interaction potential, I, may be given by sum of interactions among all pairs, ). (C9) The relativistic kinematic factors p i |p i = 2E i (2π) 3 δ(p i − p i ) emerge when i-th particle is free propagating and not involved in the interaction. At the limit of static heavy particle, M → ∞, p i → q i , and define potentials by U = 1 2M I U and V = 1 2m I V , we find ψ(q 1 , q 2 ) = 1 If all interactions are turned off except interaction between pair (13), above equation is thus reduced to ψ(q 1 , q 2 ) = 1 2E 1 which is just a relativistic two-body LS equation in Eq.(C4) with second particle as a spectator.
In non-relativistic limit, Eq.(C10) yield a familiar form, Appendix D: Solutions of heavy-light three-body system in infinite space In infinite space, the dynamics of a non-relativistic heavy-light three-body system is described by Schrödinger equation, where σ 2 = 2mE,T = ∇ 2 x1 + ∇ 2 x2 , and r = x 1 − x 2 . For scattering solutions, the wave function has the following form [58,59], Ψ = Ψ (0) + 3 γ=1 Ψ (γ) , where Ψ (0) stands for the incoming free wave, and Ψ (γ) satisfy coupled equations, Ψ (i) (x 1 , x 2 ) = dx 1 dx 2 G (i) (x 1 , x 2 ; x 1 , x 2 ; σ) and Ψ (3) (x 1 , x 2 ) = dx 1 dx 2 G (3) (x 1 , x 2 ; x 1 , x 2 ; σ) where the Green's functions are the solutions of equations, and where q = p1−p2 2 and P = p 1 + p 2 . Usually, singular terms, T (0) 's, are not convenient for numerical computation, so it is better to introduce g-amplitude by a shift, where Q (i) 's satisfy coupled equations, The analytic expressions of Q (i) 's are given by , (i = j) = 1, 2. (D17) The g amplitudes are the solutions of coupled equations that are free of δ-function type singularities, and g V (p) = g (0) where In terms of g amplitudes, the wave function of threebody heavy-light system is thus given by where ψ p (x) = e ipx + it U ( p 2 )e i √ p|x| is heavy-light twobody wave function. The first term on the right-hand side of Eq.(D21), ψ p1 (x 1 )ψ p2 (x 2 ), stands for the solution of heavy-light quark system at the limit of zero interaction between two light particles, V 0 = 0. The second term represent the disconnected contribution of two light particles interacting while heavy particle play the role of spectator. The last two terms in Eq.(D21) produce diffraction effect that are generated by the rescattering effect among all the pairs. Integral equations of g amplitudes can be solved with standard matrix inversion method, see Fig.12 for examples of numerical solutions of Eq.(D18-D19).