Relativistic $N$-particle energy shift in finite volume

We present a general method for deriving the energy shift of an interacting system of $N$ particles in finite volume in terms of the threshold parameters of the two- and three-particle relativistic scattering amplitudes. To this end, we use the nonrelativistic effective field theory (NREFT), and match the pertinent low-energy constants to the scattering amplitudes. Relativistic corrections are explicitly included up to a given order in the $1/L$ expansion. We apply this method to obtain the ground state of $N$ particles, and the first excited state of two and three particles. We use these expressions to analyze the $N$-particle energy shift in the complex $\varphi^4$ theory.


Introduction
Recent years have witnessed a substantially increased interest in the study of many-body dynamics (three particles and more) on the lattice. On the one hand, available computing resources already enable one to carry out calculations of the three-meson spectrum in QCD at a reasonable accuracy, so that results of lattice measurements can be used for the extraction of parameters of many-body interactions in infinite volume [1][2][3][4][5][6][7]. On the other hand, there has been substantial progress in the development of the framework that enables one to analyze the lattice data.
These developments have followed different patterns. First, Rayleigh-Schrödinger perturbation theory has been used to calculate the ground-state energy shift of N identical bosons in a finite box of size L in a systematic expansion in powers of 1/L (modulo logarithms). This problem has been attracting attention for decades already (see, e.g. [8][9][10]) but in the context of lattice calculations, Lüscher [11] was the first to propose using finitevolume two-body energy levels for the extraction of scattering lengths. In the three-particle sector, perturbative calculations have culminated in Refs. [12,13], where the result is given up-to-and-including order L −6 and L −7 , respectively (see also Ref. [14]). On the other hand, starting from 2012 (see Ref. [15]), the exact three-body quantization conditionnot based on the perturbative expansion-has been formulated in three different settings, which can be referred to as relativistic field theory (RFT) [16,17] non-relativistic effective field theory (NREFT) [18,19] and finite-volume unitarity (FVU) [20] approaches, see also Ref. [21] for a review. 1 In the following years, extensive exploratory studies and theoretical extensions have been undertaken, addressing various aspects of all three approaches . The equivalence of different approaches has been discussed in Refs. [19,21], and demonstrated for RFT and FVU in Ref. [36]. It should be especially stressed here that the perturbative expansion of the quantization condition in powers of 1/L reproduces the results of Rayleigh-Schrödinger perturbation theory for the ground state [26]. The expansion of the first excited level can be also obtained in this manner [39]. Finally, alternative approaches to study many-body dynamics are provided by the HAL QCD method [48,49], the optical potential method [50], a method of the spectral density functions [51] and the variational method [52][53][54]. The volume dependence of bound states has also been studied in Ref. [55]. Moreover, multi-particle interactions in lower-dimensional systems have been addressed in Refs. [56,57].
It should be stressed that for the extraction of the many-body interaction parameters (effective couplings) from the measured energy levels on the lattice, the perturbative expansion in powers of 1/L provides a very convenient tool, whereas for the study of resonances, an exact quantization condition should be used.
In our opinion, despite the recent progress in the field, there still remain issues requiring further scrutiny. Namely: i) Using the exact quantization condition for the derivation of the perturbative expression for the energy shift is a rather complicated enterprise. On the contrary, the Rayleigh-Schrödinger perturbative expansion provides a straightforward path to a final result. However, the inclusion of relativistic corrections, up to now, has been only performed in the ground state. Since at order L −6 several effects (relativistic corrections, two-body effective-range corrections, three-body force) contribute to the ground-state energy, it would be important to extend this to excited states in order to separate these effects from each other.
ii) The NREFT three-particle formalism provides a very simple derivation of the quantization condition. Thus, it would be interesting to systematically study the size of the relativistic effects in finite-volume energy levels.
iii) All the above derivations implicitly assume that the corrections, exponentially suppressed at large L, are very small and can be safely neglected. In reality, the situation is more nuanced. The three-body coupling, which should be extracted from data, enters first at order L −6 , i.e., it is suppressed by three powers of L with respect to the leading-order contribution. Thus, making L very large will render the extraction more difficult, whereas the exponential corrections are not sufficiently suppressed for smaller values of L. Hence, one has to balance between two opposing trends. It would be interesting to estimate the order of magnitude of exponential corrections, in order to establish, whether a window in L exists for which the extraction can be carried out with a small systematic error.
iv) The energy shift of the N -particle state (both the ground state and the excited ones) at O(L −6 ) depends only on three independent dynamical parameters: the two-body scattering length and the effective range, as well as the three-body non-derivative coupling constant. Consequently, a global fit to the data with different N enables one to put more stringent constraints on these parameters that well result in the substantially improved precision.
v) It would be very interesting to test the extraction method in models where the answer is known from the beginning (e.g., where the perturbation theory is applicable).
In the present paper, we address these problems in complex scalar ϕ 4 theory on the lattice, which has been already used for the similar purpose in [58]. A big advantage of this model as compared to lattice QCD is that it can be simulated with much less computational effort. It allows one to carry out lattice calculations at many different values of L in the sectors with different number of particles. Moreover, the model is perturbative for small values of the coupling constant, so that the values of the parameters, extracted on the lattice, can be confronted with results of direct perturbative calculations. In the present, follow-up work, we, in addition to Ref. [58], calculate the ground-state energies of the 4-and 5-particle states. We also discuss relativistic effects, exponentially-suppressed polarization effects, and compare to perturbative results in finite volume. Finally, as the present work was already in progress, Ref. [59] appeared. There, relativistic corrections were perturbatively included at O(L −6 ) in the case of N identical particles, albeit only for the ground state.
The paper is naturally divided in two parts, containing the derivation of the theoretical framework (Sections 2 and 3), and the analysis of the results of numerical calculations (Section 4). More specifically, in Section 2, we carry out the construction of the effective non-relativistic Lagrangian and the matching of its couplings in great detail. Further, in Section 3, the evaluation of the consecutive terms in Rayleigh-Schrödinger perturbation theory order by order is considered. The energy shift (both the ground state and the first excited state) is calculated up to O(L −6 ), and relativistic corrections are taken into account at the same order. The calculations on the lattice (at present, only ground-state levels in the sectors with different N ), as well as the extraction of parameters through the fit to the measured energy shifts are presented in Section 4. Finally, Section 5 contains our conclusions.

Matching in the NREFT Lagrangian
Below, we shall present the derivation of the energy shift to the N -body state of identical spinless particles with mass m. The derivation proceeds in two steps. In the first step, we perform the matching of the two-and three-body effective Lagrangians, taking into account relativistic corrections. In the second step, we use these effective Lagrangians to perform the calculation of the energy shift in the framework of the Rayleigh-Schrödinger perturbation theory. It will be seen that this procedure allows one to get the final result with a surprising ease and elegance.
In this section, we will focus on the first step of this process. The starting point is the non-relativistic effective Lagrangian containing all two-and three-particle interactions that are relevant in the calculations up to order L −6 . In a general moving frame, it takes the form: It can be checked straightforwardly that four-and more particle interactions do not contribute to the energy shifts at order L −6 . In general, the matching of relativistic and non-relativistic theories proceeds as follows. Let M N denote the relativistic N -particle scattering amplitude. The corresponding nonrelativistic amplitude will be denoted by T N . The matching condition states that with w i = m 2 + p 2 i for incoming particles and, similarly, w i = m 2 + p i 2 for outgoing particles. More precisely, it is understood that both sides of Eq. (2.3) are expanded in a Taylor series in the three-momenta p i , p i around threshold, and the coefficients of this expansion are set equal to each other up to a given order in p. Then, the matching condition expresses the non-relativistic effective couplings g 1 , g 2 , g 3 , . . . in terms of the relativistic amplitudes. Note that both relativistic and non-relativistic amplitudes contain also terms that are not analytic in external momenta p i , p i . However, the structure of such terms is exactly the same in both theories, and they automatically drop from the matching condition without imposing additional constraints on the effective couplings. An explicit example of such a cancellation will be given below. Note also that due to particle number conservation in the non-relativistic theory, the sectors containing different number of particles do not talk to each other. For this reason, the matching in the two-and three-particle sectors can be done without considering other sectors. For a detailed review on the matching of non-relativistic and relativistic theories, see, e.g., Ref. [60].
In the remainder of this section, we will express the effective couplings, g 1 , g 2 , g 3 and η 3 , in terms of observable quantities, such as the scattering length, a 0 , the effective range r 0 and the (subtracted) three-particle scattering amplitude at threshold.

