An alternative scheme for effective range corrections in pionless EFT

We discuss an alternative scheme for including effective range corrections in pionless effective field theory. The standard approach treats range terms as perturbative insertions in the T -matrix. In a finite volume this scheme can lead to singular behavior close to the unperturbed energies. We consider an alternative scheme that resums the effective range but expands the spurious pole of the T -matrix created by this resummation. We test this alternative expansion for several model potentials and observe good convergence.


I. INTRODUCTION
In nuclear physics, there is a hierarchy of Effective Field Theories (EFTs) which all describe nuclear phenomena at a certain resolution scale (for reviews see, e.g., Refs. [1][2][3].) Pionless EFT describes the interactions of individual nucleons at momenta small compared to the pion mass [4][5][6][7][8]. Apart from electroweak interactions, the effective Lagrangian contains only short-range contact interactions between non-relativistic nucleons. It can be understood as an expansion around the unitary limit of infinite scattering length. The breakdown scale of pionless EFT is set by the pion mass, M high ∼ M π , while the typical low-energy scale is M low ∼ 1/a ∼ k. For momenta k ∼ M π , pion exchange can no longer be treated as a short-range interaction and has to be included explicitly. This leads to chiral EFT whose breakdown scale M high is set by the chiral symmetry breaking scale Λ χ [9,10]. The pionless theory exploits the large scattering length but is independent of the mechanism responsible for it. Thus it can be applied to a variety of systems ranging from ultracold atoms to hadrons and nuclei.
At leading order (LO), one needs to resum a momentum-independent contact interaction in order to describe the large scattering length physics. This resummation is conveniently implemented using dibaryon or dimer fields [11]. At next-to-leading order (NLO) the twobody ranges have to be included perturbatively. In the dimer framework this requires one insertion of the dimer kinetic-energy operator between LO amplitudes. At higher orders, the procedure of perturbative range insertions becomes tedious, and a direct calculation of the corrections requires fully off-shell LO amplitudes. To avoid this, range corrections can be resummed by including the effective range in the denominator of the dimer propagator. Early on it was noted that this resummation introduces a spurious pole in the deuteron propagator [12]. Located at a momentum scale of roughly 200 MeV, it is outside the range of validity of the EFT and thus in principle is an irrelevant UV artifact. However, in threeand higher-body systems it can limit the range of cutoffs that can be used in the numerical solution of the scattering equations. In the three-nucleon system, this is especially true in the doublet S-wave of neutron-deuteron scattering (triton channel) unless measures are taken to remove the pole. In the quartet S-wave, due to the Pauli principle, the solution is not sensitive to this deep pole and the cutoff can be made arbitrarily large.
In Ref. [13] it was proposed to partially re-expand the resummed propagators and to use terms up to order n for a calculation at N n LO. Using these "partially resummed" propagators generates all desired terms at a given order, but still retains some higher-order corrections, which have to be small. 1 The first strictly perturbative NLO calculation of nd scattering in the doublet S channel was carried out in [14], implementing the procedure suggested in [15]. Ji et al. [16,17] extended these calculations to N 2 LO and pointed out that an additional three-body term enters at NLO when the scattering length is varied. This is particularly relevant for applications in ultracold atoms and quark mass extrapolations. Finally, Vanasse [18] developed a scheme that avoids the numerically expensive determination of full off-shell amplitudes made in previous perturbative calculations. Overall, he obtains nd phase shifts at N 2 LO which are in good agreement with the empirical behavior up to laboratory energies of 24 MeV. In this paper, we revisit the problem of range corrections in the three-body system from the perspective of the three-body quantization condition in a finite volume, following the formalism developed in Refs. [19,20], see also [21][22][23] for alternative formulations. For simplicity, we focus on the three-boson system, which is known to have the same qualitative features as the neutron-deuteron doublet S-wave channel. We expect the approach of [18] to be problematic numerically in a finite volume. Indeed, in a finite box of size L, the S-wave dimer propagator gets replaced by [19,20]: .
Here, k, k * denote the total three-momentum of a dimer and the magnitude of the relative momentum of two particles, constituting a dimer, in their center-of-mass frame. Furthermore, δ(k * ) denotes the pertinent phase shift and the quantity S(k, k * 2 ) stands for the infinite sum where E is the total energy of the particle-dimer system in the rest frame. 2 In the infinite volume, the sum turns into the integral that can be easily evaluated, leading to a well-known result.
The problem with expanding the finite-volume dimer propagator in a manner proposed in Refs. [14][15][16][17][18] is related to the singularities of the denominator. Namely, from Eqs. (1) and (2) it can be immediately seen that, in a finite volume, the propagator has an infinite tower of poles above the elastic threshold, corresponding to the finite-volume energy spectrum in the two-particle subsystems. In the infinite volume, these poles condense and form an elastic cut. Next, we note that, in a finite volume, the expansion will not work in the vicinity of these poles, producing denominators that are more and more singular. Bearing this fact in mind, we aim at an alternative procedure for removing the spurious poles, which is not based on such an expansion and, hence, high powers of the energy denominator never appear. Below, we shall demonstrate, how this goal can be achieved.
The paper is organized as follows. In In Sect. II we set up the EFT framework which allows one to study the three-particle problem in a systematic manner. In Sect. III we formulate a method that allows one to consistently remove a spurious subthreshold pole from the dimer propagator. In Sect. IV this method is numerically tested within a toy model. The convergence of the approach, as well as the applicability of the power counting is discussed in detail. Finally, Sect. V contains our conclusions.

