Two- and three-body interactions in $\varphi^4$ theory from lattice simulations

We calculate the one-, two- and three-particle energy levels for different lattice volumes in the complex $\varphi^4$ theory on the lattice. We argue that the exponentially suppressed finite-volume corrections for the two- and three-particle energy shifts can be reduced significantly by using the single particle mass, which includes the finite-size effects. We show numerically that, for a set of bare parameters, corresponding to the weak repulsive interaction, one can reliably extract the two- and three-particle energy shifts. From those, we extract the scattering length, the effective range and the effective three-body coupling. We show that the parameters, extracted from the two- and three-particle energy shifts, are consistent. Moreover, the effective three-body coupling is significantly different from zero.


Introduction
In the last decade, the use of the (multi-channel) Lüscher equation [1][2][3][4][5][6][7][8][9][10][11] has become a standard tool for the extraction of the infinite-volume two-particle scattering phase shifts and inelasticities from the energy spectrum on the lattice [12][13][14][15][16][17][18][19][20][21][22]. Further, the Lellouch-Lüscher formalism, as well as its generalizations to the case of the coupled two-particle channels, have been successfully used for the measurement of various matrix elements with two-particle final states, as well as for the calculation of the resonance form factors [7,[23][24][25][26][27][28][29]. Very recently, the focus in the development of the formalism has been shifted from two-to three-particle systems, where a significant amount of work has been done in the course of last few years . In particular, the formalism of Refs. [46,47,49] provides a systematic fit strategy to the data, even if the latter consists of only few data points with sizable errors (a conceptually similar approach was proposed in Ref. [50]). The output of this fit determines the low-energy a e-mail: fernando.romero@uv.es constants in the three-particle sector (three-particle force), which can be further used, as an input, in the infinite-volume scattering equations, in order to determine all S-matrix elements in the three-particle sector. It is not immediately clear, however, whether this ambitious program can be realized in practice or, in other words, whether the three-body force can be reliably extracted from the lattice data. The present study intends to explore this gap and test the viability of the theoretical framework.
The strategy laid out in Refs. [46,47,49] could be briefly summarized as follows: the interactions in the effective Lagrangian, describing three particles, are parameterized by two sets of couplings. The first set of couplings describes the low-energy interactions in the two-particle sub-systems. These can be related to the usual effective range expansion parameters: the scattering length a, the effective range r , and so on. The second set corresponds to the tower of the threebody local couplings, which are multiplying the operators with an increasing number of space derivatives acting on the fields. At lowest order, there is just a single constant describing the contact non-derivative interaction of three particles at one point. If the particle-dimer picture is used, these two sets of couplings are uniquely translated into the sets describing the dimer-two-particle vertex and the contact particle-dimer interactions. The details are given in Refs. [46,47].
Next, note that the couplings in the two-body sector can be independently fitted to the two-particle phase shift, extracted from the lattice data by using the Lüscher method. On the contrary, the finite-volume spectrum in the three-particle sector, which is determined by the quantization condition derived in Ref. [47], depends on both sets of couplings, with the former already fixed from the fit to the two-particle scattering data. Thus, the fit to the three-particle energy levels allows one to determine the first few couplings in the three-particle (particle-dimer) sector. Using the values of the couplings, determined in this way, one may then calculate all infinite volume S-matrix elements.
Note also that the presence of the three-body force in the equations (both in infinite and finite volume) is crucial. These equations contain the UV cutoff and, without the threebody force, one could not ensure UV cutoff-independence of the observable quantities, obtained by solving the equations. Namely, the coupling constants should show the log-periodic dependence on the cutoff (see, e.g., [52,53]). 1 The aim of our investigation is, 2 in particular, to answer to the following questions: (i) Is the theoretical approach, proposed in Refs. [46,47,49], practically viable? Could one unambiguously extract the effect from the genuinely three-body (particledimer) coupling(s) from data with the present computational ressources? (ii) What is the best fit strategy? At which volumes should the data be taken to be most sensitive to the three-particle coupling(s)? (iii) Does one have all finite-volume artifacts under sufficient control in order to extract the three-body coupling(s) with a reasonable systematic error?
In order to answer these questions we have used our framework from the beginning to the end to analyze the data in the toy model described by the ϕ 4 Lagrangian. 3 On the one hand, this model is chosen for the illustration of the fact that the three-body calculations can be systematically carried out even without committing extraordinary computational capabilities. On the other hand, we use this model as an example to discuss the general issues, which were outlined above. For simplicity, we will from now on assume that the threeparticle force is described by the single non-derivative coupling. It will be seen, a posteriori, that the accuracy of lattice data is not yet sufficient for extracting higher-order couplings as well.
It is important to realize that the role of the three-body force in the following situations is physically very different: (i) There exists a three-particle bound state: in nature this scenario is realized e.g. in a particular three-nucleon systems (triton). In this case, one can directly extract the pertinent three-body coupling by extrapolating the measured bound-state energy to infinite volume. (ii) There is no three-particle bound state: in nature this corresponds to e.g. the final-state interactions in the K → 3π decay. In large volumes, the free three-particle energy levels are displaced (such a case was considered, e.g., in Refs. [37,38,45,[58][59][60][61]). The displacement contains the three-particle coupling only at N 3 LO in the perturbative expansion in 1/L, where L is the size of the box and, hence, it is legitimate to ask, whether this coupling can be extracted from data at any precision. Further, making L smaller in order to enhance the relative contribution of the three-particle force, one is necessarily confronted with the question about the exponentially suppressed finite volume effects. Does one control these effects to a sufficient precision, to ensure a clean extraction of the three-particle coupling?
We stress here once more that in nature one encounters both scenarios, so it is useful to address them both in the model study. Whereas the extraction of the three-particle coupling in scenario (i) is relatively straightforward, the scenario (ii) should be given more scrutiny.
In particular, we would like to point out the issue of consistency of the infinite-volume and the finite-volume descriptions of the same system in the framework we are considering. Suppose that scenario (ii) is realized. Then, the three-body coupling appears only at order L −6 in the perturbative expansion of the ground-state energy level, which starts at O(L −3 ) (to be more precise, a certain cutoff-independent combination of this coupling and the UV cutoff should appear in the expression of the finite-volume energy). Owing to this fact, the extracted value of this coupling will potentially have large errors. This could lead to problems, if the same coupling would emerge at the leading order in the infinite-volume scattering amplitude. It can be seen, however, that this is not the case -the pertinent contribution to the scattering amplitude is suppressed in perturbation theory.
Finally, note that the quantization condition given in Ref. [47] determines the energy levels implicitly as the solutions of the secular equation. Perturbative expansion in the vicinity of the unperturbed energy levels (both for the ground state and the excited states) is possible [62]. As expected, after this expansion, the results of Refs. [37,38,45,[58][59][60][61] for the ground state are reproduced 4 -there is nothing more about the quantization condition in this case. Therefore, in order to fit the three-particle energy levels, we shall further use the expressions given in Ref. [45]. A detailed comparison of the threshold expansion given in Ref. [45] and the expansion emerging in our approach, es well at the explicit result for the first excited three-particle state is forthcoming [62].

Description of the model
We now turn to the description of the model which will be used in numerical lattice simulations. We choose the complex ϕ 4 theory, because the charge quantum number protects the mixing of one-, two-and three-particle states. The continuum Euclidean Lagrangian reads and the discretization is performed as follows: and this way, the discretized action reads with The interpolating operators used for the one-, two-and threeparticle states are: respectively. In principle, one must subtract the vacuum expectation value (vev) to an operator before constructing the correlation function. However, the transformation is a symmetry of the action and the operators in Eq. (7) have a vanishing vev: This way, the correlation functions for the one-particle state are built as follows: and, similarly, for the two-and three-particle states. That is, we extract the one-, two-and three-particle energies in the center-of-mass (CM) frame.
The energies can be calculated by fitting the correlation functions to their theoretical form, which can be determined by using the Transfer Matrix formalism: and calculating explicitly in terms of all possible states: This way, the different correlation functions look like: The additional terms in C 2 and C 3 come from thermal pollution, i.e., antiparticles propagating backwards in time and crossing the periodic temporal boundary. They are relevant, since the matrix elements A 1→1 ∝ ϕ † |Ô 2ϕ |ϕ and A 1→2 ∝ ϕ † |Ô 3ϕ |ϕϕ are different from zero and, in general, they are of the same order of magnitude as A 2 and A 3 , respectively. In contrast to C 2 , the pollution term in C 3 is time-dependent. Both of these additional terms vanish in the limit of T → ∞ but, for finite T , they have to be taken into account in the analysis. For C 2 , this can be done by using the discrete derivative of the correlation function, the so-called shifted correlation function: This has also advantage in terms of less correlation among the different time slices [63]. Using Eq. (16) turns the functional dependence on t in C 2 (t) into sinh. For C 3 , we take the timedependent pollution into account explicitly in our fitting. In doing so, we are using the single particle mass M and the two particle energy E 2 as an input, determined from C 1 and C 2 , respectively. Our ensembles are generated using the Metropolis algorithm with simultaneous updates in the even and odd parts of the lattice. We choose m 2 0 = −4.9 and λ c = 10.0 implying λ = 0.253308 and κ = 0.159156. We calculate the spectra in 15 different volumes, L = [5,18] and L = 20 with T = 24 and two volumes L = 14, 24 with T = 48. The number of independent configurations is always in the range 7500-100,000 and the errors are calculated through jackknife resampling. A summary of all ensembles can be found in Table 2. For the reasons given in the previous paragraph, we fit the corresponding sinh functional form to our data for the correlation functionsC 1 ,C 2 andC 3 . The excited states are expected to be strongly suppressed because we use charged operators and because of the small value of the renormalized coupling λ c . In practice, there seems to be no sign of excited states in the correlation functions, and thus the fits are performed to all time slices. This can be seen in the effective mass (m eff ) plots compiled in Appendix B, where m eff is defined through: For C 2 , the relation looks identical. For C 3 , the functional form is more complicated, but m eff can be defined likewise.

The spectrum in finite volume
In the ϕ 4 theory in the symmetric phase, which is studied in the present paper, the vertices with the odd number of the field ϕ are barred. This situation resembles the one in chiral perturbation theory (ChPT) and we can straightforwardly use the perturbative expression, obtained in Ref. [64] (see also Refs. [65,66]), to describe the volume-dependence of the single-particle mass. This expression has the following form where M = M ∞ and K 1 (z) denotes the modified Bessel function of first kind. Taking into account the asymptotic expression of the Bessel function for a large value of the argument, we get The threshold expansion of the two-and three-particle energies can be taken from the literature. For convenience, we generally stick to the notations used Ref. [45] (we use lattice units in the formulae below). The two-body energy shift E 2 reads: where a and r denote the two-body scattering length and the effective radius, respectively, and c 1 , c 2 , c 3 are known numerical constants. The last term is the relativistic correction, which vanishes at M → ∞.
The three-body energy shift of the ground-state level E 3 takes the form: where d 1 , d 2 , d 3 are the known numerical constants and the parameter D is a sum of different terms including the three-body coupling constant we are looking for (or, the divergence-free three-body scattering amplitude M 3,thr in the notations of Ref. [45]). The individual terms in this sum (including the three-body coupling) depend on the ultraviolet cutoff -only the sum of all terms in D is cutoff-independent. 5 In addition, the quantity D depends on the choice of the scale in the logarithm, which is present in Eq. (21) -there, it is chosen to be equal to M, which is the only available natural scale in the theory. Choosing any other scale around this natural scale leads to an additive contribution to D. Consequently, our goal can be re-formulated as follows: we want to show that D is clearly non-zero beyond the statistical uncertainty, when the scale in the logarithm is of order M. This fact would be interpreted as a detection of the contribution of the threeparticle force, which also emerges at O(L −6 ). Any further refinement of the argument seems not to be possible, because the individual contributions are scale-and cutoff-dependent. Before fitting the above formulae for the energy shift to the lattice data, the following important question is in order.
As we see, in order to reliably extract the coefficient of order L −6 , we have to go to not so large L. In this case, the exponentially suppressed corrections in L, which were neglected in Eqs. (20), (21), might become non-negligible (and, as we shall see, they really do). Then, what is the systematic error imposed on the extracted value of the coefficient of order L −6 by neglecting such terms? We shall show that one may treat the leading exponential corrections in a relatively simple fashion, greatly improving the accuracy of the Eqs. (20) and (21).
In the infinite volume, the two-and three-particle thresholds are located at the energies 2M and 3M, respectively. In a finite volume, this shifts to 2M L and 3M L . This exponentially suppressed shift, albeit not large, may become statistically relevant, if low L values are included in the fit. Below we shall argue that, up to next-to-leading order in perturbation theory, using E 2 = E 2 (L) − 2M L and E 3 = E 3 (L) − 3M L already captures the bulk of these exponentially suppressed corrections -the remaining terms should be suppressed by two powers of L and by one power of the coupling constant λ c with respect to the leading correction. 6 In order to study the exponentially suppressed terms, one has to abandon the non-relativistic effective field theory, which provides a very comfortable framework to derive equations like (20) or (21), and turn to the relativistic description of the system. In the following, we work in Minkowski space and consider the two-particle case first. It is convenient to start with the Bethe-Salpeter equation for the two-body scattering amplitude T in the center-of-mass (CM) frame. In the infinite volume, the equation takes the form: 7 Here, the 4-momenta p, k, p denote the relative momenta of the two-particle system with total momentum P = (E, 0). Further, K ( p, p ) is the kernel of the Bethe-Salpeter equation (the sum of all two-particle irreducible diagrams), and is the product of two dressed one-particle propagators, including all self-energy insertions. To the order we are working, however, these selfenergy insertions do not contribute and G can be replaced by the free relativistic propagator.
In a finite volume with periodic boundary conditions, the counterpart of Eq. (22) is written as: where L stands for the spatial size of the (cubic) box (the temporal extent of the box is assumed infinite). The threemomenta p, p are discretized as well. Further, in analogy with Ref. [67], we single out the positive-energy contribution 6 It is conceivable that the argument can be extended to all orders in perturbation theory -at least, we were not able to find an example of a diagram in higher orders that would invalidate it. A full proof, however, seems to be a rather challenging affair. Since we are dealing here with the perturbative case, where the coupling constant λ c is small, we were tempted to restrict ourselves to a statement, which is valid up to the next-to-leading order in λ c only. 7 We implicitly assume that the ultraviolet divergences are cured, e.g., by introducing some cutoff in the integrals. These divergences are not relevant in the discussion of the finite-volume effects. in the two-particle propagator by writing In the above equation,Ĝ 2L (k) contains in general the contributions from the negative-energy states (anti-particles), as well as the contributions coming from the many-particle intermediate states in the spectral representation. Defining further and we get: If one now singles out the term with k = 0, it is possible to definê Then, it can be shown that there exists an algebraic relation between T L andT L : From the above expression, it is seen that the exact energy shift of the ground state is determined from the equation Note that the quantityT L (0, 0) itself depends on the energy E = E 2 (L), so the solution of the above equation can be written down perturbatively, expanding order by order in the energy shift (E − 2M L ). Let us now find the solution of this equation up to second order in the coupling constant λ c . We shall be interested in the exponentially suppressed terms only, since the power-suppressed terms will obviously coincide with those in the Lüscher equation. At lowest order in λ c , we get T L (0, 0) = −24λ c = const, so there are no exponentially suppressed terms at all. These appear first at order λ 2 c . The Bethe-Salpeter kernel at this order is given by where K t,u L are the t, u-channel one-loop diagrams of the type shown in Fig. 1. The Feynman "integral" in a finite volume, which corresponds to this diagram, is given by 8 This result can be most easily obtained by noting that Eqs. (18), (19) are nothing but the contribution from the tadpole graph to the one-particle self energy. Then, the leading asymptotic behavior in Eq. (32) where Using Poisson's summation formula, one can single out the leading exponential correction to the infinite-volume result: 8 We have replaced M L by M everywhere, since the difference is exponentially suppressed and does not matter at the order we are working.
where the ellipses stand for the more suppressed terms. Shifting now the integration contour in k 1 into the complex plane: k 1 → k 1 + i M (we again replace M L → M), and rescaling the integration variables according to k 1 → k 1 /L , k 2,3 → k 2,3 / √ L, we obtain: Here, w(k) = √ M 2 + k 2 and k ⊥ = (k 2 , k 3 ). It is seen that the leading exponentially suppressed terms here have an extra suppression factor L −1/2 , as compared to the t, u-channel contributions.
The calculation of I L proceeds analogously, with the difference that the denominator is now singular -so, one has to carry out subtractions in the numerator first, in order to be able to use Poisson's formula. The final result is given below: where q 2 0 = E 2 /4 − M 2 . The last two terms in the above expression give the familiar power-suppressed contributions to the Lüscher formula, which are already contained in Eq. (20) (in particular, the last term is responsible for the relativistic correction). Putting things together, we finally conclude that the leading exponentially suppressed corrections in T L (0, 0)/L 3 are of order λ 2 c e −M L /L 7/2 , whereas in the mass M L these contribute at order λ c e −M L /L 3/2 . Thus, calculating the shift E 2 = E 2 (L) − 2M L instead of E 2 (L) − 2M one gains the suppression factor λ c /L 2 for the exponential contributions at no extra cost. Considering the three-particle energy we remark that the above result for the two-particle energy shift could be interpreted in terms of the non-relativistic effective Lagrangians, with the coupling constants having exponentially suppressed contributions. to the data in the fit range [4,24] a L = a + O(e −M L /L 1/2 ), and so on. One can use these volume-dependent constants in Eq. (21), in order to find the leading exponential terms. Further, in this language, the leading exponential correction to the three-particle coupling constant arises from the diagram, shown in Fig. 2. In analogy with Eq. (32), it is straightforward to derive that the leading exponential corrections coming from this diagram are suppressed by the factor e −M L L 1/2 . Using this information in Eq. (21), we arrive at the same conclusion as in the two-particle case, namely, that replacing E 3 (L) − 3M by E 3 (L) − 3M L , one gains the suppression factor λ c /L 2 for the exponential contributions.

Numerical results
We now address the analysis of the lattice spectrum in our model, which is shown in Table 2 in Appendix A. Moreover, as an illustration we also show in Appendix B the plots for the effective mass and the correlation function for the case L = 18. First, we study the volume dependence of the single particle mass, which is given by Eq. (18). The results are rangeindependent and seem to describe the data very well. This is immediately seen in Fig. 3, which displays the fit of Eq. (18) to the data on M L in the whole range of volumes from L min = 4 to L max = 24. The fit yields indicating a very good description of the data. The best fit parameters for the fits with other fit ranges are compiled in Appendix D. Next, we consider the two-particle energy shift E 2 (L), for which data from L = 4 to L = 24 are available. We fit Eq. (20) to the data, varying the lower end of the fit range in L and the number of the fit parameters. The upper end in the fit range is kept fixed at L max = 24. For large L min , one expects that only the scattering length is important and that in Eq. (20) one can neglect all terms of order L −6 or smaller. Adding the effective range r as a fit parameter, it should be possible to include smaller values of L in the fit, while reproducing the results of the fit with only the scattering length as a free parameter.
In the left panel of Fig. 4, we show an example of a fit of Eq. (20) to the data including all terms up to order L −6 , i.e. with a and r as fit parameters. L min = 9 was chosen for this particular fit. In the right panel of that figure, we show the p-value of the fit of Eq. (20) to the data as a function of L min . Red circles correspond to the one-parameter fit up to order L −5 , and blue squares to the two-parameter fit up to order L −6 . For the former, we observe the p-values around 0.5 from L min = 11 and larger, for the latter from around L min = 8.
From these L min values on, we also observe stable values for the fit parameters a and r within errors. This is shown in Fig. 5, where we plot a in the left panel and r in the right panel, both as functions of L min . In addition, as one infers from the left panel, the values for a agree within errors between the single-and the two-parameter fits from the aforementioned L min -values onward. As a result, we quote here the best fit parameters from the two-parameter fit with L min = 9, which can be found in Table 1.
Further, each data point for E 2 can be translated into one phase shift at a certain value of the scattering momentum. For that, we calculate the S-wave phase shift as [1] cot δ = Z 00 (1, where Z 00 (1, q 2 ) is the generalized Lüscher zeta-function, q = Lk 2π and 9 E 2 = 2 k 2 + M 2 L . Note that we used here the mass M L , taking into account the discussion in Sect. 3. The results can be seen from Table 3 in Appendix C. The phase shift obeys the inequality |δ| < 2 • , guaranteeing that with the current model parameters we are in the perturbative regime. The smallness of the phase shift also supports the validity of the arguments in Sect. 3. In Fig. 6, the phase shift is shown as a function of the scattering momentum k. We are now able to check, whether the parameters obtained from the fit to E 2 also describe the phase shift data properly. For this purpose, we use the effective range expansion: with the best fit parameters compiled in Table 1. In Fig. 6, we show the phase shift (in degrees), determined by using Eq. (39), as a function of the scattering momentum k (in lattice units). In addition, we show δ(k), determined from Eq. (40) by using the best fit parameters from the fit to E 2 (the solid line in this figure). It can be observed that both data analysis methods match nicely. Turning to E 3 now, we have checked first that the data for the ratio is close to 3, as one expects. This is shown in the left panel of Fig. 7. Next we fit Eq. (21) to the data for E 3 , including only the scattering length a as a single fit parameter. Therefore, we include only the terms up to the order L −5 in the fit. If everything is consistent, this should allow us to reproduce the value of the scattering length, determined by the fits to E 2 for sufficiently large values of L min . The result of this exercise is shown in the left panel of Fig. 7 and the right panel of Fig. 8. We plot a and the p-value of the fit, both as functions of L min . We observe reasonable p-values from L min = 11 onward and the best fit values for a lie between 0.4 and 0.5. Within errors, this is in agreement with the value quoted in Table 1, which was obtained from the fit to E 2 . We hence conclude that, for sufficiently large values of L, the data for E 3 and E 2 are reasonably well described by the same value of the scattering length a. We cannot perform the same test, including also r , because at that order also the three-body term contributes.
As a next step, we will use the results for a and r from the fit of Eq. (20) to the data for E 2 as priors for the fit of Eq. (21) to the data for E 3 , including terms up to the order L −6 . To this end, we introduce an augmented χ 2 function with χ 2 defined as usual, P r and P a being the priors for a and r , respectively, and a and r denoting the corresponding standard errors. In the fit, we minimize the augmented χ 2 aug . The actual values for P a , P r , a and r are the best fit parameters with errors, which can be found in Table 1. The goal here is to see whether we get a significant contribution from the three-body term in Eq. (21) or not. An example of a fit with L min = 9 and L max = 24 is shown in the right panel of Fig. 8.
The scattering length a and the effective range r are highly constrained by the prior included in the fit. This can be observed from Fig. 9, where a (left panel) and r (right panel) are plotted as functions of L min . The fit to E 3 appears not to be sensitive to the effective range r , with the prior value always reproduced. This is expected, because the term proportional to r interferes with the term proportional to D. Thus, the fit will always choose the prior value for r to minimize the augmented χ 2 . For a, we again observe that, for sufficiently large L min , the value from the fit to E 2 is reproduced within errors. From this, one may infer that L min ≥ 9 should be chosen.
This conclusion is also supported by the p value, shown in the left panel of Fig. 10 as a function of L min . For L min = 8 and 9 we observe reasonable p values. For L min = 10 or larger the p-values become close to 1, indicating at too many fit parameters for too few data points. Finally, in the right panel of Fig. 10, we show the three-body parameter D as a function of L min . In the relevant region of L min we observe a plateau within errors at a value significantly different from zero. Actually, D is different from zero to many σ for all considered values of L min , in particular 4σ for L min = 9. Final results of the fit to E 3 are shown in Table 1.
Eventually, we can use the parameters from Table 1 to compute the ratio E 3 / E 2 , using Eqs. (20) and (21), both to the order L −6 . The result is shown as a solid line in the left panel of Fig. 7. The agreement is satisfactory within error bars.
We close this section with two remarks: first, we have also included an effective term of order L −7 in the fits and observed that our results are stable under such a change. Second, a direct fit to the data for E 3 / E 2 appears not to be sensitive to the fit parameters we are interested in.   Table 1, using Eqs. (20) and (21). Right: scattering length a as a function of L min determined fits of Eq. (21) up to order L −5 to our data for E 3 , i.e. with a the only fit parameter

Summary and outlook
In this paper, we have investigated the two-and three-particle interactions in complex ϕ 4 theory on the lattice. In particular, we have carried out a theoretical study of the exponential finite-volume corrections to the one-, two-and three-particle energy levels. To next-to-leading order in perturbation theory it was shown that, using the finite-volume particle mass M L instead of the infinite-volume one M in the definition of the two-and three-particle threshold energies allows one to significantly reduce the exponential corrections to the Lüschertype formulae for the two-and three-particle ground-state energy shifts. For example, in case of the model considered, these corrections get a suppression factor λ c /L 2 , when M L is used. A proof of this statement to all orders in perturbation theory seems challenging 10 but worth trying, since it could be interesting for many lattice QCD applications. We plan to take up this challenge in our future publications. The big advantage of ϕ 4 theory as compared to QCD lies in a much reduced computational complexity. This fact allowed us to study a large set of ensembles with different finite volumes at fixed values of bare parameters. At the first stage, we study the finite-volume corrections to the single-particle  21) to the data for E 3 . Both a and r were included in the fit with priors from the fit to E 2 (the values can be found in Table 1) mass on these ensembles. We show that these corrections are well described by the corresponding expression known from the effective field theory. Next, for each of the aforementioned ensembles, we determine the two-particle and threeparticle energy shifts in a finite volume, E 2 and E 3 . Then, we extract the scattering length a and the effective range r from the L-dependence of E 2 by applying the Lüscher formalism. We also show that, depending on the range of values of L, used in the fit, either only a or a and r together can be determined, with consistent results for a in both cases. As a next step, we show that we can extract the scattering length for sufficiently large values of L from E 3 as well, and its value turns out to be consistent with the value extracted from the fit to E 2 . Hence, we conclude that our data for E 2 and E 3 are consistent. Now, we follow the idea of Refs. [46,47] and use a (and r ), determined from E 2 , as an input for the analysis of E 3 . In this case, the (higher order) formula for E 3 includes the three-body effective coupling as well. On the basis of the analysis of data, taken at many different volumes, we find a statistically significant threebody contribution D, which is definitely away from zero. This is a central result of our work. 11 11 One should not be confused by the fact that the three-body force emerges irrespectively to the absence of the ϕ 6 vertex in the Lagrangian. The notion of this force cannot be separated from the framework used to describe the energy levels in a finite box -the non-relativistic effective field theory, in our case. In this theory, such a force is present from the beginning. In other words, we could simply state that our analysis is sensitive to the terms of order L −6 in the three-particle energy, where the three-particle force starts to contribute.
In this simple theory, there are still many issues to be explored. For example, one could study the multiparticle states with more than three particles and try to figure out, how their energy depends on the volume. Moreover, one may include excited three-particle states in the rest frame or moving frames for a more sophisticated study of the three-body coupling. In addition, it would be interesting to investigate ϕ 4 theory at parameter values where a the three-body bound state is present, since the volume-dependence of such an energy level is qualitatively different from the one seen it this work. For this, one would have to further explore the parameter space of the theory.    Table 3.