Matching in the two-body sector
Let us start with the matching in the two-body sector. It should be pointed out that it is not enough to perform matching in the center-of-mass (CM) frame of two-particles, because in many-body systems two-particle subsystems have, in general, a non-zero CM momentum.
Let p be a generic three-momentum of a particle, and |p| m. As seen from Eq. (2.1), the kinetic term contains the expansion of the relativistic one-particle energy up to and including O(p 4 ). Further, the terms up to O(p 2 ) are retained in the expansion of the two-body interaction Lagrangian. Note that the term with g 2 is absent if the matching is carried out in the CM frame. When working in an arbitrary moving frame this term must be included. As it will become clear later, the coupling g 2 is not an independent physical quantity-it is uniquely fixed by the requirement of relativistic invariance.
Let us now recall that the expression of the relativistic S-wave amplitude 2 in all frames is given by the Lorentz-invariant expression: Here, s is the usual Mandelstam variable related to the relative three-momentum in the CM frame as k = (s/4 − m 2 ) 1/2 . Further, δ(k) denotes the S-wave phase shift and a 0 , r 0 are the pertinent scattering length and the effective range, respectively. In addition, the scattering K-matrix can be expressed in terms of the amplitude given in Eq. (2.4): As noted above, the matching condition for the scattering amplitudes contains both analytic and non-analytic terms at small momenta. Only the former are relevant for carrying out the matching, and the latter should match automatically. For this reason, it would be technically convenient to "purify" the matching condition from these non-analytic terms (we shall check a posteriori that the non-analytic terms are indeed the same in the relativistic and in the non-relativistic theories and cancel in the matching condition). The goal can be easily achieved in the two-body sector, where the sole source of the non-analytic contribution is rescattering in the s-channel (the so-called bubble diagrams). One could remove this contribution by considering the quantity K 2 and its non-relativistic counterpart V 2 which, in contrast to the pertinent amplitudes, are analytic at threshold. In terms of the K-matrices, the matching condition is written down in a form identical to Eq. (2.3): Namely, in dimensional regularization, which will be used throughout this paper, V 2 is given only by the tree-level diagrams, produced by the non-relativistic effective Lagrangian (2.1). In order to do the matching, one has to perform the non-relativistic expansion of the relativistic variable k 2 : where P = p 1 + p 2 = p 1 + p 2 , Expanding the right-hand side of the matching condition in Eq. (2.6) up to and including terms of order p 2 , we get: Figure 1: Diagrams contributing to the non-relativistic scattering amplitude T (p 1 , p 2 ; p 1 , p 2 ) up to and including O(p 3 ). The shaded squares stand for vertices with the non-derivative coupling g 1 , the filled circles for the vertices with derivative couplings g 2 , g 3 , and the cross corresponds to the relativistic insertion p 4 /(8m).
On the other hand,in the non-relativistic effective theory we obtain: (2.11) From the above two equations, one can readily fix the values for the constants g 1 , g 2 , g 3 : (2.12) Note also that three effective couplings g 1 , g 2 , g 3 are expressed in terms of two dynamical quantities a 0 , r 0 only. This confirms the statement that was made above, namely that one of the couplings is not independent and is fixed through relativistic invariance.

Loops in the two-body sector
Now, we shall explicitly demonstrate that the matching of the couplings, given in Eq. (2.12), ensures the matching of the scattering amplitudes up to and including O(p 3 ) (i.e., the loops that contribute to the amplitude, are automatically the same in the relativistic and non-relativistic theories). This is because additional couplings first emerge at O(p 4 ). All diagrams, contributing to the non-relativistic amplitude at O(p 3 ), are depicted in Fig. 1. Apart from the tree-level contributions, shown in Fig. 1a, at this order one has the one-, two-and three-loop diagrams with the non-derivative vertex g 1 shown in Fig. 1b,c,d, the one-loop vertex with the insertion of the relativistic correction (Fig. 1e) and the one-loop vertex with the insertion of the derivative vertices (Fig. 1f). We shall calculate all of them separately. The tree-level contribution coincides with the quantity V 2 given in Eq. (2.11): Further, using dimensional regularization, we can calculate the one-loop diagram in Fig. 1b: where P µ = (p 1 + p 2 ) µ is the total four-momentum of the pair and D is the number of the space-time dimensions (D → 4 at the end of calculations). Similarly, we have: The insertion of the relativistic correction in the internal line gives: Finally, Adding everything together up to and including O(p 3 ), we get: Now, let us expand the quantity relativistic amplitude up to and including O(p 3 ): Taking into account the relation between the quantities k 2 andq 2 0 , which is given by Eq. (2.7), one can finally verify that the expansion of the expression in Eq. (2.19) exactly coincides with Eq. (2.18) term by term. This is a nice check of the matching procedure, carried out in a moving frame and including leading-order relativistic corrections.