A. Non-relativistic Lagrangians
We consider the system of three identical non-relativistic bosons with a mass m, described by the field ψ. In this system a non-derivative three-body interaction is required for renormalization already at leading order [24]. The Lagrangian takes the form (only S-wave contributions are shown explicitly): The couplings C 0 , C 2 , describe the interactions in the two-particle sector and can be related to the S-wave scattering length a and effective range r, respectively. D 0 and D 2 correspond to three-body interactions with zero/two derivatives. Higher-order terms with more derivatives are not shown explicitly.
To describe the three-body systems, it is convenient to work in the particle-dimer formalism. The dimer can be introduced as an auxiliary integration variable in the path integral. In this manner, it is obvious that the theory with dimers leads to the same Green functions. The particle-dimer Lagrangian takes the form 3 Here, the ellipses stand for the terms that contain more space derivatives or higher partial waves, d denotes the dimer field, and the sign σ = ±1 determines the sign of the effective range.
In the examples discussed below, we have σ = −1. The two Lagrangians (3) and (4) describe the same physics, so the couplings can be matched to each other. This matching has been considered in the literature many times (see, e.g., Refs. [25,26]) and we do not repeat it here once more. Note only that two couplings C 0 , C 2 (or, the scattering length and the effective range) can be traded for two parameters ∆, f 0 , whereas two other couplings D 0 , D 2 can be expressed through h 0 , h 2 .
In the dimer picture, the three-particle amplitude is expressed through the particledimer amplitude in a closed form. The latter obeys an integral equation (the Faddeev or Skorniakov-Ter-Martirosian equation), which can be readily obtained, considering the diagrammatic expansion of the amplitude. Note that the dimer need not correspond to a physical particle. Within this approach, it is just a useful mathematical tool that makes the bookkeeping of various diagrams extremely simple. In the numerical study that follows, however, we shall adjust the parameters so that the dimer is a stable particle, and use parameter values from the two-nucleon systems. The on-shell particle-dimer scattering amplitude then has a direct physical interpretation.
B. Faddeev equation for the particle-dimer scattering As already mentioned, the particle-dimer scattering amplitude in the non-relativistic effective theory obeys the Faddeev equation where E is the total energy of the particle-dimer system in the center-of-mass (CM) frame, and τ (k; E) denotes the two-body amplitude. It is always assumed that E has an infinitesimal positive imaginary part E → E + iε. As in the Lagrangian (4), we have included only S-wave two-body interactions. Higher-partial wave interactions contribute only beyond the order considered here. The S-wave two-body amplitude in Eq. (5) is given by: where δ(k * ) denotes the S-wave phase shift, and k * is the magnitude of the boosted relative momentum. In the non-relativistic kinematics, Here, m stands for the particle mass. Further, for small momenta, the effective-range expansion can be carried out: where a and r stand for the scattering length and the effective range, respectively. The kernel in the Faddeev equation consists of the one-particle exchange contribution and a tower of the polynomial terms with the increasing powers of momenta, which are obtained from the particle-dimer interaction Lagrangian: where the parameters H 0 , H 2 , . . . can be expressed in terms of the effective couplings in the Lagrangian h 0 , h 2 , . . .. Further, H 0 , H 2 , . . . depend on the cutoff Λ so that the scattering amplitude M (p, q; E) is Λ-independent at a given order in the low-energy expansion.
Carrying out a partial-wave expansion in the Faddeev equation and projecting onto Swave results into: where Z(p, q, E) = 1 2pq ln where the subscript = 0 has been dropped in all amplitudes. Note that this has been done only in order to keep the formulae simple and transparent. If needed, the formalism can be easily extended to include higher partial waves (see, e.g., Ref. [20]). Further, as shown in Ref. [13], introducing a trimer auxiliary field in the Lagrangian along with the dimer field, is it possible to simplify the Faddeev equation. In the kernel of the transformed equation, the three-momenta are traded for the total energy E: ln where γ = √ mE d and E d denotes the binding energy of the dimer. The amplitude, which is a solution of the equation with the transformed kernel, is equal to the original amplitude up to the higher-order terms. It is slightly easier to use the transformed kernel in numerical calculations and we shall stick to this option in the following.
In the presence of a stable dimer, the on-shell amplitude M is related to the particle-dimer scattering phase, according to: This phase is real below the dimer breakup threshold at E = 0.