Matching in the three-particle sector
At the order we are working, it suffices to include a single operator that describes the shortrange non-derivative three-particle interactions. Note also that the coupling η 3 should be ultraviolet divergent in order to cancel the divergences emerging from the loops. When the dimensional regularization is used, then η 3 contains the term proportional to (D − 4) −1 and a finite piece. The matching condition that allows one to determine this finite piece, is given by: Again, both left-and right-hand sides of the above equation contain non-polynomial contributions, which should be separated before the matching for the regular part is carried out. In the two-particle case, one has achieved this goal by rewriting the matching condition in terms of the K-matrix, which does not contain non-polynomial contributions at all. Here, the situation is more complicated. It is important to realize that there is a certain ambiguity in the definition of the regular part of the three-body amplitude at threshold, in contrast to the two-body sector, where the parameterization of the amplitude through the effective-range expansion is commonly accepted. Below, we give a simple prescription that operates on the on-shell amplitudes only. The advantage of this prescription is that such amplitudes are physical observables that can be defined both in relativistic and non-relativistic theories. 3 One should also note here that, once the relation between η 3 and the threshold amplitude is established, one can use either of these in the fit of the three-particle energy levels. From this point of view, the matching procedure described below is superfluous. Still, we shall present it here in order to shed light on the physical meaning of the threebody coupling η 3 (a three-body force, in the terminology of the non-relativistic scattering theory).
Let us start by identifying the most singular piece in the relativistic amplitude. It is given by the (infinite) sum of all one-particle reducible diagrams, shown in Fig. 2. In total, there are 9 diagrams of this type, differing by permutations of the initial p l , p m , p n and final p i , p j , p k momenta. The contribution of all such diagrams can be written down in the following form: The two-body scattering amplitude, M 2 , entering the above expression depend on three four momenta as q = p i + p j − p n = p l + p m − p k . These amplitudes are in general offmass-shell, since (p i + p j − p n ) 2 = m 2 and (p l + p m − p k ) 2 = m 2 . However, the singular contribution at q 2 = m 2 contains only on-shell amplitudes. In order to define the on-massshell amplitudes, say, for the first amplitude in the numerator of Eq. (2.21) we perform a Lorentz boost to the center-of-mass (CM) frame of p i , p j and denote the boosted momenta p i * , p j * , p n * and q * = p i * +p j * −p n * . In this frame, the amplitude depends of four variables.
Two of them are related to the total energy in the initial and final state: and two are unit vectors: The partial-wave expansion of the two-body amplitude takes the form: The on-shell two-body amplitude is then defined at s = s : The second amplitude in the numerator of Eq. (2.21) can be treated similarly. It is evaluated at s = (p l + p m ) 2 . Finally, taking into account the identity where only the first term is singular, one could finally define the most singular (pole) piece of the three-body amplitude as: ) . (2.27) The advantage of this definition is that the singular part is defined solely in terms of the observable on-shell two-body amplitudes. Note also that at the order we are working, it suffices to retain only the S-wave contribution in the partial-wave expansion, Eq. (2.24). Consequently, the numerator in Eq. (2.27) simplifies to M 2 (s )M 2 (s), with the on-shell S-wave amplitudes given in Eq. (2.4).
is still singular. Moreover, the singularity depends on the configuration of the initial and final three-momenta. In order to define the threshold amplitude in general, one should first perform the partial-wave expansion and then consider the limit where the magnitudes of external momenta tend to zero. However, to invoke such a cumbersome framework for matching a single coupling η 3 seems to be an overkill, since this goal can be achieved much easier. For example, we could agree to perform the limit p → 0 for a particular configuration of the initial and final momenta. Namely, let e x , e y , e z denote unit vectors in the direction of axes in the momentum space. Then, choose the configuration, for example, as: and p i = −p i . This means that become functions of a single variable λ. In the vicinity of the three-particle threshold, their difference can be expanded in a Lorraine series: where ellipses stand for the terms that vanish as λ → 0, and hence do not contribute to the matching condition at threshold at this order. Here, we have arbitrarily chosen m to be a scale appearing in the logarithm. Any other choice would change the definition of the regular part of the threshold amplitude M Diagrams that contribute to the non-vanishing three-particle scattering amplitude at threshold. The large shaded blobs in the diagram a depict the two-body scattering amplitude at O(p 2 ), the shaded squares correspond to the vertices with the non-derivative coupling g 1 , the cross in the diagram b depicts the relativistic insertion p 4 /(8m), and the large filled circle in the diagram f corresponds to the three-particle coupling η 3 .
Next, we turn to the non-relativistic amplitude and carry out an extraction of the singular pieces in a similar fashion. It is important to note that if the dimensional regularization is used, only a finite number of diagrams, shown in Fig. 3, will contribute to the regular part of the threshold amplitude. In the following, we shall examine all these diagrams one by one.

Singular diagrams in the non-relativistic theory
We start from the diagrams Fig. 3a and Fig. 3b. In fact, one could resum all relativistic insertions in the internal line-this obviously results in replacing the non-relativistic energies in the energy denominator by their relativistic counterparts. Thus, we get: 5 The calculations are performed in the CM frame of three particles, i.e., p 1 + p 2 + p 3 = p 1 + p 2 + p 3 = 0. The quantity T 2 denotes the two-point amplitude in this particular kinematics. For our purposes, it suffices to evaluate this amplitude at order p 2 , since the higher-order terms will lead to the contributions that vanish at threshold (we remind the reader that the denominator in the above equation is of order p 2 ). Using Eq. (2.18) at this order, we get: The amplitude T 2 (p k , q; p l , p m ) is defined in a similar way. The one-loop diagram, shown in Fig. 3c, is given by whereas the two-loop diagrams in Fig. 3d,e take the following form: and Finally, the diagram in Fig. 3f gives:

Evaluation of the singular diagrams
We start with the diagrams in Fig. 3a,b. As already mentioned, there is a subtlety here, which consists in the fact that T 2 , given by Eq.
For the on-shell amplitude, the Mandelstam variable s = (p i + p j ) 2 . Since both particles i, j are on the mass shell, we get also s = 4m 2 + (p i − p j ) 2 , but s = 4m 2 + (p n − q) 2 . Thus, the on-shell amplitude is given by: For the second amplitude, we have the similar expression. Further, in the matching condition, we have the factor w(p i )w(p j ) 1/2 for the outgoing particles and the similar factor for the ingoing particles. In the off-shell kinematics, the second factor will not be equal to w(p n )w(q) 1/2 anymore. Our aim is to establish, which factor should be present here. To this end, let us start from the outgoing state. In the CM frame of the pair (i, j), the above factor is replaced by w(p )w(−p ) 1/2 , where p is the relative momentum of the pair. Performing now the boost to the frame where the total momentum of the pair is P ij = p i + p j , we get where u denotes the unit vector in the direction of P ij . Let now the initial pair have the relative momentum p on shell that means that |p | = |p|, but the direction can be arbitrary. Consequently, under the Lorentz transformation, Taking into account that we are dealing with the identical particles, i.e., w(p)w(−p), we finally get: The brief summary of this a bit lengthy discussion is simple: at O(p 2 ), the matching condition for the on-shell amplitudes looks as follows Now, the evaluation of the finite part in T is straightforward: , (2.45) and The second term in Eq. (2.44) can be simplified by expanding both the numerator and denominator in momenta. The denominator cancels, and we get: As one sees, the difference is a regular function at threshold, albeit both T a and T (a+b) 3 scale as 1/λ 2 at small momenta. The rest of calculations is straightforward, because one can consider non-relativistic two-body amplitudes at O(p 0 ) (i.e., replace these with the scattering length). The diagram in Fig. 3c is ultraviolet-finite in the dimensional regularization. Its explicit expression is not needed. It suffices to note that, in a given configuration of external momenta, so that the constant piece is absent. Further, the diagrams in Fig. 3d,e are ultravioletdivergent: where µ denotes the scale of the dimensional regularization, and factor 9 comes from summing all permutations. The finite parts of the above diagrams δ (d) , δ (e) are pure numbers, given by the integrals over Feynman parameters, see Appendix A.
Using now the matching condition for the coupling g 1 , it can be seen that the threebody coupling should contain the following divergent part that cancels the divergences from the loops: with η r 3 (µ) being the renormalized coupling. Choosing, for simplicity, the scale µ = m in this low energy constant, we may rewrite the matching condition as: whereas at an arbitrary scale Here, the threshold amplitudeT is defined as: Finally, using Eq. (2.30), one may expressT through the relativistic threshold amplitude: with δ (d) , δ (e) as in Appendix A. In principle, Eqs. (2.52) and (2.55) solve the problem completely: they express the coupling η r 3 (m) in terms of physical observables in the twoand three-particle sectors. A non-trivial statement, which is implicit in these equations is that threshold singularities of the three-particle amplitude are produced by two-body rescattering processes only-that is, by subtracting diagrams that describe exactly these processes, one has to arrive at a non-singular expression. Mathematically, this is equivalent to the statement that the limit λ → 0 in Eq. (2.54) exists.