A. Spurious states
The hard scale M high in the two-body interactions is set by the effective range r. To make further discussion as transparent as possible, let us assume that true dynamics of a system, which at small momenta is described by the non-relativistic effective Lagrangian, is such that no deeply bound two-body states with √ mE 2 |r| −1 emerge. The effective field theory setting in the present form could not be used to consistently describe such states anyway, and we merely discard them (in the two-particle sector, the presence of such states at small momenta will show up only indirectly, through their contributions to the effective couplings). Only shallow bound states with √ mE 2 |r| −1 will be allowed. In particular, in the following, we shall tune our parameters so that only one shallow bound state -a dimer -with the binding energy E d > 0 exists. Hence, the two-body scattering length a must be large and positive, a |r|. After this introduction, let us formulate the problem. If in the Faddeev equation (5), the integration momentum |k| runs from 0 to Λ, the quantity k * varies from k * = √ −mE to k * √ 3 2 Λ (if E < 0, the quantity k * is always real). Thus, the subthreshold amplitude at large momenta enters the equation. In the effective theory, all that can be done is to approximate k * cot δ(k * ) by means of the effective range expansion, which does not make sense at large momenta. One would argue that the behavior at large momenta should not really matter and can be taken care of by an appropriate renormalization prescription. Hence, it would be harmless to extend the integration to high momenta. In reality, however, the situation is more subtle. Let us retain only the first two terms in the effective-range expansion. Then, if r > 0, the two-body amplitude τ (k * ) develops a spurious pole at large momenta: where It is obvious that k 1 and k 2 correspond to the physical dimer and to a spurious deep pole, respectively. Such a spurious pole emerges, because effective range expansion is applied in a region where it is not valid anymore. Including higher orders in the expansion will generate even more spurious poles. An immediate consequence of the emerging spurious pole is that the integration contour hits a singularity where, originally, there was no singularity. It should be understood that the presence of the singularity is not a problem per se: in fact, in a theory where physical deeply bound states are present, there are also singularities and one has to handle them by deforming the integration contour or otherwise. On the other hand, the fact that such a spurious pole contributes to the unitarity relation is a true problem. If in reality there is no such state, there should be no such contribution at all. Even worse, in the case relevant to the two-nucleon problem, a > 0 and r > 0, the spurious pole has a residue with a wrong sign, leading to the negative probabilities. Indeed, as can be seen from Eq. (14), the residues at two poles have opposite signs and, since the dimer corresponds to a true bound state, the second pole has to correspond to a ghost. In addition, it is not immediately clear, how such contributions can be removed by changing the renormalization prescription for the effective couplings, which are presumed to be real.
In the literature, one encounters different prescriptions for treating such spurious singularities. For example, one may keep the cutoff Λ low enough, so that the spurious poles do not appear on the integration path. The shortcomings of this approach, both conceptual and practical, are obvious. First of all, one cannot remove the cutoff and ensure the independence of the results on the regularization. Moreover, the upper bound of the cutoff depends on the order one is working, and on the values of the effective-range expansion parameters. Hence, setting up a universal upper bound is not possible in general.
The power counting of pionless EFT stipulates that the effective range corrections in the three-body system are perturbative, since |a| |r| [15]. This approach is implemented in Refs. [16,18]. It is reminiscent of the threshold expansion of Beneke and Smirnov [27] (see also Refs. [28,29]), and the heavy baryon expansion in Chiral Perturbation Theory [30][31][32]. This approach is based on the observation that the Taylor-expansion of the propagators alters only high-momentum contributions in the Feynman graphs -exactly those, which are responsible for the trouble. Namely, following Refs. [16,18], one may expand the quantity τ (k * ), given by Eq. (14), in series in the effective range r and include the contributions in strict perturbation theory. The energy denominators, (−1/a + k * ) −n , obtained in a result of this expansion, do not produce spurious poles. The resulting Faddeev equation can be readily solved -the solution is written down as a series in powers of the effective range parameter r. The method is very appealing, successful, and fully consistent. However, using this method in a finite volume, following the approach of Refs. [19,20], is not very convenient numerically, since the denominator in a finite volume becomes very singular (cf. the discussion in the introduction). For this reason, in this paper we propose an alternative approach to this problem, where only the spurious pole contribution is expanded. In this manner, high powers of energy denominator never appear. In addition, in our opinion, this method could be even simpler in the applications.