Energy shifts in the ground and excited states
In this section, we aim to derive energy shifts of the ground and excited states of multiparticle systems. Following Ref. [12], this calculation can be easily done within the framework of Schrödinger perturbation theory.
Let us first set up the notations. The energy levels of N identical particles in a free, non-relativistic theory in finite volume are given by: where p i is the momentum of the i-th particle. The set n = {n i } labels the unperturbed levels. For example, the ground state corresponds to n i = 0, i = 1, . . . N .
We shall further see that in a regime, in which |a 0 /L| 1, and mL 1, the effect of the interactions on the finite-volume spectrum in the relativistic theory can be treated as a perturbation to the free non-relativistic energy levels. Using the canonical procedure, one can construct the Hamiltonian of the system from the Lagrangian density, given in Eq. (2.1): where the space integral is taken over a finite cube with side L, and In the above expressions, H 0 denotes the free Hamiltonian, whereas In a finite volume, free fields can be expanded in a Fourier series over the creation and annihilation operators a † (p), a(p): Here, the momentum p runs over discrete values p = 2πn/L , n ∈ Z 3 , and the creation/annihilation operators obey the following commutation relations: The properly normalized unperturbed eigenstates are given by: These eigenstates obey the eigenvalue equation: One may also define the matrix elements of the interaction Hamiltonian between the unperturbed states: The Rayleigh-Schrödinger perturbation theory enables one to calculate the energy shift of a state n using the interaction Hamiltonian H I . A standard formula, which can be found in the quantum mechanics textbooks (see, e.g., [61]), reads: Note that we decided not to display here a (rather cumbersome) fourth-order formula, which is also needed in the calculations at order L −6 .

Two-particle levels
As a warm-up example, we consider the well-known case of two particles in the center-ofmass (CM) frame. The 1/L expansion of the energy shift for the two-particle ground state has been known for a long time up to and including O(L −5 ) [9,11]. At this order, the effects of relativity are irrelevant and one can work in a non-relativistic setup. The next order was worked out in the NREFT approach [12], and in a fully relativistic setup [26,27]. As pointed out in Ref. [26], these two calculations differ by a term that vanishes in the limit m → ∞, and can be referred to as "relativistic corrections". Recently, in Ref. [59], relativistic corrections to the ground state were derived in the NREFT with a method, which is identical to the one used here. Below, for illustrative purpose, we shall rederive 1/L expansion of the energy shift up to and including O(L −6 ) by using Rayleigh-Schrödinger perturbation theory.
The matrix element to the interaction Hamiltonian between the two-particle unperturbed states is given by: (3.10) Calculating the matrix element, one gets: where the sum in the last term runs over the two-element permutation group S 2 for the indices i, j. The relation between g 1,3 and a 0 , r 0 is given in Eq. (2.12). The term with g 2 does not contribute in the CM frame and is therefore omitted here. Now, using Eq. (3.9), we can calculate the energy shift of the ground state order by order in perturbation theory. It is easy to verify that, at lowest order, only terms with the non-derivative coupling g 1 from the interaction Hamiltonian, H I , contribute: In order to evaluate the next-order result, we use Eq. (3.11) and get: Substituting this expression into the second-order energy shift, Eq. (3.9), and using E 0 = 0, (3.14) Here, one has used p = 2πn/L in order to display the powers of L explicitly. It is already clear that one has a double perturbative expansion in powers of 1/L. First, each order in Rayleigh-Schrödinger perturbation theory brings a factor L −3 associated to the momentum sum and a factor L 2 , emerging from the propagator, L −1 in total. The pertinent expansion parameter here is a 0 /L. Further, each derivative vertex is suppressed by powers of L, since the momentum p scales as 1/L. In this case, the pertinent expansion parameter is 1/(mL). Consequently, in order to obtain an expression of the energy shift to a given order in 1/L, it suffices to retain finite number of terms in the effective Lagrangian and truncate the perturbative series at a fixed order. Furthermore, the momentum sums emerging in Eq. (3.14), are formally ultraviolet divergent. In order to be consistent with the matching condition, this divergence should be subtracted in the same way as in the infinite volume, i.e., by using the MS scheme in dimensional regularization. In this particular case, it is very simple. Namely, one might rewrite the first sum as: Here, we have added and subtracted the term with n = 0 in the last sum. Now, the first term is perfectly convergent, and the last term can be calculated by means of the Poisson summation formula, applying the MS scheme to the |n| = 0 term: This way, one gets: π |n| e −2π|n| −8.9136329. (3.17) Acting in the same way, one can easily see that 18) and finally, using the matching condition for g 1,3 , one arrives at the last line of Eq. (3.14). Continuing in the same manner, one can calculate corrections up to and including order L −6 . The result is given by the following expression: This, of course, agrees with the result of Refs. [26,27] as well as Ref. [59]. Moreover, relativistic corrections are easy to identify. They emerge from the contributions from H 3 and H r . If they are dropped, the result in Ref. [12] is recovered. Next, we turn to the relativistic energy shift for the first excited state in the twoparticle sector. We start by defining the unperturbed states. In the non-interacting theory and for identical particles, this state is threefold degenerate in the relativistic, as well in the non-relativistic case and has unit back-to-back momentum in the units of 2π/L, directed along the x, y or z-axis. These three states can be classified in irreducible representations of the cubic group. There is a one-dimensional irrep A + 1 , and a two-dimensional irrep E + . The two-particle A + 1 state is given by: where p x = 2π L (1, 0, 0), and so on. The basis vectors of the irrep E + are given by: Because of the cubic symmetry, any of these two states can be used for the calculation of the energy shift.
The calculation of the energy shifts proceeds analogously to the case of ground state. We do not display the intermediate steps here. Note that the energy shifts are calculated with respect to the relativistic unperturbed energy of the first unperturbed state E 1 = 2 m 2 + (2π/L) 2 − 2m. It is easy to see that the leading contribution to the energy shift in the irrep E + arises from the d-wave interaction and is beyond the accuracy we are calculating: The shift in the irrep A + 1 at order L −5 was calculated in Ref. [11]. In the present paper this result is extended up to and including order L −6 , where the method, used in Ref. [11], does not directly apply for the calculation of the relativistic corrections. The final answer reads: Here, π |n| e −2π|n| −1.2113357 , (3.25)

Ground state of N particles
If we are dealing with more than two particles, two new effects arise. First, the two-particle subsystems can be in the moving frames, even if the total momentum of the N -particle system vanishes. Therefore, H 2 in Eq. (3.2) starts to contribute (still, as we will see, this contribution is only non-vanishing for excited states ). Second, one has to take into account the contribution from the non-derivative three-particle interaction, which is parameterized by the coupling constant η 3 . Four-and more particle interactions do not contribute at the order we are working. The expression in Eq. (3.11) can be easily generalized to more than two particles. For instance, for three particles, we have: where i = j = k. This expression is easy to understand. In the first line, there are three pairwise interactions, in which k labels the spectator particle; the second line stems from the three-particle contact interaction; and the third line corresponds to the kinetic term at O(p 4 ). Note that, in this case, the indices in the last sum run over the three-element permutation group S 3 . For more than three particles, the generalizations are straightforward at this order and include only additional combinatorial factors. In particular, there exist N 2 pairwise interactions, N 3 three-particle interactions, and N ! permutations in the kinetic term.
Further, we shall take an advantage of the fact that the expression without relativistic corrections has been already derived previously exactly with the same method [12]. The remainder can be calculated relatively easy. The final result for the N -particle ground state, expressed in terms of the threshold amplitudeT , has the form: (3.27) Here, the quantities I, J , K were defined above, whereas the quantities Q −102.1556055 and R 19.186903 are the quantities defined in Ref. [12] 6 . The factor (Γ (1) + ln 4π) emerges, because the quantities Q, R in Ref. [12] are defined in the MS renormalization scheme rather than MS, which is used in the present paper. Further, we confirm that the above equation agrees with Eq. (42) of Ref. [59], if expressed in terms of the renormalized coupling η r 3 . In order to compare with the result of Refs. [26,27], one has to carry out the matching of the threshold amplitude, defined in those papers, to the quantityT . We do not do this (rather straightforward but tedious) exercise here.