B. Method
Let us in the beginning assume that we work below the dimer breakup threshold, E < 0. The argument is then crystal clear. We start by splitting off two poles in Eq. (14) from each other: Here, the first/second term contain the dimer/spurious poles, respectively. Note now that the second term is, in fact, a low-energy polynomial -since k 2 is a quantity of order of a heavy scale, k 2 ∼ M high , it can be expanded in Taylor series in k * 2 . Doing this, one gets rid of the spurious pole. It should be however demonstrated that the change in the amplitude, which results by replacing the deep pole by its Taylor expansion, can be indeed accounted for by adjusting the effective couplings. Below, we shall demonstrate this by explicit calculations at one loop and interpret this adjustment physically. It is convenient to introduce the notations: as well as and so on. In other words, from the term corresponding to the spurious pole, we subtract its Taylor expansion, up to some order. Further, writing down τ the Faddeev equation can be rewritten in the following form: Note now that in the first equation of the above system, which determines the amplitude M one is looking for, the spurious pole is replaced by its Taylor expansion. Consequently, the culprit has been removed. The question remains, however, whether the effective potential W , which is determined by the second equation, has the same properties as Z, i.e., is given by a sum of the one-particle exchange diagram and a low-energy polynomial. In this case, one could forget about the second equation at all, since the difference between W and Z could be accounted for by a change of the renormalization prescription.
In the following, we expand the quantity W in the Born series in order to study the structure of each term separately. In particular, considering a couple of simple examples at the second order, we verify that W (2) has indeed the structure which was conjectured from the beginning. The generalization to higher orders is clear.
Let us now start with the calculation of W (2) . The quantity Z, displayed in Eq. (9), contains an infinite number of terms, and hence W (2) will contain infinite number of all cross products. To illustrate our statement, we pick out a single term. The simplest choice is the one, proportional to H 2 0 : Note that the sign of iε follows from the prescription E → E + iε. The imaginary part of I 00 is a constant, which depends on the energy E: We assume here that the cutoff Λ is chosen large enough, so the pole is inside the integration region -otherwise, the imaginary part would vanish. Further, the real part is also a lowenergy polynomial: It can be seen that the real part can be removed by altering the renormalization prescription. The sole subtle point is that the counterterms depend on (are low-energy polynomials of) the total three-particle CM energy E which, in the Lagrangian, translates into time derivatives on both the particle and dimer fields. The following discussion demonstrates, how one could circumvent this problem. First, if one is interested only in the on-shell particle-dimer scattering matrix, one could directly use the equations of motion (EOM) in the particledimer Lagrangian, trading the time derivatives for space ones. In the description of the generic three particle processes, however, the dimers may go off-shell. In this case, one should first integrate the dimer field out and then use the EOM for the particle fields that leaves the three-body S-matrix elements unchanged. Applying the same procedure to the imaginary part leads, however, to a conceptual inconsistency, since the counterterms, which are needed to remove it, should be complex. The problem with the spurious poles shows up exactly at this place. Note, for example, that if the cutoff Λ is chosen so small that the integration contour does not hit the pole, then the problem does not arise, since the imaginary part vanishes. It is also clear that one could circumvent the problem, which originates from the use of the effective-range expansion beyond the range of its applicability, by merely dropping the imaginary part by hand (because, in the exact theory, there are no poles and thus no imaginary part).
As a side remark, this discussion also shows how physical deep bound states should be treated. The corresponding poles are physical and cannot be eliminated from the theory. On the other hand, it would be inconsistent to treat them in the present setting explicitly, because their binding energy is determined by the hard scale M high . According to the above discussion, such a deep bound state pole will show up indirectly, through the contribution to the effective couplings, which become complex. In contrast to the case of spurious poles, the imaginary part corresponds to the contribution of the physical deep bound state to the unitarity relation and cannot be discarded. The potential W becomes now a kind of "optical potential" [33], in which the shielded states manifest themselves through the imaginary part. It should be also mentioned that the contribution from the physical states to the imaginary part always comes with a correct sign, in accordance with unitarity.
Next, we shall consider another contribution to the quantity W (2) that will allow us to have a closer look on its structure at small momenta. We shall namely single out the term where both Z are replaced by the one-particle exchange contribution: and ρ 2 = 4 3 (k 2 2 + mE). The integral I pole is ultraviolet-finite, and hence the cutoff Λ can be taken to infinity. Using the Feynman trick, it can be written in the following form: where The integral over the variable y can be performed, yielding: The first term is again a low-energy polynomial (with complex coefficients) and can therefore be discarded, while the second term is not. Expanding the numerator in the integrand in a Taylor series, we get: Next, consider the subtraction integral I subtr = n I (n) subtr /k 2n 2 . The leading term is given by: It is immediately seen that the leading-order term I (1) subtr cancels the leading-order nonpolynomial piece in I pole that emerges from the first term in the expansion in Eq. (28). The higher-order terms such as have the same property. The integral cancels against the next-to-leading order nonpolynomial contribution, emerging from the second term in Eq. (28), and only the polynomial contribution is left at this order. The role of the higher-order subtraction terms is similarthey merely remove the non-polynomial contributions at the pertinent order, leaving only the polynomial parts (as it should indeed be). The general pattern becomes crystal clear already from these examples, and there is no need to consider higher-order terms. To summarize, the quantity W is indeed a lowenergy polynomial up to an order fixed by the order of the subtracted polynomial. The coefficients of this polynomial are energy-dependent and complex. The energy-dependence can be eliminated through the use of the EOM. The imaginary parts, arising from the spurious poles, are artifacts of the use of the effective-range expansion for large momenta. Our prescription consists of dropping these artifacts since, in the full theory, there are no poles leading to the complex potential. Thus, one may finally assume that W = Z, modulo the change in the renormalization prescription.
Final remarks about unitarity are in order. The un-expanded two-body amplitude, which still contains the spurious pole, obeys exact two-body unitarity by construction, whereas this property is lost after expansion. However, the violation is small in the physically relevant region of small momenta, because k * 2 /k 2 2 ∼ M 2 low /M 2 high is a small parameter there. Moreover, the violation of unitarity in this region can be systematically reduced, including higher-order terms in the Taylor expansion. Further, our argument can be extended for the energies above the breakup threshold, E > 0. In this region, it is no longer true that the contributions to the imaginary part of W come solely from the spurious subthreshold pole. In fact, they can emerge also from the denominators, corresponding to the particle exchange between the dimer and spectator particle. This contribution to the imaginary part is physical and should be retained. Note, however, that this contribution emerges exclusively from the region of small integration momenta, where the quantity k * is small. In this region, the quantity f (k * ) is also small (it converges to zero in the Taylor expansion in k * 2 /k 2 2 ). Hence, the corresponding contribution to the imaginary part of W should be small. It can be systematically reduced by including higher-order terms in the Taylor expansion. Thus it can be safely neglected.
It should also be mentioned that the relation of the amplitude to the phase shift is modified along with the unitarity relation, if the subtraction is done. In particular, instead of Eq. (13), one now has: Note that Eq. (31) reduces to Eq. (13) in the limit r → 0, as it should.

C. Order of the subtraction polynomial
It is natural to ask how large the order of the subtracted polynomial in f (k * ) should be. Is it so that, if one subtracts more terms, the accuracy of the method increases? The answer to this question is obviously no. Recall that one has to compensate the subtraction by adjusting effective couplings in the Lagrangian. If one does not have enough couplings H 0 , H 2 , . . ., a further subtraction does not lead to an improved accuracy.
Since the problem is highly non-perturbative, it is difficult to establish the order of the subtraction polynomial a priori without a non-perturbative calculation. We stress that the requirement to promote the three-body interaction to leading order in Ref. [24] was also established by explicitly investigating the cutoff dependence of numerical solutions of Eq. (10). Alternatively, one can analyze the asymptotic behavior of non-perturbative solutions [34,35]. In order to get a first idea on optimal number of subtractions, we start with a perturbative analysis of Eq. (10), being well aware of the shortfalls of this approach.
It is convenient to consider the effective potential W , rather than the amplitude M . It is straightforward to establish counting rules for W in perturbation theory. Indeed, assume that one is using dimensional regularization to tame ultraviolet divergences in this quantity (the use of any other regularization, say, the cutoff regularization, will alter only the polynomial part of W that can be compensated by a choice of the renormalization prescription). Let us now consider the perturbative expansion of the potential W , given by Eq. (20). Each consecutive term in this expansion contains one additional factor Z, f (k * ) and d 3 khence, the power in p increases at least by one, when one goes to higher-order terms. Hence, the most stringent constraint on the number of subtractions arises from the term W (2) . At lowest order, one has to replace f (k * ) by f 1 (k * ). Then, W (2) counts at O(p 3−2+2−2 ) = O(p) according to our power counting. Of course, this counting concerns the non-analytic piece of W (2) only. Furthermore, taking f 2 (k * ) instead of f 1 (k * ), we get the non-analytic piece starting at O(p 3 ), and so on.
Imagine now that we have only one coupling H 0 at our disposal that counts at O(p 0 ). Adjusting this single coupling, one can achieve that Re W (2) = O(p) if f 1 (k * ) is used since the non-analytic piece starts at O(p 3 ). If f 2 (k * ) is used, the non-analytic piece starts only at O(p 3 ) and the leading contribution comes from the analytic piece at O(p 2 ), i.e., Re W (2) = O(p 2 ). Using f 3 (k * ), . . . in the calculations does not lead to the further improvement, since we do not have the H 2 counterterm at our disposal to remove the O(p 2 ) piece. By the same token, using f 3 (k * ) should be optimal in case of two constants H 0 , H 2 . In this case, Re W (2) = O(p 4 ) can be achieved.
Finally, we reiterate that the above discussion should be taken with a grain of salt as it is based on perturbation theory. Hence, the counting rules, given above, can provide only a hint about the optimal number of subtractions in the non-perturbative case. We therefore conclude that it is important to numerically check the expectation, based on the above power counting, in non-perturbative calculations. This goal will be accomplished in the next section.

IV. NUMERICAL TEST
In this section, we shall test the approach described above using explicit nonperturbative calculations. In these calculations, a quantum-mechanical system of three identical bosons, interacting pairwise through some model potential, will play the role of an exact underlying theory. The underlying theory, by definition, does not contain spurious poles. These appear, when one replaces an exact two-body amplitude in the Faddeev equations by the effectiverange expansion. Thus, one may check, whether the results, obtained in our scheme, do indeed converge to the known (exact) result, and estimate the rate of this convergence. We will consider a Yamaguchi potential first and then repeat this analysis for a Gauss potential.

A. Yamaguchi model
As mentioned above, we consider a toy model with three bosons of a mass m, interacting through the Yamaguchi potential [36], as an exact theory. This potential is given by: Here, λ denotes the strength of the potential, and β is related to its range. To connect the parameters of the Yamaguchi potential to the scattering length a and the effective range r, we calculate the two-body scattering amplitude: The on-shell amplitude takes the form: where E p = p 2 /m + iε. Expanding this amplitude and comparing the result to the effectiverange expansion, we obtain: In the three-body sector, the model does not contain a three-body force. This is a valid choice since all integrals are convergent at the upper limit (the parameter β plays a role of the ultraviolet cutoff). The equation for the particle-dimer scattering amplitude M Y (k, p, E) takes the form (see, e.g., the textbook by Schmid and Ziegelmann [37]): where the dimer propagator τ Y (q, E) is given by: with γ = −mλβ 3 /(8π)−β = √ mE d , and the convention E → E+iε is implicit everywhere in Eq. (37).
The one-particle exchange potential Z Y (p, q, E) in the Yamaguchi model can be written down in the following form: The calculation of the amplitude M Y (k, p, E) can be carried out by using standard numerical procedures. We namely use a large momentum cutoff Λ = 1500 MeV to approximate the integral (the presence of the momentum cutoff is not critical since, as said above, the integral converges even in the absence of the cutoff). 5 In the model, the particle-dimer scattering phase shift δ Y (p) is defined, according to: As already mentioned above, below the dimer breakup threshold E < 0, the phase shift δ Y (p) is real, according to the unitarity. Note also that, in order to ease notations, we did not choose the same normalization for the amplitudes M and M Y . This does not cause a problem, since the particle-dimer phase shifts are compared, which are independent on the normalization chosen.