Energy shift for the three-particle excited states.
In this section we discuss the calculation of relativistic corrections to the excited state levels in the three-particle state. To the best of our knowledge, this issue has not been addressed so far in the literature and appears here for the first time. Note also that the combinatorics in the excited states becomes rather cumbersome. For this reason, we first focus on three particles.
Let us start with the definition of the state vectors. In the three-particle system in the CM frame, the first excited unperturbed state contains one particle at rest, whereas two other particles have back-to-back momenta of unit magnitude, in units of 2π/L. Again, the first excited states contains only the A + 1 and E + irreps of the cubic group. The pertinent basis functions are: where α = 1, 2 and the basis vectors |(p i , p j ) A + 1 and |(p i , p j ) E + , α with back-to-back three-momenta p i and p j are given by Eqs. (3.21) and (3.22), respectively.
As a simple check, let us calculate the energy level shift at leading order: , It is seen that the leading-order result of Ref. [39] is readily reproduced. One might continue along the same path and evaluate all contributions up to and including O(L −6 ). However, since the result in the non-relativistic limit is already known [39], one can achieve the same goal with a significantly less effort, evaluating only the relativistic corrections to this result. Below we shall demonstrate, how this can be done. Let us start from the energy shift in the irrep E + , which does not contain the threebody coupling η r 3 . The non-relativistic result, given in Ref. [39], reads: Here, where the numerical constants y 1 = 2.984094, y 2 = 3.001706, y 4 = −28.89478538 are complicated combinations of the infinite sums. The quantity I 2 is defined as: where n 0 = (0, 0, 1) is a unit vector [39]. Relativistic corrections to this result emerge from three different sources: 1. shifting the effective range r 0 → r 0 − 1/(a 0 m 2 ), according to Eq. (2.12): (3.33) 2. contributions, containing the coupling g 2 , which vanishes in the non-relativistic limit. This contribution describes the CM motion of the two-particle subsystem in the threeparticle system at rest: 3. contributions, emerging from the O(p 4 ) kinetic term H r : Bringing everything together, we get: where and ∆E E + 1 (non-rel) as in Eq. (3.30). Next, we turn to the energy shift in the irrep A + 1 , which, in contrast to the E + , contains contribution from the three-body coupling. We shall again be using a shortcut, calculating only the relativistic corrections to the result given in Ref. [39]. To this end, comparing first the expressions for the ground-state energy shift in the non-relativistic limit, we read off the relation between the three-body coupling η r 3 (m) and the threshold amplitudeM from Ref. [39]: where A collects all numerical factors that account for the different regularization and conventions, used in Ref. [39] and the present article: 1278.49567 . (3.40) Here, δ 887.65392 andQ −105.0597779,R −32.60560475 are the counterparts of Q, R in the cutoff regularization. All these quantities have been defined in Ref. [39].
Next, let us consider the expression for the energy shift of the first excited state in the A + 1 irrep in the non-relativistic theory, given in Ref. [39]: where where y 1 0.279070, y 2 8.494802 and y 4 −172.001650 are numerical constants, defined in Ref. [39]. Using now the relation (3.39), one may rewrite last of the above equations as: Now, we can readily evaluate the relativistic corrections to this result from three different sources, listed above: 1. shifting the effective range r 0 → r 0 − 1/(a 0 m 2 ): 2. the CM motion of the two-particle subsystems: (3.45) 3. the kinetic term: Finally, expressing the three-body coupling in terms of the threshold amplitudeT , we arrive at the following relation: The equation (3.48) represents our final result for the energy shift in the first excited state, irrep A + 1 . The last two terms in the brackets are the relativistic corrections.

Leading contribution to N -particle excited states
To conclude this section, we focus on the leading contribution to the first N -particle excited states. For this, we need the N -particle generalization of the A + 1 and E + states, which must contain all possible pairs with back-to-back momentum in one specific irrep: Now, using the N -particle generalization of Eq. (3.26), we calculate the leading effect to the N -particle energies with n 2 = 1: The extension to higher orders in the previous equation is straightforward but rather tedious.

Application: N -particle energy levels in ϕ theory
In this section, we investigate the extraction of the scattering parameters on the lattice in the complex ϕ 4 theory. For the extraction, Eq. (3.27) will be used. The main goal is to explore strategies and challenges in the determination of the three-particle coupling from lattice simulations, using the fit to the ground-state energy spectrum in the sectors for two, three, four and five particles. Note that, to this end, in Ref. [58] only the two-and threeparticle energies were used. In addition, we study the impact of including relativistic and exponentially-suppressed corrections. Moreover, a comparison between the obtained results and perturbative predictions is made. In this manner, we confront the lattice description of the ϕ 4 theory to the continuum perturbation theory.

Perturbative results for the ϕ 4 theory in the continuum
Before turning to the lattice simulations, we discuss here perturbative results for the ϕ 4 theory in the continuum. This will provide useful check for the scattering parameters, obtained from lattice simulations in this work. The starting point is the continuum Lagrangian of the complex ϕ 4 theory: Here m 0 denotes the bare particle mass and λ 0 the dimensionless coupling constant. The first quantity of interest is the two-particle scattering amplitude. Up to one loop, the contributing diagrams are shown in Fig. 4. The amplitude is then given by with the choice of momenta as in Fig. 4, and the loop integral being We shall use the dimensional regularization and MS scheme for renormalization: In the two-particle sector, the two quantities of interest are the scattering length and the effective range. In order to compute them, we just need to compute the expansion up to O(p 2 ) in Eq. (4.3): (4.7) Carrying out the renormalization, λ 0 is replaced by the renormalized coupling λ r . The two-body scattering amplitude is then given by: Now, by comparing to the expansion of Eq. (2.4): we can obtain the one-loop values of the scattering length and the effective range: (4.10) Note here that in the running coupling λ r (m) the scale µ is set equal to m. We now turn to the computation of the subtracted three-particle scattering amplitude at threshold,T . At leading order, the only diagram that contributes is like diagram (a) in

Exponentially-suppressed corrections
The Eq. (3.27), which is used for the extraction of the scattering parameters from lattice simulations, contains only terms that vanish like powers of L. There are, however, exponentially-suppressed effects, which might be relevant. In order to take them into account, one can apply perturbation theory in the continuum and evaluate the corrections in the parameters that enter Eq. (3.27). In this manner, it is possible to explicitly include the leading exponential corrections into the fit.
We start from the single-particle mass, where the leading-order result of Ref. [62] will be used: (4.12) In this formula, M ϕ (L) denotes the finite-volume mass and m the infinite-volume one. Further, c M , c M are constants whose exact value is not of interest at this stage, and K ν (z) are the modified Bessel functions of second kind. The leading exponentially-suppressed finite-volume correction to the scattering length can be calculated, following the strategy explained in Ref. [58]. The result is given by: (4.13)