B. Matching of the EFT framework
As stated before, our aim is to compare the solution of the Faddeev equation M Y (k, p, E) with the solution of the Eq. (19), where W (p, q, E) = Z(p, q, E) is assumed. In the calculations, again the hard cutoff is imposed, and two values Λ = 250 MeV and Λ = 600 MeV are used. Note that, in this case, the cutoff plays a crucial role as regulator, since the momentum integrals are otherwise divergent.
Owing to the initial choice of the parameters, both propagators τ Y and τ have a pole at the deuteron energy E d = k 2 1 /m, corresponding to k 1 46 MeV. For the given choice of parameters, the quantity τ exhibits a second, spurious pole at k 2 179 MeV as well, whereas in M Y (k, p, E), such a pole is absent.
In order to apply our method, we define the subtracted propagators τ i (k * ) = τ (k * )−f i (k * ), where i denotes the number of subtractions. Thus, Note also that, for the remaining (shallow) pole in τ i (k * ), the prescription k * → k * − iε is implicit in all above expressions. This corresponds to E → E + iε. Various approximations, which can be constructed within our approach, differ by a) order in the effective range expansion and the number of the three-body couplings H 0 , H 2 , . . . used, and b) the number of the retained terms in the expansion in τ i (k * ). The calculations are done for leading order (LO), next-to-leading order (NLO) and next-to-next-to-leading order (N 2 LO) in pionless EFT. According to the standard power counting in the two-and threebody sectors, the following parameters appear:  Next, we briefly discuss the matching of the low-energy couplings H 0 , H 2 . If there is only one three-body coupling present, as at LO and NLO, it is most convenient to determine it from matching at threshold. For technical reasons, we perform matching of the particledimer scattering phases in two theories p cot δ Y (p) and p cot δ(p) at small, but non-zero value of the momentum p = 0.001 MeV. When the second coupling H 2 is present (N 2 LO), it would be natural to match in addition the first derivative of the function p cot δ(p) at threshold. Equivalently, one could match the value of the function p cot δ(p) at a some value of p above threshold. We have opted for the second option, because it is easier to implement  in our numerical algorithm, and have chosen the value of the second matching momentum p = 10 MeV, which is still quite close to threshold. Below, we shall discuss the matching condition briefly. First, note that the values of the couplings H 0 and H 2 , in addition to the cutoff Λ, depend on the number of the retained terms in the Taylor-expansion of the spurious pole (this latter dependence is not present in LO, because there are no spurious poles at this order). Further, it is seen that the results of the matching for H 0 do not depend on, whether H 2 is included or not. This follows from the fact that the contribution from H 2 is multiplied by a factor (mE + γ 2 ) (see Eq. (12)), which exactly vanishes at the particle-dimer threshold. This is seen in the Table II,  To begin with, we calculate the particle-dimer scattering phase shift δ in the toy model with the Yamaguchi potential, and in the effective theory, amended by our prescription for treating the spurious poles. As mentioned above, we have in fact to deal with two different expansions: the EFT expansion (i.e., including more derivative terms in the Lagrangian that are accompanied with the independent couplings), and the Taylor-expansion of the spurious pole. The convergence of these expansions need to be investigated separately.
Since it turns out to be the most efficient choice, we use the subtracted propagators τ 1 (k * ) and τ 2 (k * ) in the calculations at NLO and N 2 LO respectively. Remember that at LO, no subtraction is needed. Note also that this choice differs from our perturbative estimate in Sec. III by one order. The other possible choices of τ i (k * ) at NLO and N 2 LO, including the one based on perturbation theory, are discussed below. The real part of the results of these calculations are shown in the left part of Fig. 1. It is seen that LO is precise only at small momenta, whereas NLO can describe data at much higher values of p. The situation further improves at N 2 LO, albeit this improvement is very small (practically not visible by a bare eye). In the left part of Fig: 1 the imaginary part of δ is shown. It can be seen that the NLO and N 2 LO results describe the model better than LO, while the N 2 LO results are clearly improved compared to NLO. The errors of the EFT calculation for p > 1/a can be estimated as (p/Λ) n+1 at N n LO. A more detailed evaluation of the EFT errors is presented in the discussion of possible choices for τ i (k * ) below.
Up to now, everything follows the standard EFT pattern. However, in order to answer the question, whether a systematic improvement is achieved in higher orders, as well as to address the subtraction of the spurious pole, a more elaborate study of the problem is necessary. To this end, it is convenient to use the so-called Lepage plots, which will be considered below.

D. Lepage plots and consistency assessment
Lepage [38] has proposed a method, which allows one to check, how well the data are described by an EFT. The method makes use of certain double-logarithmic plots, known as the Lepage plots. Grießhammer [39] has suggested to verify the internal consistency of an EFT along the similar pattern. In the following, we shall adapt these methods for the problem we are working on.
Let us consider an EFT, describing the fundamental theory up to order n. The corrections are of the order [(k typ , p)/Λ b ] n+1 , where k typ ∼ 1/a is a typical momentum in the reaction and Λ b is the breakdown scale of an EFT. For an arbitrary observable, and, in particular, for the three-body phase-shift p cot(δ), we have: This means that Here, c, c and c stand for some constants. The quantity η describes the corrections due to the denominator. It is also assumed that k typ p, this is discussed below. Hence, the slope in a double-logarithmic plot gives the order n of the neglected term. To determine this slope, a linear function can be fitted to the numerical results.
Further, one may check the internal consistency of an EFT without comparing to data at all [39]. Instead, one can compare the results of calculations within the same EFT, at two different values of the ultraviolet cutoff Λ 1 and Λ 2 , Here c(Λ 1 , Λ 2 , k typ , p, Λ b ) is a slowly varying function of k typ and p. Further, the parameter η describes the dependence of p cot(δ EF T (Λ 2 ) ) on p at LO and will be determined from the fit at the LO. The slope in a double logarithmic plot is, approximately, n + 1 − η. Note that the η in the consistency assessment and the Lepage plots does not have to be the same. Since k typ is not uniquely determined and the double expansion in k typ /Λ b and p/Λ b complicates the analysis, it is very useful to stick to the region, Moreover, we choose the cutoff Λ of the order of the breakdown scale Λ b to simplify the analysis. In this region, termed as the "window of opportunity", the dependence on k typ should disappear (recall that, in our case, k typ = 1/a). On the other hand, one cannot use too large values of the variable p, of the order of the hard scale M high of the theory, determined by the effective range and/or ultraviolet cutoff. Hence, ensuring that one can reliably determine the slopes from the fits in the "window of opportunity" is a non-trivial exercise. For example in Fig. 2 we see that around 80 MeV a spike appears. This is due to Re[δ] = 0 (compare to Fig. 1) in the denominator. This spike will change the slope in this region, so the "window of opportunity" is restricted to be below this. With this we choose the window between 42 MeV and 55 MeV for the δ-slopes. MeV for all orders (gray shaded region). The spike (zero of δ (Fig. 1)) around 80 MeV limits us to low energy regions. Note that the LO result does not predict this zero, therefore the spike is not visible in the consistency assessment at LO, for the Lepage plot the results are divided by the Yamaguchi results, therefore the spike can be seen at all orders. As expected the slope is increasing by approximately one order by order. The deviant value for N 2 LO τ 2 is due to the accidental zero around 30 MeV (change in the sign), compare to [39].  III. Results for the slopes for the particle-dimer phase shifts δ fitted in the "window of opportunity" for the Yamaguchi model. The uncertainty in the slopes is about 10%. Left for the Lepage plot right for the consistency assessment. The value with the asterisk * is unnaturally large due to a accidental zero, compare to Fig. 2 (right).
We start with the slope fits, using the subtracted propagators τ 1 (k * ) and τ 2 (k * ) for NLO and N 2 LO respectively, we analyze the results for the real part of the particle-dimer phase shift Re[δ]. The plots are shown in Fig. 2. The slopes are increasing order by order as expected, for the Lepage-plots (left) as well as for the consistency assessment (right). The exact value of the increase should be one per order. In the left part of Table III the slopes are shown for the Lepage plot, in the right part the slopes for the consistency assessment. The slopes for other choices of τ i (k * ) are also included. By varying the "window of opportunity" slightly, we estimate the uncertainty in determination of these slopes from the fit at about 10%. Note that the value for η can not be predicted [39], it is determined by the slope for the LO results. 6 It can be seen, that all results agree with the predicted increase approximately. Note that the result for N 2 LO τ 2 (k * ) in the consistency assessment is an exception due to the accidental zero (compare to the discussion in Fig. 2). The values for N 2 LO using τ 1 (k * ) or τ 3 (k * ) are close to the expected value of 5, the corresponding graphs do not exhibit the accidental zero.
Taking into account the 10% uncertainty in the determination of the slope, the results in Table III show that using τ 2 and τ 3 at NLO leads to no significant improvement of the slope compared to τ 1 . This provides a justification for our choice of using τ 1 at NLO. Since we have one more constant, H 2 , at our disposal at N 2 LO, one more subtraction can be accommodated. This motivates our use of τ 2 instead of τ 1 at N 2 LO despite the insignificant improvement in the slope.