Lattice simulation
For the numerical approach to complex scalar ϕ 4 theory, the discretization of space-time is required. This can be realized by introducing a finite hypercubic lattice of a volume T · L 3 .
Here T denotes the Euclidean time length of the lattice and L is the spatial extension. Such a four-dimensional lattice can be described by means of discrete points x µ , separated by the lattice spacing a. In addition, periodic boundary conditions are assumed in all four directions. These are implemented by requiring where ϕ (x µ ) ∈ C is the complex scalar field and e µ denotes a unit vector in the direction of axis µ.
For the numerical simulation of field configurations, a discretization of the action is mandatory. This is realized by replacing the four-dimensional integral by a sum, as well as the partial derivatives by finite differences. Furthermore, introducing a hopping parameter κ, the fields are rescaled, according to Similarly, the dimensionful mass m is removed from the action by converting the terms, depending on absolute values of ϕ x to a binary quadratic form. For that purpose, we use the following relations: Taking this into account, the discretized expression for the action becomes: In the simulations, particle mass and coupling constant are set to m 2 = −4.9 and λ 0 = 10.0, respectively. This leads to λ = 0.253308, and κ = 0.159156 in the discretized action, Eq. (4.17). This choice-identical to Ref. [58]-has the advantage that it leads to a small renormalized coupling constant, which strongly suppresses excited states. By means of Metropolis-Hastings simulations, ensembles are generated with constant T = 48 and varying L in the interval L ∈ [6,24] (from now on, L and T will be given in lattice units). More extended lattices are not studied, because of the increased numerical cost. Likewise, smaller lattices with L < 6 are not taken into account, as they are dominated by the boundaries, and marginally incorporate physical information. The different volumes and the corresponding number of produced configurations, as well as the volume-dependent number of Metropolis updates per thread are summarized in Tab Table 1: Summary of simulations used in this work. T = 48 is kept throughout. n conf labels the number of field configurations. The number of updates per thread, ud/thr, is provided as well. In order to decrease correlation, only every 1000 th configuration is saved.

Effective potential
The ϕ 4 theory has two distinct phases: the symmetric and the broken one. Before performing simulations, one wants to know, which phase is realized by a given choice of parameters. Our framework, used for the analysis of the multi-particle energy levels, can be used in the symmetric phase only, where the particle, described by the field ϕ, has a non-zero mass.
In order to answer the above question, note that measuring the effective potential allows one to distinguish between the two phases: a monotonically growing potential is associated with the symmetric phase, whereas a "Mexican-hat" potential characterizes the broken phase.
The effective potential of the ϕ 4 can be computed by including a term in the actionwith the coefficient J-that explicitly breaks the U (1) global symmetry. The continuum Lagrangian in the presence of the symmetry-breaking term has the form: As can be seen, the shape of dV d|ϕ| is that of a positive monotonically growing function. This implies a monotonically growing effective potential, V (ϕ). Therefore, we can confirm that the set of parameters, chosen for this work, corresponds to the symmetric phase of the ϕ 4 theory.

The operators and correlation functions
With the configurations from Tab. 1, the one-to five-particle correlators are computed. The errors are determined from the corresponding bootstrap samples (standard deviation).
Then, particle energies are extracted by fitting the analytic form of the respective correlation function to the data. The N -particle correlation function can be written as: (4.24) We restrict ourselves to N = 1, . . . , 5 henceforth. In Euclidean time, the temporal evolution of an operator generally readsÔ Considering this definition, with the source operator at the initial time t in = 0, and the sink operator at the final time t out = t, Eq. (4.23) can be expressed as The trace occurring here is due to the periodic boundary conditions. It is the origin of the undesirable thermal pollutions of the correlation functions, that is, effects from particles propagating across the boundary backwards in time. In order to find the analytic form of the correlation function, Eq. (4.23) has to be expanded in an infinite sum of eigenstates, resulting in: Keeping only the leading terms, the explicit expressions of the one-to five-particle correlation function follows are: The first term in each of the previous equations includes the energy to be determined. Except for the single-particle correlator C 1 all functions exhibit thermal-pollution terms in addition to the "cosh" signal. In case of C 2 , this term is time-independent, whereas the remaining ones contain time-dependent thermal effects. Although these vanish in the large T -limit, they will taken into account, as the analysis is restricted to the finite lattice volumes. For instance, time-independent terms can be removed by shifting the affected functions: This also allows to reduce the correlation between different time slices. The analytical form of the shifted correlators will be as in Eqs. (4.28)-(4.32), with cosh replaced by sinh.

Fit strategies and parameter extraction
For the parameter extraction, the correlated fits are performed. When the correlation functions C N (t) with N > 2 are examined, energy levels corresponding to lower numbers of particles, contribute to the time-dependent thermal pollution terms, cf. Eq. (4.28) to (4.32). Those terms will be taken into account in the fit. In order to increase the fit sensitivity to the N -particle energy E N , lower lying energy levels, already determined in previous fits, may be used as priors for constraining the fit. This procedure is shown in the following scheme (4.34): Hence, the best fit parameters obtained previously-shown on top of the arrows-enter the current fit as priors. The χ 2 -function has to be modified accordingly: (4.35) Here, P p i denotes the prior applied to constrain the i th parameter p i and ∆p i the corresponding standard error. Finally, the energy shifts with respect to the free energies are computed as: including the corresponding bootstrap samples. Note that the volume-dependent mass is used for the computation of the shift, as it has been shown to cancel the leading exponentially-suppressed effects [58].
The first quantity to be studied is the single-particle mass M ϕ (L) and its volumedependence. This is relevant because the mass enters in the threshold expansion in Eq. (3.27). We will use the infinite-volume mass, m, obtained by fitting Eq. (4.12) to the M ϕ (L) data. In practice, using m or M ϕ (L) in Eq. (3.27) induces only higher-order exponentially-suppressed effects.
Once the infinite-volume mass is determined, different strategies are carried out to fit the energy shifts: 1. Fixed-N fits. Here the energy shifts of N particles are fitted separately. For N > 2, we also employ a 0 and r 0 as priors. For this, the a 0 and r 0 , extracted from the twobody sector, are used as input for the analysis of the three-to five-body sector. Priors in r 0 are crucial to disentangle r 0 andT , as the two parameters enter at the same order. This also increases the fit sensitivity for the threshold amplitudeT . According to the prescription in Eq.(4.35), the modification of the χ 2 -function is given by (4.37) 2. Unconstrained global fits. This fit includes all data with N = 2−5. It is based on the determination of the values of {a 0 , r 0 ,T } that minimize the sum of all contributing χ 2 -functions, i.e., Here,p denotes a vector containing all fit parameters. Furthermore, d labels the data, L i denote the sizes of different boxes used in the calculations, and D represents a complete set of ensembles.
3. Prior-constrained global fits. In this case, the data with N = 3 − 5 is considered, and we include priors in a 0 , r 0 from the ∆E 2 -fits.
Finally, it is very important to check the convergence of the 1/L expansion of the energy shifts. For this purpose, we study the phase shift, extracted from the Lüscher formalism [63], and compare it to the scattering parameters, obtained from the various fits to Eq. (3.27). According to [63], the S-wave phase shift is determined by and the Lüscher zeta function Z 00 . The momentum k can be calculated as Furthermore, the effective-range expansion at order k 6 reads where P denotes the S-wave shape parameter. Starting from this expansion, four different models are set up, as summarized in Tab. 2. They are fitted to the phase-shift data, calculated from Eq. (4.39) using bootstrap samples. This allows to extract both a 0 and r 0 independently from the N -particle ground-state fits. In case of model (IIa) and (IIb), the S-wave shape parameter P is obtained in addition. The "a" and "b" models should be equivalent, and they allow us to explore the systematics of different fitting techniques. Table 2: S-wave phase-shift fit models.

Results
In this section, the numerical results are presented. The material is organized as follows. First, the volume-dependence of the single-particle mass is analyzed. Then, the energy shifts are investigated at order L −5 . After that, the results, obtained from different types of fits at O(L −6 )-described in Section 4.4-are presented. The effects of relativistic and exponentially-suppressed corrections are studied separately. Finally, the results, obtained from the lattice simulations, are compared to the corresponding perturbative predictions.