slope
LO NLO N 2 LO fit [39] 1.9 2.9 4.8 our fit, no sub. 1.8 our fit, τ 1 2.8 4.6 our fit, τ 2 2.9 6.1* our fit, τ 3 2.8 3.6 Additionally, we have repeated the same analysis for k cot δ instead of the phase shift δ (the observable considered in Ref. [39]). The extracted slopes in Table IV are again consistent with our choice τ 1 (k * ) and τ 2 (k * ) in the calculations at NLO and N 2 LO, respectively.
To summarize, solving the scattering equation for the particle-dimer amplitude in EFT, while treating the spurious pole as proposed above, we have explicitly demonstrated that the numerical solution systematically converges to the exact result, obtained in the Yamaguchi model, which does not contain spurious poles. Moreover, the pattern of this convergence, in general, follows the theoretical predictions. Hence, the theoretical construction of Sect. III has been verified. In the last subsection we have focused on the consistency and model description of the EFT-expansion. We have provided some evidence for our choice τ 1 (k * ) and τ 2 (k * ) in the calculations at NLO and N 2 LO, respectively, based on the behavior of the slopes in Lepage and consistency plots. In this subsection the optimal order of the subtraction polynomial is investigated further, providing additional justification for the choice done earlier. Namely, the numerical calculations discussed in the last sections are repeated for different orders, which means different choices of τ i (k * ) as defined in equation (40). In the left part of Fig. 3 the EFT results at NLO for different τ i (k * ) are compared with the Yamaguchi model for the real part of δ. It becomes clear that τ 1 (k * ) describes the model the best. Further subtractions do not improve the reproduction of the model, they actually make it worse. This means one subtraction seems to be optimal. The right part of Fig. 3 shows the corresponding imaginary part, here a tiny improvement from τ 1 (k * ) to τ 2 (k * ) is visible. However, this is only true for very large values of the momentum p and the improvement is well below the expected EFT accuracy. The results for τ 1 (k * ) agree with the Yamaguchi model everywhere within the EFT uncertainty. There is no improvement from τ 2 (k * ) to τ 3 (k * ) at all. To conclude, at NLO the phase shift can be described most accurately using τ 1 (k * ). This choice is consistent with the slopes for the Lepage plots and the consistency assessments shown in Table III and  Table IV, where all slopes for NLO agree within 10%. Therefore we choose the minimal amount of subtractions in the following, which means using τ 1 (k * ) at NLO. In Fig. 4 the phase shift is shown for N 2 LO. The results appear to be similar to the situation for NLO, the real part is described the best by τ 1 (k * ), while the difference between τ 2,3 (k * ) and the model is larger. However, the effect is small in N 2 LO, as all three choices of τ i (k * ) agree within a power-counting estimation of the EFT uncertainty in the wide interval including the opportunity window. 7 For the imaginary part, however, the increase from τ 1 (k * ) to τ 2 (k * ) is large. The reproduction of the imaginary part of the model is better for τ 2 (k * ) than τ 1 (k * ). But again from τ 2 (k * ) to τ 3 (k * ) no improvement can be seen. Since the differences for the real part are not significant and the imaginary part clearly indicates to use τ 2 (k * ), we choose τ 2 (k * ) for the N 2 LO calculations. The results in the last chapters show a zero for the phase shift δ = 0 around p = 80 MeV for the Yamaguchi model. As discussed above this makes the determination of the slope difficult, and limits the "window of opportunity" to low energies. To test our method for higher values of the window, a different choice for the effective range is investigated. We choose r = 0.8768 fm(= 0.5r) and the same a = 5.4194 fm as before. This moves the zero outside the considered energy region. Besides the unphysical pole is shifted to k 2 = 410.149 MeV. The corresponding Yamaguchi parameters are given by λ = −0.000049 MeV −2 and β = 622.5 MeV. The values for the three-body forces are summarized in Table V. The results for the quantity k cot δ are shown in Fig. 5. The pole around p = 80 MeV is not present anymore. Everything else follows the pattern described for the Yamaguchi model for r = 1.7536 fm. The description becomes better with increasing orders of the EFT.
Since the spike is shifted we are able to choose higher values for the "window of opportunity". We choose it to be between 75 MeV and 125 MeV. In Fig. 6 the Lepage plot is shown. The different orders of the EFT can clearly be distinguished. The slopes increase order by order (see Table VI), but the increase in the Lepage plot from LO to NLO is slightly larger than expected. With the spike shifted, the slopes are stable under small changes to the "window of opportunity". The general pattern of the method behaves as expected and strengthens the assumptions.

G. Gauss model
To further check our results, we perform the same analysis for an additional model potential, namely a Gauss potential. For the Gauss model the regulator is given by Similar to the Yamaguchi model (compare to equation (34)), this leads to (for E < 0) To connect the EFT parameters to the model, we choose λ G to fulfill For the effective range a = 5.4194 fm and r = 1.7536 fm, this results in λ G = 359.134 MeV. In equation (46) the parameter E d is a input value to the Gauss model, so the two parameters λ G and E d are described by the two parameters a and r in the EFT. However E d is equivalent to the position of the root of d −1 G (E) and therefore is the value of a two-body bound state. So for the chosen values of a and r it can be identified with the binding energy of the deuteron, E d ≈ 2.22 MeV. The dimer-propagator τ (q, E) is given by In the numerical calculations we use the un-expanded equation for d G (z) (first line of equation (46)). The one-particle exchange in the Gauss model, Z G (p, q, E), is given by a formula similar to Eq. (38). To avoid numerical difficulties regarding the poles of the angular integral, this is calculated partly analytically and partly numerically. For more details regarding this see the appendix. The values for the three-body forces are fine-tuned to reproduce the Gauss results at p = 0.001 MeV for H 0 and at p = 10 MeV for H 2 . The values can be seen in Table VII.  In Fig. 7 numerical results for the Gauss model and the EFT at different orders are shown. It can be seen, that NLO and N 2 LO are clearly better in describing the Gauss model than LO, with N 2 LO being also better than NLO. In the right part the imaginary part is shown, the EFT results also improve order by order. It is useful to note, that since for the Yamaguchi model in the last chapter and the Gauss model here the parameters λ Y , β and λ G , E d are fine-tuned to give the same a and r, the models both exhibit a pole (zero in δ) around 80 MeV. This results in the same problems for the EFT describing the Gauss as before. The "window of opportunity" is chosen between 42 MeV and 55 MeV.
In the Lepage plot in Fig. 8   results differ, not only is the increase of slope the slope form NLO to N 2 LO larger than expected, 8 also the values are larger than in the Yamaguchi case. This can be explained by the accidental zero at p = 33 MeV. Similar to the results for the Yamaguchi model shown in the consistency assessment in Fig. 2 (right) for N 2 LO using τ 2 (k) the sign of the difference is changing. In the consistency assessment in Fig. 8  Results of the slopes for the real part of the quantity k cot δ for the Gauss model fitted in the "window of opportunity". Left for the Lepage plot right for the consistency assessment. All results for the Lepage plot for N 2 LO, marked by an asterisk, exhibit a accidental zero and are therefore unexpected large, compare to the Fig. 8.
To conclude, the results using the presented method to deal with the unphysical pole k 2 can also be used to describe the Gauss model. The description is improved order by order. The obtained slopes increase as expected as well in the Lepage plot as in the consistency assessment, where the discussed deviations are not caused by the method.

V. SUMMARY AND CONCLUSIONS
In this paper, a novel procedure for removing the contribution from spurious poles in the three-body Faddeev equation for pionless EFT has been proposed. These poles emerge in the two particle scattering amplitudes, which enter the three-body integral equation. Albeit the spurious poles appear below threshold, at energies where the EFT treatment is no more applicable, they still influence the low-energy behavior of the particle-dimer (three-particle) amplitudes. In the three-body integral equation the two-particle amplitudes are evaluated at large negative energies because an integration over all momenta is carried out. Furthermore, the residue of these poles can have either sign, leading to the problems with the three-particle unitarity at low energies.
In the literature, there exist different methods for treating spurious poles. The most popular one is based on a strictly perturbative expansion of the two-body amplitude in the range parameter(s) [14][15][16][17][18]. It will be, however, difficult to use this approach in a finite volume for the extraction of the three-body observables from lattice data. The reason for this is that the expansion diverges in the vicinity of the two-particle energy levels in a finite box, leading to more and more singular expressions at higher orders.
i) In the present paper, we propose a method which enables one to circumvent this problem, expanding only the part of the two-body amplitude that contains spurious poles. Such an expansion can be systematically carried out. Furthermore, in perturbation theory, the counting rules in the underlying EFT are closely linked to the above-mentioned expansion -at a given order in the EFT counting, only first few terms in this expansion should be retained (the number is determined by the order in the EFT expansion). Adding more terms in the expansion does not lead to an increased accuracy. However, due to the non-perturbative character of the three-body integral equation, the above counting can be regarded merely as a rule of thumb, and the optimal number of subtractions should be determined in actual calculations.
ii) The proposal has been tested in numerical calculations in a toy model, using Yamaguchi and Gauss potentials in the two-body sectors. The results of the exact calculations have been confronted with the results, obtained within EFT, matched to the model parameters in the two-and three-body sectors. Moreover, the consistency assessment has been carried out, comparing the EFT results in different orders. In a result of these studies, a clear pattern emerges. The agreement with the exact calculations systematically improves at higher orders. Already at N 2 LO, the exact results are reproduced very well. Moreover, expanding the spurious pole part in the two-body amplitude, it is seen that, after few steps, the accuracy does not further increase when more terms are subtracted. This is fully in line with our expectations. The optimal number of the subtraction terms is slightly lower than the expectation from perturbation theory. This is not entirely surprising, bearing in mind the non-perturbative character of the three-body problem at hand.
iii) It would be extremely interesting to reformulate the three-body quantization condition in a finite volume as given, e.g., in [19][20][21][22][23] along similar lines. We leave this application for a future publication.