Volume dependence of the mass
The results from the fits of the mass to Eq. The corresponding curve is shown in Fig. 6, suggesting that a rough approximation M ϕ (L) ≈ m is valid for L ≥ 20, or equivalently, M ϕ (L)L ≥ 4. In the forthcoming analysis, the numerical result in Eq. (4.42) will be used, when fitting to Eq. (3.27 Table 3: Results, obtained from the single-particle mass fits. Single-particle mass fit Figure 6: Volume dependence of the single-particle mass. The lattices with temporal extension T = 48 are considered. The blue curve corresponds to the fit in the range [8,24] for L.

N -particle ground state to order L −5
We start by studying the energy shifts at order L −5 . At this order, all energies are solely described by the scattering length. First, the ratios of energy shifts are explored. Because of the binomial coefficients, occurring in the N -particle threshold expansion, Eq. (3.27), these ratios behave as: The results are shown in Fig. 7, where it can be seen that the values are in a right ballpark. The accuracy of the energy shift ratios, however, does not suffice for reliable parameter extraction. For this reason, Fig. 7 should only be considered as a rough quality test of the simulation. Next, we perform fits to the threshold expansion up to the same order. To this end, we truncate Eq. (3.27) at O(L −5 ), and include only the data at larger values of L, for which this approximation should be valid. We separate fits for N = 2, 3, 4 and 5, and obtain a 0 . The best fit range is determined by varying the lower bound of the fit interval, while the upper one is kept fixed at L = 24. The best choice turns out to be [14,24], for which the results are summarized in Tab. 4. We also perform a fit, including the exponentiallysuppressed corrections (EC) of Eq. (4.13). As expected, the results with and without EC match reliably. This means that, for the values of L, considered here, the EC are weak, and vanish within the range of statistical errors. Moreover, the fits for different values of N agree nicely within the uncertainty.  Table 4: Comparison of the results obtained from the energy shift fits to order L −5 in the interval [14,24]. We provide the results with as well as without EC

Fixed-N fits to order L −6
In order to extract information on the three-body interactions-parameterized by the quan-tityT -the fits at order O(L −6 ) should be performed. We start with the fixed-N fits as described in Section 4.4. We first fit ∆E 2 , and use a 0 , r 0 as priors for N > 2. The selected results come from the interval [9,24]. We shall try various fit models to explore the effects of EC and relativistic corrections (RC  The results, obtained from the fixed-N energy shift fits at order L −6 in the interval [9,24] with and without EC and RC. The results for a 0 and r 0 obtained from the ∆E 2 -fit serve as priors in the remaining fits.
cf. Tab. 4. As can be seen, the fit quality-represented by χ 2 /ndof and p-value-is not affected by the RC. Yet, it slightly improves when EC is included. On the other hand, both RC and EC modify significantly the best fit parameters. Since the change is comparable to-or even larger than-the statistical error, It is important to include these effects in the fit.

Global N -particle ground-state fits
In order to increase the fit sensitivity forT , it is appropriate to try another strategy, namely, a combined fit of Eq. (3.27) to the energy shifts ∆E N =2,3,4,5 . Furthermore, global fits to the ∆E N =3,4,5 data with prior-constrained scattering length (or both a 0 and r 0 ), determined from the ∆E 2 -fit to order L −6 , are conducted. Again, the lower fit interval limit is varied, while the upper one stays fixed at L = 24. The best results are obtained from the fit interval [10,24]. They are summarized in Tab. 6. For illustrative purposes, a global fit with EC as well as RC without any prior knowledge is shown in Fig. 8.
The results without priors are highly consistent with those from the previous fixed-N fits to order L −6 , cf. Tab. 5. Similarly to the fixed-N fits, EC slightly improve the fit quality, and both RC and EC affect the best fit parameters beyond the size of the statistical error. If only a 0 is used as prior, the fit is consistent but errors remain comparable. Yet, the inclusion of priors in a 0 and r 0 reduces notably the statistical errors, while changing the best fit parameters only by order the statistical error. Remarkably, in this last fit we  Table 6: Comparison of the results, obtained from the global energy shift fits to order L −6 in the interval [10,24] with and without using priors, EC and RC. The four results above the horizontal line correspond to the fits without priors. The last two results are fits, for which the parameter(s) a 0 (a 0 and r 0 ) obtained from the corresponding ∆E 2 -fit within the interval [10,24] serve(s) as prior(s). observe thatT = 0 with 7σ considering only the statistical error.

Phase-shift fits
In order to check the convergence of the 1/L expansion, we will compare the previous fit results to the ones from the full Lüscher formalism-see Eq. (4.39). All fit models, used in the phase-shift analysis, are listed in Tab. 2. In general, the best results are obtained from the fits to the complete set of data points. This is particularly important for the fit models that include the shape parameter, P . The results are summarized in Tab. 7.
As it can be seen, the results from fits, which are using model (Ia) and (Ib), are in good agreement with the previous findings, which are obtained from the perturbative expansion of the the N -particle ground-state energies, cf. with Tabs. 4, 5 and 6. In contrast, the results from the fits with (IIa) and (IIb) do not seem to match those from the energy shift fits. A potential explanation could be provided by exponentially-suppressed effects that have been seen to be relevant before-see e.g. Tab. 6-and are not considered in the phase-shift data points. It is important to note that we must include points with large k 2 to fit the parameter P reliably. However, these points originate from simulations with L < 8, and exponentially-suppressed effects can be very significant. To check this, in Fig. 9 we plot the results from the fits (Ib) and (IIb) along with the predictions, based on the results for a 0 and r 0 , obtained from the global energy shift fits in Tab. 6. If EC are not included in the expansion of the energy levels, the resulting line describes the red data points very well. From this, one can conclude two things: (i) the convergence of the 1/L is not of a concern here, and (ii) EC produce a visible impact. A side observation is that RC are also noticeable in the phase-shift curve.

N -particle ground-state fits with perturbative constraints
As explained in Section 4.1, in the ϕ 4 theory, the different scattering parameters are determined by a single coupling constant. A rough estimate for r 0 andT in perturbation theory, which can be obtained with m ≈ 0.2 and a 0 ≈ 0.4, are: r 0,pt ≈ −50 andT pt ≈ 57000. (4.44) It is however surprising to see that the perturbative results are ∼ 9σ away from the best global fit in Tab. 6. As an attempt to understand the reason for this discrepancy, one could enforce the constraints of Eqs. (4.10) and (4.11) in the threshold expansion [Eq. (3.27)]. This way, the fits would have a reduced number of free parameters. In that context, different ansätze for a modification of Eq. (3.27) are studied. We can either simultaneously constrain both r 0 andT by their one-loop values, or just one of them. The best results, stemming from global fits with L ∈ [12,24] are shown in Tab. 8. There we also show the perturbative expectations for r 0 andT , labeled r 0,pt andT pt . If only r 0 is constrained, the best fit result forT is about 1.5σ away from the perturbative prediction. Alternatively, ifT is constrained, r 0 is also in ∼ 1σ agreement with perturbation theory. Constraining both also produces a similar result. All these fits are seemingly reasonable. Yet, if no constraints are applied, the fit quality-see the p-value-moderately improves, while the  ; ×]-fit. In case of combined fits without a 0 priors, ∆E 2 contributes, in contrast to global fits with priors. All fits are performed in the interval [12,24] for L, and they include EC as well as RC. The perturbative predictions r 0,pt andT pt are calculated from Eqs. (4.10) and (4.11), using the corresponding values of a 0 .
best fit results forT and r 0 strongly differ from the perturbative predictions. This remains an open question, that will be addressed in the next section.

Final assessment of results
We have extracted the two-and three-particle scattering parameters using the N -particle ground state of simulations in complex ϕ 4 theory in the symmetric phase. Reviewing the results discussed in the previous sections, one can see that the different fit strategies seem to agree with each other. More specifically, it can be stated that a 0 is consistent, independently of the applied fit method. Nevertheless, we have noticed a concerning tension of ∼ 9σ between our unconstrained fit results and the predictions of the perturbation theory. This deviation affects two parameters, r 0 andT , whose contribution is very subleading-they appear first at O(L −6 ). To illustrate this discrepancy, a comparison between the phase shift and the perturbative prediction is shown in Fig. 10. In this plot, the perturbative prediction uses the value of a 0 from the fit to ∆E 2 at order L −5 (with and without EC), and r 0 from Eq. (4.44), corresponding to the same a 0 . At this stage, it is obvious that the strong deviations cannot emerge from the statistical fluctuations.
In Section 4.5.5, we have also checked that the 1/L expansion seems to be converging properly. In addition, we have explicitly accounted for the leading exponentially-suppressed effects. Besides this, we note that the impact of EC in the perturbative prediction is minimal-see Fig. 10. Even though subleading EC effects could contribute, it seems unlikely that they are sizable.
At this point, there remain very few explanations for this discrepancy. One possibility could be discretization effects. In order to study this, one could perform a Symanzik analysis [64]. The idea is to write down the effective Lagrangian of the theory, including higher-dimensional operators that account for the cutoff effects. The contribution to the Lagrangian in the Euclidean space from the dimension-6 operators is 7 As can be seen, both c 1 and c 2 contribute at tree-level to the two-particle scattering amplitude, and particularly to a 0 and r 0 . Moreover, all three coefficients c 1 , c 2 , c 3 contribute at the three level toT . The crucial point is that the discretization effects are of the order (aM ) 2 ∼ 0.04. However, a measure of the strength of the interactions, M a 0 ∼ 0.08, is of the same order. Thus, it is plausible that cutoff effects could be significant, and so the continuum perturbation theory would not be an appropriate description of the lattice ϕ 4 theory with our set of parameters.
To gain further insight into the above discrepancy and the fit systematics, we study the shape of the χ 2 function for the previously discussed energy shift fits to order L −6 . For this, we use heat map plots to localize emerging minima and identify potential quasiflat directions. Furthermore, the best fit parameters can then be easily compared with perturbative predictions. For a two-dimensional visualization in the r 0 ,T plane, a 0 is fixed to the best fit result. We show the heat map plot for a global fit with EC and RC in Fig. 11. There, one can see that the χ 2 function grows much slower in one particular direction in the r 0 ,T plane. This happens, because in Eq. : χ 2 heat map that corresponds to the global fit at order L −6 without priors, taking both EC and RC into account. The fit interval is [10,24], and a 0 is fixed to 0.439. The black data point with error bars marks the corresponding fit results for r 0 andT , cf. Tab. 6. The black star marks the perturbative prediction.
of measuringT with statistical significance. In the Appendices B and C we discuss further heat map plots. Since the inconsistency between perturbation theory and the unconstrained fit could not be resolved even after many different tests, we conclude that the perturbative results do not describe the lattice ϕ 4 theory properly. Beyond that, Fig. 11 displays the challenges that arise, when extracting low energy scattering parameters on the lattice. Namely, is clearly demonstrates the difficulty of extracting subleading effects and separating two-and three-body effects in the spectrum.

Conclusions and outlook
In the present paper, we discuss the derivation of the formula for the finite-volume energy shift of N identical relativistic bosons in detail, up to and including O(L −6 ). Moreover, in case of two-and three-particle systems, the derivation is extended to the first excited state, which is contained in the irreducible representations A + 1 and E + of the cubic group.
The derivation is performed in the framework of a non-relativistic effective theory, and relativistic effects are included perturbatively. Moreover, if needed, using the same effective Lagrangian, the derivation can be straightforwardly extended to higher excited states, and to excited states of systems containing more particles. This paper contains the necessary prescription in order to carry out such calculations.
We have also calculated finite-volume energy levels in the sectors with one, two, three, four and five identical bosons in the symmetric phase of the complex ϕ 4 theory on the lattice for many different values of the box size L. The wealth of the available lattice data enables us to carry out fits using different fit strategies, aiming at finding an optimal way to extract infinite-volume observables from the measured spectrum.
The main goal of lattice calculations is to establish few parameters (effective couplings), which describe the interactions in multiparticle systems. Only three parameters contribute at order in 1/L we are working. These are: the two-body scattering length a 0 and the effective range r 0 , as well as the non-derivative three-body coupling that can be traded for the threshold amplitudeT . The latter is of our primary concern. It turns out that a 0 can be reliably determined in all fits, whereas r 0 andT are correlated in the fixed-N fits. Only performing a global fit, involving sectors with different N , enables one to fix these parameters at a reasonable precision. The three-particle amplitude,T , turns out to be different from zero at 7σ, taking into account only statistical error. Still, the χ 2 heat map in the parameters r 0 ,T represents a rather elongated ellipse that indicates a strong correlation. It is important to stress also that this property stems merely from the way these two parameters enter the expression of the ground-state energy shift of N -particles at order L −6 . In this sense, it does not depend on the dynamics and will show up in QCD as well.
From this point of view, the study of the excited states becomes particularly interesting. As it is clear from the explicit expressions that are contained in the present paper, in excited states the relativistic and the effective-range corrections contribute already at order L −5 , whereas the contribution fromT emerges again first at order L −6 . Moreover, the contribution fromT is present in certain irreducible representations and is absent in others. All this gives hope that the study of the excited states might help one to break the spell of large correlations between extracted parameters. It would be certainly interesting to study this issue in the future.
The study of the role of relativistic effects, as well as leading-order exponentiallysuppressed polarization effects was one of the main aims of the present work. We have shown that at the values of L considered, both effects are substantial in the sense that including them in the fit leads to the changes in some observables beyond one standard deviation. This circumstance must be kept in mind when analyzing data from lattice QCD.
We have explicitly observed that finite lattice spacing effects might have a substantial impact on the values of the extracted parameters. In our example, the continuum theory badly fails to predict the values of r 0 andT , which are determined on the lattice. This fact should be also remembered when extracting the corresponding observables in lattice QCD that contribute at the subleading order(s) in the 1/L expansion. B χ 2 heat map for the ∆E 2 -fit to order L −6 Fig. 12 depicts the χ 2 heat map, corresponding to the ∆E 2 -fit to order L −6 in the interval [9,24], including EC as well as RC. Here, we follow the visualization convention of Fig. 11. As expected, a distinct minimum that can be observed, is in agreement with the parameters obtained from the corresponding fit (black data point with error bars, cf. Tab. 5). The black star marks the corresponding perturbative prediction. Similarly to other fits in this work, the perturbative prediction is far from the best fit results. : χ 2 heat map that corresponds to the ∆E 2 -fit to order L −6 taking both EC and RC into account. The fit interval is [9,24], and a 0 is fixed to 0.432. The black data point with error bars marks the corresponding fit result for r 0 , cf. Tab. 5. The black star corresponds to the perturbative prediction.
C χ 2 heat map for the ∆E N -fits to L −6 In this appendix, we show some additional heat map plots in order to explore the impact of priors in our fits -see Fig. 13. As already explained, r 0 andT cannot be disentangled in a fixed-N fit with N > 2, because they enter at the same order in the threshold expansion. This is seen very well in Figs. 13a, 13c and 13e, as a flat direction with constant χ 2 . If priors from the ∆E 2 -fit are used for r 0 , the flat direction disappears and r 0 andT can be obtained -see Figs. 13b, 13d, and 13f. f) Heat map for ∆E 5 -fit to order L −6 with priors, cf. Tab. 5; a 0 is fixed to 0.433. Figure 13: χ 2 heat maps a) to f) that correspond to the fits with EC as well as RC in the interval [9,24] for L. The black data points with error bars mark the respective results for r 0 andT from prior constrained fits, cf. Tab. 5. The black stars correspond to the perturbative predictions, calculated from the fixed a 0 values by using Eqs. (4.10) and (4.11).