1 Introduction

Hadronic interactions at low and intermediate energies are typically characterized by a combination of elementary and composite particle features. While at long distances hadrons behave as elementary particles and their interactions can be described in terms of purely color singlet degrees of freedom, at short distances their composite character becomes manifest in terms of quark and gluon fields in the fundamental and adjoint representations of the color group, respectively. The relevant scale separating between this dual description marks the onset of a confinement scale and we expect it to be of the order of the hadron size, which generally is found to be about 1 fm. While the hadronic dynamics can be organized quite often as a long distance perturbative hierarchy with an increasing number of exchanged particles, it is by itself incomplete; some further either ab initio or phenomenological information reflecting the underlying quark-gluon structure is needed to provide a full description of the scattering process.

The way how this separation is visualized in the complex energy plane is not completely straightforward. Traditionally, and within a genuinely hadronic picture, one appeals to Mandelstam analyticity [1], i.e. the assumption that a scattering amplitude can be expressed by double dispersive integrals in terms of double-spectral density functions, where the integration ranges extend over those regions in the Mandelstam plane where the corresponding double-spectral functions have non-vanishing support [2]. This viewpoint is ultimately grounded in the Mandelstam conjecture, which holds in lowest order in the coupling constant in quantum field theory [1, 3] or to all orders within a non-relativistic context in potential scattering [4], and, which, in the \(\pi \pi \) scattering case, has been rigorously proved in a finite domain [5, 6]. It is noteworthy that under this same assumption an equivalent local and energy dependent optical potential of non-relativistic form was derived many years ago by Cornwall and Ruderman [7, 8]. For a balanced review on these issues at the textbook level see, for e.g., [9, 10]. The existence of a finite analyticity domain suggests in turn the very existence of a finite cut-off on a purely hadronic basis but without an explicit reference to the underlying quark-gluon dynamics and in particular to the confinement scale, so that the cut-off may be determined phenomenologically from data.

Pion-pion scattering is the simplest reaction in QCD mediated by strong interactions involving the lightest hadrons. Tight theoretical constraints based on analyticity, crossing, unitarity, chiral symmetry and Regge behavior can be imposed (see for e.g. [11] for an early review). The machinery of effective field theories (EFT) [12] and in particular its implementation in Chiral Perturbation Theory (\(\chi \)PT) [13] has enabled as a consequence, the most precise theoretical extraction of the \(\pi \pi \) S-wave scattering lengths to date with about an order of magnitude more precision than the experiment [14,15,16,17,18,19,20,21,22], an unprecedented case in strong interactions, where invariably just the opposite situation happens. A historic overview is given in [23]. Along these lines, the most precise \(\pi \pi \)-scattering analyses to date have been obtained in [16, 18, 22]. The latter corresponds to a \(\pi \pi \) description up to \(\sqrt{s}=1.42\) GeV, obtained by fitting the available experimental data from \(\pi N \rightarrow \pi \pi N\) and \(K_{e4}\) decays while imposing as further constraints Roy and Roy-like equations, and with statistical uncertainties satisfying the necessary normality requirements of the residual distributions [24], (see for e.g. [25, 26] for reviews). We stress that despite all these tight mathematical constraints, most of its non-perturbative setup rests upon the validity of the Mandelstam conjecture [1, 3], a result which, as already mentioned, has not yet been rigorously proven since it was first proposed in 1958. This tacit assumption will also be made throughout our work.

In the present paper we invoke the equivalent local and energy dependent optical potential approach suggested long ago in [7, 8] to describe \(\pi \pi \) scattering in coordinate space. In order to do so, we consider a relativistic Schrödinger equation and define a potential to describe the \(\pi \pi \) interaction by matching the field theoretical result to an equivalent quantum mechanical problem in perturbation theory. Phenomenological precursors of \(\pi \pi \) scattering analyses in coordinate space were prompted in [27, 28] within the boundary condition model of strong interactions [29]. Equivalent coordinate space potentials using the Mandelstam representation or the Bethe–Salpeter equation as a starting point were also proposed to all orders in [30,31,32]. As it will become clear below, it is remarkable if not surprising that so little work on \(\pi \pi \) scattering has been conducted within this approach as compared to more popular momentum space methods. Our work fills this gap by implementing Wilsonian ideas inspired by recent developments in the NN case [33,34,35]. These NN investigations had as a consequence a selection of the largest np and pp database up to energies about pion production threshold of \(3 \sigma \) mutually consistent data. Our present investigation within \(\pi \pi \) is in a sense of exploratory character and it pretends also to provide some training playground with an eye put on the more compelling \(\pi N\) case, where the selection of the currently existing database is largely needed (see for e.g. [36,37,38] and references therein).

At short distances, where the interaction is non perturbative, we will assume a complete ignorance of the strong interaction behavior and consider a coarse graining of the interaction instead, very much in the spirit of the work done in [33,34,35] for the NN case. The basic idea is to separate the \(\pi \pi \) interaction into an inner and outer region at a given separation distance, \(r_c\), located at about some elementarity radius. This radius is defined so that at larger distances pions behave effectively as point-like particles. We will assume that in this long distance regime their interactions are ruled by chiral symmetry, and hence they become calculable within \(\chi \)PT. Thus, for \(r>r_c\), we will construct a chiral potential with the correct low-energy analytic properties by matching both quantum mechanical and field theoretical scattering amplitudes in perturbation theory.Footnote 1 On the contrary, the inner region, \(r<r_c\), is regarded as unknown and sampled with the minimal de Broglie wave length determined by the maximum energy we want to describe. This corresponds to a coarse graining of the short range piece and, in its simplest realization, the inner potential will be written as a superposition of equidistant delta-shell interactions. A key issue is to confidently determine the numerical value of the separation scale \(r_c\), since, as noted in [41, 42] and we will see below, the combination \(p r_c \), with p the CM momentum, will fix the total number of independent fitting parameters. The longest range interaction corresponds to a \(2\pi \)-exchange which is \({\mathscr {O}} ( e^{-2 m_\pi r}) \), so that a naive estimate suggests \(r_c \sim 2/(2\,m_\pi )\sim 1.4 \) fm,Footnote 2 a number which will be corroborated by our numerical analysis.

While the potential approach has been explained in great detail in previous works within the NN context (see for instance [42]), it is unconventional within the \(\pi \pi \) scattering folklore. Thus, we will assume no previous knowledge from the side of the reader and for the sake of completeness we will briefly go through all the important issues along the paper. Moreover, \(\pi \pi \) scattering is characterized because at resonance energies relativistic effects cannot be ignored. For instance, for the prominent case of the \(\rho \)-meson \(\sqrt{s}=m_\rho \gg 2 m_\pi \). Unlike the NN case, a new important aspect in the discussion is related to crossing symmetry, which actually intertwines the s, t and u channels.Footnote 3 In addition, the current extraordinary precision achieved theoretically in extracting the S-waves scattering lengths or the lightest \(\pi \pi \) resonance pole parameters [43,44,45,46] provides a great confidence on the theoretical ideas supporting these benchmarking extractions. The fact that the coarse graining approach works for NN scattering in a regime where relativistic and inelastic effects become important, such as pp scattering up to \( \sqrt{s}~\sim 2\) GeV [42], suggests extending the method to other hadronic reactions under similar operating conditions.Footnote 4

Finally, for the sake of completeness let us mention that lattice calculations are naturally formulated in coordinate space. These calculations attack the problem on the finite lattice spacing and the finite volume in two different fashions: either an (energy dependent) potential is determined and the Schrödinger equation is solved subsequently in the continuum or, alternatively, the energy level shifts are determined on the lattice and converted into phase-shifts by means of Lüscher’s formula [47]. Actually, in a pioneering work [48], the \(I=2\) \(\pi \pi \) S-wave scattering phase shifts from Lattice QCD have been determined. Later on, \(\pi \pi \) scattering has been studied from \(N_{\mathrm {f}} = 2+1\) and \(N_{\mathrm {f}} = 2+1+1\) flavors in [49] and [50], respectively. Connected and Disconnected Contractions have also been analyzed in [51]. In addition, the \(\pi \pi I=2\) channel has also been studied within the potential approach [52]. A comparison between potential and Lüscher’s approaches has been undertaken in [53] for the \(I=2\) case, with rather similar results. We remind that both methods have potential drawbacks. On the one hand, the potential method uses interpolating fields which may distort the physics at short distances, and we will explicitly show that in a chiral expansion such potential presents a short distance singularity, which evades the conventional solutions of the Schrödinger equation. On the other hand, the current applicability of this Lüscher’s method [47] requires the interaction to sharply vanish at the edges of the volume (in the relative coordinate), a fact that has been often ignored in momentum space treatments (see for e.g. [54, 55]) but needs to be established for the \(\pi \pi \) case. Our analysis below supports this assumption.

The paper is organized as follows. In Sect. 2, we provide a general and brief field theoretical overview of \(\pi \pi \) scattering to fix our notation in a way that our problem can be easily formulated. In Sect. 3, we show our choice for a quantum mechanical description in terms of a complex local and energy dependent optical potential. We analyze the long-range contributions within \(\chi \)PT in Sect. 4, where an expression for the potential is obtained from the discontinuities of the \(\pi \pi \) scattering amplitude in the t-channel. This requires introducing a short distance cut-off to handle the strong short distance power divergences of the chiral potential, an issue which we discuss at length in Sect. 5. In Sect. 6, we analyze the concept of effective elementarity in order to display in two examples how the elementarity radius depends on the particular process. In Sect. 7, we address the problem of coarse graining \(\pi \pi \) interactions with and without the long-range contributions. The analytical properties of the scattering amplitude and the relation of our approach with the N/D method is discussed in Sect. 8. The implementation of inelasticities within a coarse grained perspective is explained in Sect. 9. We also analyze some aspects concerning low energy constants and the number of parameters in Sect. 10. Finally, in Sect. 11 we summarize our main results and provide some outlook for future work.

2 Formalism for \(\pi \pi \) scattering

We start by summarizing the relevant formulae for \(\pi \pi \) scattering to fix our notation and to provide a proper perspective of our subsequent analysis. A comprehensive presentation at the textbook level can be seen in [11] and also in the lecture [56]. More recent upgrades can be consulted in [25, 26]

2.1 Kinematics

For a pion state \(\varphi _\alpha \) with \(\alpha =\{\pm ,0\}\), the \(\pi _ \alpha (p_1)+ \pi _\beta (p_2) \rightarrow \pi _ \gamma (p_1')+ \pi _ \delta (p_2') \) relativistically invariant scattering amplitude can be written as

$$\begin{aligned} T_{\alpha \beta ; \gamma \delta }= & {} (\varphi _\gamma ^* \cdot \varphi _\delta ^*) (\varphi _\alpha \cdot \varphi _\beta ) A(s,t,u) \nonumber \\&\quad +\, (\varphi _\gamma ^* \cdot \varphi _\alpha ) (\varphi _\delta ^* \cdot \varphi _\beta ) B(s,t,u) \nonumber \\&\quad +\, (\varphi _\delta ^* \cdot \varphi _\alpha ) (\varphi _\beta \cdot \varphi _\gamma ^* ) C(s,t,u) \, , \end{aligned}$$
(1)

with \(s=(p_1+p_2)^2\), \(t=(p_1-p_1')^2\) and \(u=(p_1-p_2')^2\) the standard choice of Mandelstam variables. If we take \(\varphi _{\pm } = ( \phi _1 \pm i \phi _2) /\sqrt{2} \) and \(\varphi _{0}= \phi _3\), with \(\phi _a \cdot \phi _b = \delta _{ab}\), in the Cartesian basis we obtain

$$\begin{aligned} T_{ab;cd}= & {} A(s,t,u) \delta _{ab} \delta _{cd} + B(s,t,u) \delta _{ac} \delta _{bd}\\&\quad +\,C(s,t,u) \delta _{ad} \delta _{bc}\, , \end{aligned}$$

where A(stu) stands for the \(\pi ^+\pi ^-\rightarrow \pi ^0 \pi ^0\) amplitude. This amplitude is the only independent one thanks to isospin, crossing and Bose-Einstein symmetries, \(B(s,t,u)=A(t,s,u)\) and \(C(s,t,u)=A(u,t,s)\). Denoting \(T_{I}(s,t,u)\) as the isospin combination with well defined isospin I in the s-channel, one has

$$\begin{aligned} T_{I=0}(s,t,u)= & {} 3A(s,t,u) + A(t,s,u)+ A(u,t,s) \, ,\nonumber \\ T_{I=1}(s,t,u)= & {} A(t,s,u)- A(u,t,s) \, ,\nonumber \\ T_{I=2}(s,t,u)= & {} A(t,s,u)+ A(u,t,s) \, . \end{aligned}$$
(2)

For the normalization, we will use here the conventions in [57,58,59]. The partial-wave decomposition in the s-channel becomes

$$\begin{aligned} T_{I}(s,t,u)= & {} 16\pi \sum _{J=0}^\infty \left[ 1+(-1)^{J+I} \right] \nonumber \\&\times (2J+1) t_{IJ}(s) P_J (z)\, , \end{aligned}$$
(3)

where \(z = 1 + 2 t/(4m_\pi ^2 -s)\) is the s-channel scattering angle, \(m_\pi =139.57\) MeV the pion mass, \(P_J( \cos \theta )\) the Legendre polynomials and \(t_{IJ}(s)\) is the partial-wave projection of the \(\pi \pi \) scattering amplitude with isospin I and total angular momentum J. Thus, for waves fulfilling the relation \((-1)^{J+I}=1\) one has

$$\begin{aligned} t_{IJ}(s)= & {} \frac{1}{64\pi } \int ^{+1}_{-1}\, \text {d}z; P_J(z)\, T_I\left( s,t(s,z), u(s,z)\right) \nonumber \\= & {} \left( \frac{\eta _{IJ} (s) e^{2\mathrm{i}\delta _{IJ}(s)}-1}{2\mathrm{i}\,\sigma (s)}\right) \, , \end{aligned}$$
(4)

with

$$\begin{aligned} \sigma (s) = \sqrt{1-\frac{4 m_\pi ^2}{s}}\, , \end{aligned}$$
(5)

the \(\pi \pi \) phase factor and \(\delta _{IJ}\) the scattering phase shift. The in-elasticity \(\eta _{IJ} (s)= 1 \) for \(s< 16 m_\pi ^2\) and the unitarity condition for the partial wave amplitude reads in the elastic region

$$\begin{aligned} \mathrm{Im}\,t_{IJ} (s) = \sigma (s)|t_{IJ}(s)|^2\quad \mathrm{for} \quad 4 m_\pi ^2 \le s\le 16 m_\pi ^2. \end{aligned}$$
(6)

Of course, for \(s> 16 m_\pi ^2\) one has absorption \(\eta _{IJ} (s) < 1 \) and inelastic processes such as \(2\pi \rightarrow n \pi \) take place at \(\sqrt{s}=0.56,\,0.84,\,1.12\) and 1.40 GeV for \(n=4,6,8\) and 10, respectively, as well as \(K \bar{K}\) and \(\eta \eta \) at \(\sqrt{s} \sim 1\,\mathrm{GeV}\), etc.

In our discussion we will also use the quantum mechanical amplitude \(f_{IJ}(p)\) defined by

$$\begin{aligned} f_{IJ}(p)= \frac{2}{\sqrt{s}} t_{IJ}(s),\,\qquad s = 4 (p^2+m_\pi ^2) \, , \end{aligned}$$
(7)

with p the CM momentum. For elastic scattering one has \(f_{IJ}(s)^{-1}= p \cot \delta _{IJ} - i p\), so that at low energies one has the threshold expansion

$$\begin{aligned} \mathrm{Re} f_{IJ}(s) = p^{2J} \left[ a_{IJ} + b_{IJ} p^2 + \cdots \right] \, , \end{aligned}$$
(8)

with \(a_{IJ}\) and \(b_{IJ} \) the lowest threshold parameters. An equivalent way of representing the low energy behavior is

$$\begin{aligned} \frac{\tan \delta _{IJ}(s)}{p^{2J+1}}= a_{IJ} + b_{IJ} p^2 + \cdots \, , \end{aligned}$$
(9)

or by an effective range expansion

$$\begin{aligned} p^{2J+1} \cot \delta _{IJ} (s)= - \frac{1}{\alpha _{IJ}} + \frac{1}{2} r_{IJ} p^2 + \dots \, , \end{aligned}$$
(10)

where \(\alpha _{IJ}=-a_{IJ}\) and \(r_{IJ}/2=-b_{IJ}/a_{IJ}^2 \) is the effective range, which is generally positive (see below). Usually, the expansion (9) works for small scattering lengths, such as \(\pi \pi \) whereas (10) works for large scattering lengths, such as NN (see, e.g., [60, 61] for a discussion).

2.2 Anatomy of the \(\pi \pi \) interaction

The purpose of the present paper is to coarse grain the unknown pieces of the \(\pi \pi \) interaction in configuration space. It is thus important to gather some features emerging from comprehensive studies over the last decades [14,15,16,17, 19,20,21,22]. According to these findings the partial wave expansion in (3) is decomposed into two contributions: the low energy contribution described by means of a partial wave (PW) expansion to finite order and the high energy contribution assumed to be given by the leading Regge trajectories,

$$\begin{aligned} T_I = T_I|_\mathrm{PW} + T_I|_\mathrm{Regge} \, , \end{aligned}$$
(11)

which accounts for the long and short distance behavior of the scattering amplitude respectively.

A standard quantum mechanical argument based on the impact parameter provides in the semi-classical limit and for an interaction of finite range \(r_c\), the number of necessary partial waves.Footnote 5 The impact parameter is defined as \(b = L/p\) with p the CM momentum and L the orbital angular momentum, which in our case equals the total angular momentum J. The quantization condition for the angular moment yields \(L \approx \sqrt{J(J+1)} \sim (J+1/2)\) for \(J \gg 1\). For a finite range, the maximal impact parameter where scattering happens is \(b_\mathrm{max} \sim r_c\). Thus, for a maximum CM momentum \(p_\mathrm{max}\), the maximum angular momentum \(J_\mathrm{max}\) for which the phase shift is compatible with zero within uncertainties is

$$\begin{aligned} J_\mathrm{max} + 1/2 \sim p_\mathrm{max} r_c,\quad \mathrm{with} \quad |\delta _{J_\mathrm{max}}|\lesssim \Delta \delta _{J_\mathrm{max}}. \end{aligned}$$
(12)

For a maximum energy \(s_\mathrm{max}=2\) GeV\(^2\), corresponding to \(p_\mathrm{max} \sim 0.7 \) GeV, it was found in [22] that waves beyond \(J_\mathrm{max}=4\) are vanishingly small for \(\pi \pi \) scattering. Therefore, one obtains from (12) a range \(r_c \sim 1.3\) fm. This simple estimate will be explicitly exploited below as an educated guess.

Low energies close to threshold are encoded by the threshold parameters, see (8) and (10). The S-wave scattering lengths are \(\alpha _{00}=-0.3\) fm and \(\alpha _{20}=0.03\) fm whereas for the P-wave we have \(\alpha _{11}= -(0.48\,\mathrm{fm})^3\) [22]. These are unnaturally small numbers compared with our above estimate of the range of the interaction, \(r_c \sim 1.3\,\mathrm{fm}\) and the elementarity radius, \(r_e \sim 1.2\,\mathrm{fm}\) (see the discussion in Sect. 6). While the behavior of the isotensor S-wave resembles a repulsive core, with positive effective range \(r_{20}=131.4\,\mathrm{fm}\), the effective range in the isoscalar S-wave and isovector P-wave are negative, \(r_{00}=-8.08\,\mathrm{fm}\) and \(r_{11}=-5.25\,\mathrm{fm}\), respectively. For S waves the Wigner causality bound [64] (see also [65]) restricts the maximum value of the effective range by the inequality

$$\begin{aligned} r_{I0} \le 2 r_c \left( 1- \frac{r_c}{\alpha _{I0}}+\frac{r_c^2}{3 \alpha _{I0}^2}\right) , \end{aligned}$$
(13)

which for \(I=2\) implies \(r_c \ge 0.95\,\mathrm{fm}\). The positivity of the effective range is not implied by this condition, and is usually violated in the presence of resonances. This requires some unconventional shape for the S-wave potential as we will see.

2.3 Chiral Perturbation theory

The scattering amplitude can be computed perturbatively in Quantum Field Theory and in particular in \(\chi \)PT as a sum of Feynman diagrams in an expansion in 1 / f, with \(f \sim 86\,\mathrm{MeV}\) the pion weak decay constant in the chiral limit. In the partial waves basis the expansion can schematically be written as

$$\begin{aligned} t_{IJ} (s)= t_{IJ}^{(2)} (s) + t_{IJ}^{(4)} (s) + \cdots \end{aligned}$$
(14)

where \(t_{IJ}^{(-2n)}={\mathscr {O}} (f^{-2n})\). To one loop order, they were first computed in [12, 13] and the relevant non-polynomial contributions are reproduced for completeness in appendix 1. Explicit analytical expressions for the corresponding partial wave amplitudes are displayed in [40]. They obey the perturbative unitarity relation

$$\begin{aligned} \mathrm{Im}\,t_{IJ}^{(4)}(s)=\sigma (s)|t_{IJ}^{(2)}(s)|^2 \, , \qquad 4 m_\pi ^2 \le s\le 16 m_\pi ^2 \, . \end{aligned}$$
(15)

At lowest order (LO) in the chiral expansion the threshold parameters are unnaturally small, a fact naturally accommodated by \(\chi \)PT with pions coupled derivatively.

2.4 Unitarization vs Crossing

The requirement of crossing is a fundamental one which stems from the local character of Quantum Field Theories. Chiral Perturbation Theory implements this symmetry at any order in the chiral expansion. The problems with perturbation theory, however, are on the one hand the lack of exact unitarity given by (6) and on the other hand the impossibility of describing outstanding non-perturbative features such as the generation of resonances, which emerge as poles of the scattering amplitude on unphysical Riemann sheets. Within a \(\chi \)PT framework, many methods have been proposed (see for instance [40, 66, 67] and references therein) based on imposing exact unitarity while matching perturbation theory at low energies. Most of them are nothing but algebraic tricks or functional solutions to a set of a priori conditions. As such, unitarization methods are not unique but strongly driven by experimental information, which explains partly their success. The Bethe–Salpeter method discussed at length in [40] preserves an identification of Feynman diagrams but it is not free from field reparameterizations or off-shell ambiguities. In addition, they violate crossing symmetry, which, in general, is only fulfilled order by order, although these violations can be statistically not-significant [68].

In Sects. 4 and 7, we will propose yet a new method based on first defining an equivalent quantum mechanical problem and, more importantly, on coarse graining the interaction. Of course, above the inelastic threshold \(s \ge 16 m_\pi ^2\) one may wonder what condition should be imposed instead of just (6).Footnote 6 We will extend the coarse graining idea to the case with inelasticities.

3 Quantum mechanics Formalism

3.1 Relativistic equation

At the maximum CM energy we will be considering in this work \(s_\mathrm{max}=2\;\mathrm{GeV}^2\), relativity and inelasticities are crucial physical ingredients since firstly \(\sqrt{s_\mathrm{max}} \gg m_\pi \) and secondly we can produce up to \(n=(\sqrt{s_\mathrm{max}}-2 m_\pi )/m_\pi \sim 8\) pions as well as one \(K \bar{K} \) and \(\eta \eta \) pair in the final state. From a field theoretical point of view, this could be solved by using a multichannel Bethe–Salpeter equation for the several \(2\pi \), \(4\pi \), \(6\pi \), \(8 \pi \), \(K \bar{K}\) and \(\eta \eta \) coupled channels, but it would be an extremely difficult task, which has never been accomplished to our knowledge. Even in the simplest elastic case the off-shell ambiguities are present for calculations with a truncated kernel [39, 40]. In order to grasp the nature of the ambiguities, consider for instance the case of \(\pi ^0(p_1) \pi ^0(p_2) \rightarrow \pi ^0(k_1)\pi ^0(k_2)\) scattering, in the elastic regime. The Bethe–Salpeter (BS) equation reads,

$$\begin{aligned} T_P (p,k)&= V_P (p,k) \nonumber \\&\quad + \frac{i}{2} \int \frac{\text {d}^4 q}{(2\pi )^4} V_P (p,q) \Delta (q_+) \Delta (q_-) T_P (q,k) \end{aligned}$$
(16)

where \(p=\frac{p_1-p_2}{2}\), \( k=\frac{k_1-k_2}{2}\), \(P=p_1+p_2\) and \( q_\pm = P/2 \pm q\). \(\Delta (q_\pm ) = 1/(q_\pm ^2-m_\pi ^2+i0^+)\) is the free pion propagator and \(V_P (p,k)\) and \(T_P (p,k)\) stand for the two-particle irreducible kernel or potential and the scattering amplitude, respectively. The factor 1 / 2 comes from the scattering of identical particles.

While the BS equation has been the subject of extensive research for a given potential, the main point of [39, 40] was the flexible interpretation of the BS equation within \(\chi \)PT or more generally within EFT. Indeed, while the potential \(V_P (p,k)\) can be organized as a power series \(V_P (p,k)=V_P^{(2)} (p,k)+V_P^{(4)}(p,k)+ \cdots \) with reference to the same expansion of the scattering amplitude \(T_P (p,k)=T_P^{(2)} (p,k)+T_P^{(4)} (p,k)+\cdots \), it can be done only in a on-shell mass scheme, i.e. for

$$\begin{aligned} T(s,t)= & {} T_P (p,k),\quad p^2=k^2 = \frac{s}{4}-m_\pi ^2,\nonumber \\&P \cdot p = P \cdot k=0 . \end{aligned}$$
(17)

Thus, there is an inherent ambiguity in the definition and form of the potential, which has no consequences perturbatively but becomes relevant in the solution of (16) where the off-shellness enters explicitly. This was mended in [39, 40] by invoking an on-shell scheme, namely considering only on shell intermediate states, i.e. \(q^2 = s/4-m_\pi ^2\) and \( P \cdot q=0 \), so that the on-shell amplitude T(st) depends only on the on-shell potential V(st). Unfortunately, it also gives rise to pathologies in the coupled channel case producing spurious singularities due to an improper treatment of the crossed-channel exchanges [69]. The present paper pretends to address crossed-channel exchanges without invoking the on-shell scheme.

3.2 Invariant mass and equivalent Schrödinger equation

We will follow here the invariant mass formulation [70],Footnote 7 already used for NN scattering with an optical potential [42, 71]. This is the simplest way of retaining relativity without solving a BS equation but with a phenomenological optical potential that we review here for completeness. The idea is to write the total squared mass operator as

$$\begin{aligned} {\mathscr {M}}^2 = P^\mu P_\mu + W , \end{aligned}$$
(18)

where W represents the (invariant) interaction, which can be determined in the CM frame by matching in the non-relativistic limit to a non-relativistic potential \(V(\vec x)\). This yields for \(\pi \pi \) scattering after quantization \({\mathscr {\hat{M}}}^2 = 4(\hat{p}^2+m_\pi ^2) + 4 m_\pi V \), with \(\hat{p} = - i \nabla \). Thus, the relativistic equation can be written as \({\mathscr {\hat{M}}}^2 \Psi = 4 (p^2+m_\pi ^2) \Psi \), with p the CM momentum, i.e. as a non-relativistic Schrödinger equation

$$\begin{aligned} (-\nabla ^2 + m_\pi V ) \Psi = (s/4-m_\pi ^2) \Psi \, . \end{aligned}$$
(19)

This corresponds to the simple rule that one may effectively implement relativity by just promoting the non-relativistic CM momentum to the relativistic CM momentum. This minimal relativity ansatz is as good as the more fundamental one based on the Bethe–Salpeter equation as long as we use scattering data to determine the corresponding potential rather than an ab initio determination (see Ref. [40] for an in-depth discussion).

To take into account the inelasticity within the mass-squared construction, we assume a local and energy-dependent phenomenological potential, \(V(\vec r,s)= {\mathrm{Re}}\, V(\vec r,s) + i\,{\mathrm{Im}}\,V(\vec r,s)\), which could be obtained by fitting inelastic scattering data. Due to causality, the optical potential in the s-channel satisfies a dispersion relation for each CM radial distance r of the form  [7]

(20)

where \(\sqrt{s_0}= 4 m_\pi \) is the first \(4\pi \) inelastic threshold and V(r) is an energy independent component. The complete potential includes also the crossed u channel component.Footnote 8 The simple looking Eq. (19), together with the fixed-r dispersion relation (20), incorporates the necessary physical ingredients present in any theoretical approach: relativity and inelasticity consistent with analyticity.

3.3 Isospin and exchange potential

The incorporation of isospin into the game is straightforward. Rotational, isospin and particle exchange invariance requires the representation of the potential to be given by

$$\begin{aligned} V(r)= & {} \left[ V_A (r) + V_B (r)\,\vec I_1 \cdot \vec I_2 + V_C (r)\,(\vec I_1 \cdot \vec I_2)^2 \right] (1+ {\mathscr {P}}_{12})\nonumber \\= & {} V_\mathrm{D}(r) +V_\mathrm{X}(r), \end{aligned}$$
(21)

where \(V_\mathrm{D}\) and \(V_\mathrm{X}\) stand for the direct and exchange potential pieces, respectively. \({\mathscr {P}}_{12}\) is the particle exchange operator, which implements the Bose-Einstein symmetry and that can be factorized as \({\mathscr {P}}_{12} = {\mathscr {P}}_{x}{\mathscr {P}}_{I}\). Moreover, for states with a well defined total isospin \(\vec I = \vec I_1+ \vec I_2\), we can use the relation \(\vec I_1 \cdot \vec I_2 = I(I+1)/2-2\) with \(I=0,1,2\), so that \( {\mathscr {P}}_{I} = (-1)^I\). In addition, for angular momentum eigenstates \({\mathscr {P}}_{x} = (-1)^J \), so that \({\mathscr {P}}_{12} = (-1)^{I+J}\). Therefore, in the isospin basis the potential can be decomposed as

$$\begin{aligned} V = \sum _{I=0,1,2} P_I V_I (1+ {\mathscr {P}}_{12}), \end{aligned}$$
(22)

where we have introduced the projection operators

$$\begin{aligned} P_0= & {} \frac{1}{3} (\vec I_1 \cdot \vec I_2-1) (\vec I_1 \cdot \vec I_2+1), \nonumber \\ P_1= & {} -\frac{1}{2} (\vec I_1 \cdot \vec I_2-1) (\vec I_1 \cdot \vec I_2+2),\nonumber \\ P_2= & {} \frac{1}{6} (\vec I_1 \cdot \vec I_2+1) (\vec I_1 \cdot \vec I_2+2), \end{aligned}$$
(23)

fulfilling the orthogonality relations \(P_I P_{I'} = \delta _{II'}P_I\).

In the partial wave representation, the exchange symmetry of the potential is preserved by just solving the Schrödinger equation for the direct potential for the allowed IJ channels with \((-1)^{J+I}=1\). In addition, for a spherically symmetric potential we have the usual factorization of the wave function [72]

$$\begin{aligned} \Psi (\vec x) = \frac{u_l(r)}{r} Y_{l m_l} (\hat{x}), \end{aligned}$$
(24)

where \(Y_{l m_l} (\hat{x})\) are the spherical harmonics and \(u_l(r)\) is the reduced wave function, fulfilling the radial Schrödinger equation

$$\begin{aligned} -u_l'' (r) + \left[ U(r) + \frac{l(l+1)}{r^2} \right] u_l(r) = p^2 u_l (r), \end{aligned}$$
(25)

where \(U_I (r)=U_I (\vec r)\) is the central potential with isospin I. This equation is indeed regular at the origin Footnote 9

$$\begin{aligned} u_l (r) \rightarrow r^{l+1} \end{aligned}$$
(26)

and it satisfies the asymptotic scattering condition at infinity

$$\begin{aligned} u_l(r) \rightarrow \sin \left( p r - \frac{l \pi }{2}+ \delta _l \right) . \end{aligned}$$
(27)

Thus, the partial wave expansion for the quantum mechanical scattering amplitude with isospin I in the CM system is defined by:

$$\begin{aligned} f_I(p,\cos \theta )=\sum _{J=0}^\infty (2J+1)P_J(\cos {\theta }) \frac{\eta _{IJ}(s)e^{2i\delta _{IJ}(s)}-1}{2ip}. \end{aligned}$$
(28)

3.4 Inverse scattering problem

Although we will be determining the potentials from fits to phase shifts, it is worth reminding that the inverse scattering problem allows one to determine a local and continuous potential directly from scattering data by solving for each partial wave either the Gelfand-Levitan or Marchenko equations (see for e.g. [73] for a review). It can be shown that for holomorphic S-matrix functions both methods yield the same local potential. While usually the discussion is conducted within a non-relativistic setup, according to our discussion above, the analysis can directly be overtaken and interpreted at the relativistic level.

This inverse scattering approach was adopted in [74], where a holomorphic S-matrix was used to parameterize the scattering data. In that work it was found that the S and P-wave potentials have a range around 0.25 fm with strengths between \(100-200\) GeV. Quite remarkably they also found a barrier in the isoscalar S-wave and a repulsive core in the isotensor S-wave. While this is a very insightful and mathematically rigorous approach, this method requires exact knowledge of the phase shifts at all energies. In practice, a meromorphic function is fitted up to a maximum energy corresponding to a maximum momentum \(p_\mathrm{max}\). As we will see below, this puts in practice a limitation to the resolution \(\Delta r \sim 1/p_\mathrm{max}\) with which the potential V(r) may be determined, so that a suitable coarse graining makes sense.

4 The Chiral \(\pi \pi \) local potential

In this section we outline the perturbative matching procedure between quantum mechanic (QM) and quantum field theory (QFT) calculations in order to determine the local and energy dependent chiral potential. The connection between the QFT and QM scattering amplitudes is given by

$$\begin{aligned} T_I(s,t)=16 \pi \sqrt{s} f_I(p,\cos \theta ). \end{aligned}$$
(29)

The potential appearing in this equation will be determined in perturbation theory. For the quantum mechanical problem we have the Born series

$$\begin{aligned} f(p, \cos \theta )= & {} - \frac{1}{4 \pi } \int {\text {d}^3\vec r\,\, U (r,s) e^{-i \vec {q}\cdot \vec {r}}} \nonumber \\&- \int \text {d}^3 \vec r_1 d^3 \vec r_2 e^{i (\vec p' \cdot \vec r_2 - \vec p \cdot \vec r_1 )} \frac{e^{i p r_{12}}}{r_{12}} U (r_1,s) U(r_2,s) \nonumber \\&+\,\cdots , \end{aligned}$$
(30)

where \(r_{12}=\vert \vec r_1-\vec r_2\vert \), \(\vec p\) and \(\vec p'\) are the initial and final CM momenta, respectively, and \(\vec q=\vec p'-\vec p=2 \vec p\sin {(\theta /2)}\) is the momentum transfer. The potential U(rs) is directly defined from the two-particle irreducible states included in the scattering amplitude. We will define the potential through the t-channel exchanges of the amplitude, so that crossing symmetry will be incorporated exactly when symmetrizing the partial wave expansion.Footnote 10 Moreover, in a coordinate space description, contact terms are irrelevant as long as the field theoretical potential is not extended to the origin \(r=0\) since

$$\begin{aligned} \int {\text {d}q\,\,q^2\hat{P}(q^2)\,\frac{\sin {q r}}{qr}}= 0 \, , \qquad r> r_c > 0, \end{aligned}$$
(31)

with \(\hat{P}(q^2)\) a generic polynomial in \(q^2\). Thus, any polynomial part of the scattering amplitude gives a vanishing contribution to the long-range piece of the potential. Therefore, we will analyze only the effect of pion loop contributions on the t-channel.

4.1 Leading order

We will first discuss the lowest non-trivial order since it provides just contact terms. In the Born approximation, i.e. taking the first term in (30), the scattering amplitude just becomes the Fourier transform of the potential [72]

$$\begin{aligned} f^{B}(p,\cos \theta )= & {} -\frac{1}{4\pi }\int {\text {d}^3\vec r\,\, U(r,s) e^{-i \vec {q}\cdot \vec {r}}} \nonumber \\= & {} -\int \limits _0^\infty {\text {d}r\,\, r^2\,\, U(r,s)\frac{\sin {q r}}{q r}}\, . \end{aligned}$$
(32)

This equation can be inverted to give

$$\begin{aligned} U(r,s) = - 4 \pi \int \frac{\text {d}^3 q}{(2\pi )^3} e^{i \vec {q}\cdot \vec {r}} f_B (\vec q)\,, \end{aligned}$$
(33)

so in the Born approximation, the scattering amplitudes can be related to the potential by:

$$\begin{aligned} T_I(s,t)\big \vert _{B}= -4 \sqrt{s} \int {\text {d}^3\vec r\,\, U_I(r,s) e^{-i \vec {q}\cdot \vec {r}}}\,, \end{aligned}$$
(34)

where \(T_I(s,t)\big \vert _{B}\) denotes the part of the amplitude containing contact terms and t-channel exchange. In the same way, the potential (defined in spatial coordinates) is defined from this part of the amplitude by:

$$\begin{aligned} U_I(r,s)= & {} \frac{-1}{4 \sqrt{s}} \int {\frac{\text {d}^3\vec q}{(2\pi )^3}\,\, e^{i \vec {q}\cdot \vec {r}}} \, \, T_I (s,-\vec q^2)\big \vert _{B} \nonumber \\= & {} \frac{-1}{8 \pi ^2 \sqrt{s}} \int \limits _0^\infty {\text {d}q\,\,q^2\, \, T_I (s,t)\big \vert _{t=-q^2} \frac{\sin {q r}}{qr} }\,. \end{aligned}$$
(35)

Using the \(\chi \)PT lowest order amplitudes [13] we get

$$\begin{aligned} U_0^{(2)} (r,s)= & {} \frac{-1}{4 \sqrt{s}} \frac{m_\pi ^2-2s}{2 f^2} \delta ^{(3)} (\vec r)\,, \nonumber \\ U_1^{(2)} (r,s)= & {} \frac{-1}{4 \sqrt{s}} \frac{4 m_\pi ^2- 2 \nabla ^2-s}{2 f^2} \delta ^{(3)} (\vec r)\,,\nonumber \\ U_2^{(2)} (r,s)= & {} \frac{-1}{4 \sqrt{s}} \frac{s-2 m_\pi ^2}{2 f^2} \delta ^{(3)} (\vec r)\,. \end{aligned}$$
(36)

These algebraic manipulations are purely formal, and in fact is unspecified what is the meaning of solving the wave equation with these highly singular potentials. Already at this level, we can see the need of introducing a regularization.Footnote 11

4.2 Next-to-leading order

The NLO contribution becomes more cumbersome. Firstly, we take the potential to be expanded as

$$\begin{aligned} U_I(r,s)=U_I^{(2)} (r,s) + U_I^{(4)} (r,s) + \cdots \,, \end{aligned}$$
(37)

so we get the matching condition

$$\begin{aligned}&-\frac{T^{(4)}_I (s,t)}{4\sqrt{s}} = \int {\text {d}^3\vec r\,\, U_I^{(4)}(r,s) e^{-i \vec {q}\cdot \vec {r}}} \nonumber \\&\quad -\,\int \text {d}^3\vec r_1 d^3 \vec r_2 e^{i (\vec p' \cdot \vec r_2 - \vec p \cdot \vec r_1 )} \frac{e^{i p r_{12}}}{r_{12}}\nonumber \\&\quad \times U^{(2)} (r_1,s) U^{(2)} (r_2,s) . \end{aligned}$$
(38)

Due to the Dirac delta functions in the LO potential (36), we have a divergence for the real part of \( e^{i p r_{12}}/r_{12} \), albeit it can be absorbed in the real part of the NLO potential \(U^{(4)}(r,s)\). Besides, the non-polynomial pieces in \(T^{(4)}\) amplitude corresponding to the t-channel exchange can generally be written as

$$\begin{aligned} T_I(s,t)\big \vert _\mathrm{2\pi }=P(s,t)J(t)\,, \end{aligned}$$
(39)

where P(st) is a polynomial in both t and s, whose analytical expression can be read from (A3), and J(t) denotes the one-loop \(2\pi \) function. In order to integrate this amplitude, we will take advantage of the analytic structure of the loop function J(t), which is analytic in the whole complex plane but for a cut above \(4m_\pi ^2\) with a discontinuity, \(\mathrm{disc}\, J(t)=2\,i\,\mathrm{Im}\, J(t)=2\pi \,i\sigma (t)/16\pi \), with \(\sigma (t)\) the phase-space factor defined in (5). Thus, up to subtractions one finds the dispersion relation

$$\begin{aligned} J(t)= \frac{t-4m_\pi ^2}{16\pi ^2} \int \limits _{4m_\pi ^2}^{\infty }{\text {d}t^\prime \frac{\sigma (t^\prime )}{(t^\prime -t)(t^\prime -4m_\pi ^2)}} + \mathrm{C.T.}, \end{aligned}$$
(40)

where C.T. is a subtraction constant that can be fixed by setting the value of \(J(4m_\pi ^2) = 1/8\pi ^2 \). Likewise we have

$$\begin{aligned} P(s,t)J(t) = \frac{t-4m_\pi ^2}{16\pi ^2} \int \limits _{4m_\pi ^2}^{\infty }{\text {d}t^\prime \frac{ P(s,t^\prime ) \sigma (t^\prime )}{(t^\prime -t)(t^\prime -4m_\pi ^2)}} + \mathrm{\overline{C.T.}}\,, \end{aligned}$$
(41)

where \(\overline{C.T.}\) is a new subtraction constant ensuring the convergence of the dispersion relation. Thus, taking into account the Yukawa integral

$$\begin{aligned} \int \frac{\text {d}^3 q}{(2\pi )^3} \frac{e^{i \vec q \cdot \vec r }}{q^2 + \mu ^2} = \frac{1}{4\pi } \frac{e^{-\mu r}}{r} \end{aligned}$$
(42)

and the inversion of (38), the NLO potential becomes

$$\begin{aligned} U_I^{(4)}(r,s)=\int \limits _{2 m_\pi }^\infty \text {d}\mu \rho ^I(\mu ,s) \frac{e^{-\mu r}}{r}+ \mathrm{\overline{C.T.}}, \end{aligned}$$
(43)

where \(t=\mu ^2\) and the spectral function \(\rho _I(\mu ,s)\) is defined as

$$\begin{aligned} \rho ^I(\mu ,s)=\frac{-1}{128\pi ^3\sqrt{s}} P_I(s,\mu ^2)(\mu ^2-4m^2)^{1/2} \, , \end{aligned}$$
(44)

with \(P_I(s,\mu ^2)\) polynomials in s and \(\mu ^2\) of second degree (see Appendix 1) and C.T. map into contact terms, which are distributions at the origin and hence vanish elsewhere. All necessary integrals can be obtained from the general integral valid for \(r>0\),

$$\begin{aligned}&\int \limits _{2m_\pi }^{\infty }{\text {d}\mu (\mu ^2-4m_\pi ^2)^{n/2} \frac{e^{-\mu r}}{r}} \nonumber \\&= \frac{ n 2^n m_\pi ^{\frac{n+1}{2}} }{\sqrt{\pi } r^{\frac{n+3}{2}} } \Gamma \left( \frac{n}{2}\right) K_{\frac{n+1}{2}}(2 m_\pi r)\,, \end{aligned}$$
(45)

with \(K_n(x)\) the modified Bessel function of order n and \(\Gamma (z)\) the Euler’s Gamma. Polynomials in \(\mu \) can be generated from derivation with respect to r. The chiral potentials obtained directly from the spectral representation (43), read then

$$\begin{aligned} U_{0}(r,s)= & {} \frac{\left( -23 m_\pi ^5 r^2-200 m_\pi ^3\right) K_1(2 m_\pi r) }{128 \pi ^3 f^4 r^4 \sqrt{s}} \nonumber \\&+\frac{\left( -24 m_\pi ^4 r^2-m_\pi ^2 r^2 s-100 m_\pi ^2\right) K_2(2 m_\pi r) }{32 \pi ^3 f^4 r^5 \sqrt{s}}, \nonumber \\ U_{1}(r,s)= & {} \frac{\left( -13 m_\pi ^5 r^2-40 m_\pi ^3\right) K_1(2 m_\pi r)}{128 \pi ^3 f^4 r^4 \sqrt{s}} \nonumber \\&+\frac{\left( -18 m_\pi ^4 r^2-m_\pi ^2 r^2 s-40 m_\pi ^2\right) K_2(2 m_\pi r)}{64 \pi ^3 f^4 r^5 \sqrt{s}}, \nonumber \\ U_{2}(r,s)= & {} \frac{\left( -17 m_\pi ^5 r^2-80 m_\pi ^3\right) K_1(2 m_\pi r)}{128 \pi ^3 f^4 r^4 \sqrt{s}} \nonumber \\&+\,\frac{\left( -30 m_\pi ^4 r^2+m_\pi ^2 r^2 s-80 m_\pi ^2\right) K_2(2 m_\pi r)}{64 \pi ^3 f^4 r^5 \sqrt{s}}.\nonumber \\ \end{aligned}$$
(46)

From a more general point of view, these potentials play for the \(\pi \pi \)-system the role of relativistic van der Waals interactions (see [75] for a review in the atomic case) and hence display their characteristic features: they are attractive and diverge at short distances as \(\sim 1/( r^7 f^4 \sqrt{s}) \), i.e.

$$\begin{aligned} U_0 (r,s)= & {} -\frac{25}{16 \pi ^3 f^4 r^7 \sqrt{s}} + \cdots \,, \nonumber \\ U_1 (r,s)= & {} -\frac{5}{16 \pi ^3 f^4 r^7 \sqrt{s}} + \cdots \,, \nonumber \\ U_2 (r,s)= & {} -\frac{5}{8 \pi ^3 f^4 r^7 \sqrt{s}} + \cdots \, , \end{aligned}$$
(47)

and have the expected exponentially suppressed long distance behavior \(\sim e^{-2 m_\pi r}\), namely

$$\begin{aligned} U_0 (r,s)= & {} -\frac{23 m_\pi ^{9/2} e^{-2 m_\pi r}}{256 \pi ^{5/2} f^4 r^{5/2} \sqrt{s}} + \cdots \nonumber \\ U_1 (r,s)= & {} -\frac{13 m_\pi ^{9/2} e^{-2 m_\pi r}}{256 \pi ^{5/2} f^4 r^{5/2} \sqrt{s}} + \cdots \nonumber \\ U_2 (r,s)= & {} -\frac{17 m_\pi ^{9/2} e^{-2 m_\pi r}}{256 \pi ^{5/2} f^4 r^{5/2} \sqrt{s}} + \cdots \, . \end{aligned}$$
(48)

In Fig. 1 we show the threshold combination \(V_{I}(r, 4 m_\pi ^2) \equiv m_\pi U_{I}(r, 4 m_\pi ^2)\) and, as we can see, they are attractive at all distances. The energy dependence generates a repulsive effect for increasing values of s, i.e.

$$\begin{aligned} \frac{\partial U_I(r,s)}{\partial s}> 0, \qquad s > 4 m_\pi ^2. \end{aligned}$$
(49)

On the lattice, the energy dependence of the potential is generated from the Nambu-Bethe wave function [53]. As already stated in the introduction, the \(\pi \pi \) potential in the \(I=2\) channel has been computed on the lattice by the HAL QCD collaboration [53, 76] for \(a \approx 0.12\) fm on a \(163 \times 32 \) lattice and with a pion mass \(m_\pi \approx 870\) MeV. For these pion masses the value the chiral potentials in (46) become smaller than for the physical case depicted in Fig. 1. In addition, the HAL QCD lattice potential presents a repulsive core below 0.5 fm. This is a feature one can not obtain using the chiral potentials in (46) which display strong short distance singularities. This fact already suggests that they can not be used at arbitrary short distances. At this point it is worth stressing that both the lattice as well as the present approach based on chiral perturbation theory assume point-like sources, a unrealistic feature. In the next sections we analyze this topic in more detail.

Fig. 1
figure 1

Chiral \(2\pi \) exchange potentials \(V_{IJ}\) at threshold \(\sqrt{s}= 2 m_\pi \) as a function of the distance for \(IJ=00\) (dashed), \(IJ=11\) (full) and \(IJ=20\) (dotted) channels

5 Renormalization

The renormalization of non-perturbative amplitudes is a tricky matter, particularly with the highly power-divergent kernels deduced from \(\chi \)PT (see for e.g. [39, 40] for a discussion within the Bethe–Salpeter framework in momentum space). The chiral potential deduced in coordinate space by a perturbative matching procedure presents an energy dependence. In this section we show how the scattering amplitude stemming from the iteration of the two-pion exchange (TPE) chiral potentials in (46) can be renormalized from a coordinate space point of view if the energy dependence is ignored by taking, say, the threshold value \(\sqrt{s}=2 m_\pi \). Hence, we will implement as renormalization conditions the scattering amplitude at threshold. We will see that, while this is a mathematically viable approach, it fails phenomenologically. Furthermore, the consideration of energy dependence will prevent a sensible non-perturbative renormalization procedure. For large values of the coordinate space cut-off \(r_c \gtrsim 1.2\, \mathrm{fm}\), the results will not be affected by taking either \(U_I(s,r)\) or \(U_I(4 m_\pi ^2,r)\).

5.1 Discussion

One of the advantages of the energy dependent coordinate space representation of the potential is that off-shell and field reparameterization ambiguities manifest as contact interactions at the origin. Thus, they reflect the cut-structure of the amplitude, which is hence unambiguously defined. This is unlike their momentum space counterpart, where both polynomial and cut contributions are treated on equal footing [40].

On the other hand, a difficulty with the chiral \(\pi \pi \) potentials in the previous section is that they become singular at short distances. Hence, the solution of the Schrödinger equation is not well defined in a conventional sense, since the short distance behavior is not dominated by the centrifugal barrier and the regular solution given in (26) is not suitable. An early review on the subject can be found in [77]. Singular potentials are commonplace within EFT and finite solutions exist in a renormalization sense within well specified conditions, as discussed at length in the NN scattering case [78,79,80]. Applications for \(\alpha \alpha \)-scattering [81, 82] and atom-atom scattering [83, 84] are well documented by now (see e.g. [85] for a sucint and pedagogical presentation). While these renormalized solutions represent theoretically a viable solution to the problem, we will consider here a more phenomenological interpretation by introducing a short distance cut-off \(r_c\), whose value reflects short distance effects not taken into account in the derivation of the chiral \(\pi \pi \) potential.Footnote 12 This leaves undefined the short distance dynamics.

The energy dependence of the potential takes into account retardation effects. This can be seen as follows; if the potential is given as a function of the difference of two space-time causally related events \(K(x-x')\) then we have

$$\begin{aligned} \int \limits _0^\infty e^{-\sqrt{s}\,t} \text {d}t = \frac{1}{\sqrt{s}}\,, \end{aligned}$$

explaining the \(1/\sqrt{s}\) factor in (46). In the case of “heavy pions” or very low energies \( p \ll m_\pi \), we can take \(\sqrt{s} \sim 2 m_\pi \) and work within a non-relativistic approximation. A compelling consequence of causality is the verification of dispersion relations in the complex energy plane. We will dedicate Sect. 8 to prove that the right analytical properties hold for the partial wave amplitudes.

The spectral representation (40) suggests the use of a spectral regularization consisting of introducing a cut-off \(\Lambda \) at a given value of \(\mu \), so that the potential at short distances becomes regular. We find that for \(\Lambda > 5 m_\pi \) the regularization quenches the potential for \(r < 1/m_\pi \). We explore this issue in more detail in Sect. 6.3.

5.2 The short distance cut-off and boundary condition regularization

Fig. 2
figure 2

Zero momentum integrated-in wave functions (solid) from large distances using the scattering length as input and the chiral \(2\pi \) exchange potentials \(V_{\pi \pi }\) at threshold \(\sqrt{s}= 2 m_\pi \) as a function of the distance for \(I=0\) (left), \(I=1\) (middle) and \(I=2\) (right). We also draw the asymptotic zero energy wave function (dashed) for comparison

One way to analyze the range of validity of the chiral potentials is to discuss the zero momentum scattering, \(p=0\). Thus, we solve the Schrödinger equation using the asymptotic solutions at zero momentum (we discard the isospin label here),

$$\begin{aligned} u_l (r) = \frac{(2l-1)!!}{r^{l}} - \frac{r^{l+1} }{\alpha _l (2l+1)!!}\,, \end{aligned}$$
(50)

with \(\alpha _l \equiv - \lim _{p \rightarrow 0} \delta _l(p) /p^{2 l+1}\) and integrating inward. The result is illustrated in Fig. 2 where we see that the zero momentum wave function presents oscillations at short distances. This behavior can qualitatively be understood if we write the potential at short distances in the form

$$\begin{aligned} U(r)= & {} - \frac{1}{R^2}\left( \frac{R}{r} \right) ^7, \end{aligned}$$
(51)

where \(R \propto 1/f \) is the van der Waals scale, which in our case and from (47) takes the values \(R = 0.94,\,0.68,\,0.79\) fm for \(I=0,1,2\), respectively. Actually, for \(r \ll R \) the centrifugal barrier, \(l(l+1)/r^2\) , can be neglected and a semi-classical approximation holds since \(\lambda '(r) \ll 1\), where \(\lambda (r)\) is the local wavelength, \(\lambda (r) \equiv 1/k(r)=1/\sqrt{-U(r)}\), with k(r) the local wavenumber. Thus, one finds \(\lambda '(r) = 7/2\,(r/R)^{5/2} \ll 1\) for (51). In this case, the WKB wave function reads [72]

$$\begin{aligned} u(r)\approx & {} \frac{1}{\left[ k(r) \right] ^\frac{1}{4}} \sin \left[ \int k(r) \text {d}r \right] \nonumber \\\rightarrow & {} r^{7/4} \sin \left[ (R/r)^{5/2} + \varphi \right] \ , \end{aligned}$$
(52)

where \(\varphi \) is an arbitrary phase, which is fixed by the long distance solution.

According to the oscillation theorem [72], the number of nodes of the wave function at zero momentum corresponds to the number of bound states. However, as we know, there are no bound states in the \(\pi \pi \) scattering case. Thus, if we want to avoid these spurious solutions, we cannot remove the cut-off completely, but it should be larger than the outmost right node of the wave functions depicted in Fig. 2, which are located at \(r\sim 0.5,\,0.4,\,0.4\) fm for \(I=0,1,2\), respectively. In practice, we will take a cut-off slightly above the van der Waals scale. Note that there is a priori no other reason to discard the chiral potential down to these scales. In fact, phase shifts in any partial wave \(\delta _l(p)\) are convergent when the short distance cut-off goes to zero, provided the scattering length \(\alpha _l \) is fixed. This corresponds to a renormalization program already developed in previous studies  [78,79,80,81,82,83,84,85] that will not be pursued any further here. The upshot of these considerations in the \(\pi \pi \) case is that generally one needs a cut-off \(r_c > 0.5\) fm to prevent the appearance of unphysical bound states generated by the TPE potential. However, for these kind of short distance attractive singularity \(U(r) \propto -1/r^7\) finite \(r_c\) effects are minor in physical observables if the energy dependence of the potential is neglected. Take for instance the effective range defined in (10) and depicted in Fig. 3. As one can see, the chiral potential and the scattering length lead to a finite result at short distances, \(r < R\). Up to minor oscillations, it provides values which differ from the experimental ones. In Sect. 6, we will see that this turns out to be much smaller than the elementarity radius, \(r_e \sim 1.2\) fm. Moreover, there is no finite cut-off \(r_c\) which reproduces the experimental values \(r_{00}=-8.08\,\mathrm{fm}\) and \(r_{20}=131.4\,\mathrm{fm}\).

In the previous discussion the energy dependence of the potential was neglected. On the one hand, if we take the energy dependence into account, we see in Fig. 3 a quite different trend at short distances, namely the effective range is not convergent when the cut-off is removed, i.e. for \(r_c \rightarrow 0\). On the other hand, we also see that for \(r_c \gtrsim 1.2 \mathrm{fm}\) this energy dependence in the chiral potential becomes irrelevant.

Fig. 3
figure 3

Renormalized effective range for the integrated-in large-distance wave functions (solid), as a function of the distance for \(I=0\) (left) and \(I=2\) (right) It is obtained using the scattering length as input and the chiral \(2\pi \) exchange potentials \(V_{\pi \pi }\) at threshold \(\sqrt{s}= 2 m_\pi \). In addition, we also plot the energy dependence case (red-dashed), which leads to a divergent value at \(r \rightarrow 0\)

6 Effective elementarity

As we have seen, even at the perturbative level we must introduce a short distance cut-off, \(r_c\). In this section we elaborate on sensible choices of this cut-off on the light of the onset of effective elementarity and its corresponding radius \(r_e\). This is the scale above which particles interact as if they were pointlike. Thus, they can be taken as elementary, so that a hadronic field can be attached to the particle.

6.1 Hadron sizes and form factors

Hadrons have a finite size, which is usually characterized by their form factors and corresponding radii. Roughly speaking, we expect that for two hadrons of size \(r_1\) and \(r_2\) they will not overlap at relative separations above the mean average distance \((r_1+r_2)/2\), and their interaction will correspond to that of two elementary particles. However, hadronic sizes obtained from, say, electroweak form factors are specific on that particular process, as we will illustrate below. In fact, in order to describe \(\pi \pi \) interactions we are interested on the corresponding effective elementary size, \(r_e\), as seen by the interaction. Unfortunately, without a microscopic calculation, there is no way to know this size in the case of strong interactions. Nonetheless, for our discussion on the relevant scales it is important to estimate a priori the separation distance \(r_c\) between the inner and outer pieces of the potential, as this determination has an impact on the minimal number of fitting parameters. Quite generally, this separation distance should be larger than the elementarity radius, \(r_e \le r_c\), and the optimal choice would be to take both distances equal.

Fig. 4
figure 4

Effective elementarity of the pion. We show the Electrostatic (left panel) and gravitational (right panel) potentials for \(\pi \pi \) interactions and compare the point-like limit (red line) with the finite size case (blue line)

Since the pion size is around 1 fm, it is natural to expect that the inelastic non-perturbative description of the interaction will take place in the region below 2 fm.

In order to illustrate our point, let us consider the implications of elementarity for electric and gravitational interactions between pions. The electromagnetic pion form factor reads

$$\begin{aligned} \langle \pi ^+ (p') | J^\mu (0) | \pi ^+ (p) \rangle = (p'^\mu + p^\mu ) F_\mathrm{em} (q)\,, \end{aligned}$$
(53)

whereas the gravitational form factors \(\Theta _1(q^2)\) and \(\Theta _2(q^2)\) are defined by

$$\begin{aligned}&\langle \pi ^b(p') \mid \Theta ^{\mu \nu }(0) \mid \pi ^a(p) \rangle = \frac{1}{2}{\delta ^{ab}} \nonumber \\&\quad \times \left[ (g^{\mu \nu } q^2- q^\mu q^\nu ) \Theta _1(q^2) +\, 4 P^\mu P^\nu \Theta _2(q^2) \right] \,, \end{aligned}$$
(54)

where \(q=p'-p\), \(P=p'+p\). In the Breit frame, where there is no energy transfer, form factors can be interpreted as the Fourier transform of a density [86]

$$\begin{aligned} F(q)=\int \text {d}^3r\, e^{i q\cdot r}\rho (r)\,, \end{aligned}$$
(55)

so that the charge form factor determines the charge density. The gravitational form factor corresponds to the traceless spin 2 component. In the Breit frame the mass density of the pion is given by the 00 component of the energy momentum tensor in the pion state, since its integrated value yields the total mass of the state (see e.g. the discussion in Ref. [88] and the recent review [89]). We will work in the chiral limit so that the \(\Theta _2(q^2)\) contribution is neglected.

Quite generally, form factors are matrix elements of local operators between hadronic states. For the case of operators with well-defined \(J^{PC}\) quantum numbers, a generalized meson dominance is expected to work in harmony with the high energy behavior deduced from QCD counting rules (see for e.g. [87] for a thorough discussion and comparison with lattice QCD data). For example, in the electromagnetic case, it can be parameterized according to vector meson dominance as,

$$\begin{aligned} F_\mathrm{em}(q)=\int \text {d}^3 r \, e^{i q\cdot r} \, \rho (r) = \frac{m_\rho ^2}{m_\rho ^2+q^2}\,, \end{aligned}$$
(56)

where we have kept only the \(\rho \) meson with \(m_\rho = 0.77\,\mathrm{GeV}\), as we are only interested in the long distance properties.

6.2 Point-like vs extended particles interactions

According to the previous discussion, the electrostatic potential can be written as

$$\begin{aligned} V^{el}_{\pi \pi }(r)= & {} \int \text {d}^3r_ 1d^3r_2 \frac{\rho (r_1)\rho (r_2)}{\vert \vec r_1-\vec r_2 -\vec r\vert }\nonumber \\= & {} \int \frac{\text {d}^3 q}{(2\pi )^3} \frac{4\pi }{q^2}| F_\mathrm{em} (q) |^2 e^{i \vec q \cdot \vec r} \nonumber \\= & {} \frac{1}{r} - e^{-m_\rho r} \left[ \frac{1}{2} m_\rho + \frac{1}{r}\right] \nonumber \\\sim & {} \frac{1}{r} \quad \mathrm{for}\quad r> r_e\,, \end{aligned}$$
(57)

which is depicted in the left panel of Fig. 4 and reflects that for \(r>1.2\sim 1.5\) fm, pions start to interact as expected from point-like particles.Footnote 13 In the gravitational case, a good description is found with the tensor \(f_2(1270)\) meson [87] with the mass scale given by \(m_f = 1.2\,\mathrm{GeV}\), so that

$$\begin{aligned} V^{g}_{\pi \pi }(r)= & {} \frac{1}{r} - e^{-m_f r} \left[ \frac{1}{2} m_f + \frac{1}{r}\right] \nonumber \\\sim & {} \frac{1}{r} \quad \mathrm{for}\quad r> r_e\,. \end{aligned}$$
(58)

In this case, right panel of Fig. 4, the interaction becomes elementary at \(r_e \sim 1\,\mathrm{fm}\), a shorter scale. The previous discussion illustrates our point, namely effective elementarity depends on the particular process.

Before proceeding further, it should be noted that within a more microscopic point of view, for instance a cluster quark model, these are not the only possible contributions to the interactions between pions (see e.g. [90,91,92]). Actually, in a Hartree-Fock approximation they would correspond to the direct interaction (Hartree) term. In addition, one also has the exchange (Fock) term where the quarks inside different pions are interchanged. These terms are genuinely non-local at short distances, so that they are exponentially suppressed at long distances. Therefore, they are expected to contribute below the elementarity radius, \(r_e\), and hence can be regarded as finite size effects.

Fig. 5
figure 5

Ratio of spectral regularized and unregularized TPE chiral potentials \(U_\Lambda ^I(r)/U^I(r) \) for spectral cut-offs of \(\Lambda =1.4\) GeV (solid, black), \(\Lambda =1.12\) GeV (dotted, red) \(\Lambda =0.84\) GeV (dashed, blue) at threshold \(\sqrt{s}= 2 m_\pi \) as a function of the distance for \(I=0\) (left), \(I=1\) (middle) and \(I=2\) (right) channels

6.3 Spectral regularization

While we might estimate the elementarity radius for the chiral potentials from the corresponding folding of “strong” densities, we prefer to analyze instead the effect of introducing a cut-off \(\Lambda \) in the spectral function in (43), i.e.

$$\begin{aligned} U_I^\Lambda (r,s) = \int \limits _{2 m_\pi }^\Lambda \text {d}\mu \rho _I(\mu ,s) \frac{e^{-\mu r}}{r}, \end{aligned}$$
(59)

which corresponds to the two-pion invariant mass spectrum.

In Fig. 5 we show the ratio \(U_\Lambda ^I(r)/U^I(r) \) for spectral cut-offs of \(\Lambda =1.4,\,1.12,\,0.84\) GeV at threshold \(\sqrt{s}= 2 m_\pi \) as a function of the distance. This estimate yields a larger elementarity radius, which becomes \(r_e = 1.2\) fm for the largest spectral cut-off and a quenching of about a half for \(\Lambda \sim m_\rho \). This numerical exercise shows that the naive estimate \(r_e \sim 1/\Lambda \) is numerically rather inaccurate. Of course, one could consider \(\Lambda \) as a fitting parameter in the analysis of \(\pi \pi \) scattering. Nevertheless, in our view this has the disadvantage that a model dependence is introduced by the cut-off procedure above the elementarity radius.

Note also that the larger \(r_e\) the smaller the TPE potential, since \(U(r_e) = {\mathscr {O}} (e^{-2 m_\pi r_e})\), so that one may end up a the situation where a model independent treatment becomes only possible when the chiral contribution is actually vanishingly small. In the next section we will present a different strategy. Overall, our numerical results will confirm this pessimistic expectation.

6.4 Potential separation

The previous discussion suggests that we should decompose the potential for each isospin channel as

$$\begin{aligned} U^I (r) = U_\mathrm{Short}^I (r) \theta (r_c - r ) + U_\mathrm{Long}^I (r) \theta (r - r_c)\,, \end{aligned}$$
(60)

where \(U_\mathrm{Short} (r)\) is a short distance and \(U_\mathrm{Long} (r)\) stands for the long distance contribution. The natural choice is to take for the long-range part the unregularized chiral TPE potential

$$\begin{aligned} U_\mathrm{Long}^I (r) = U_{\chi }^I (r)\,, \end{aligned}$$
(61)

i.e. potential computed in \(\chi \)PT to a given order in the chiral expansion. As we have seen in the previous section, the lowest order effect in the chiral potential comes from the two pion exchange, which gives a contribution that at long distances is \(U_{\chi } (r) = {\mathscr {O}} (e^{- 2 m_\pi r} /f^4)\). For \(r_c \sim 1/m_\pi \) this contribution will in principle play a role. Of course, these most peripheral contributions to the interaction will contain perturbative corrections to all orders in \(m_\pi /f\), which are expected to modify slightly the tail of the potential. Therefore the question is what is the relative importance of the unknown short distance piece and the known long distance contribution. Finally, let us mention that while we will keep the energy dependence already present in the chiral potential we will assume for simplicity that the short distance potential is energy independent, as long as the inelasticity is negligible (see also Sect. 9).

Fig. 6
figure 6

Energy independent fits with the delta-shells potential given in (62) using \(\Delta r=0.3\) fm for the \(\pi \pi \) S0-, P- and S2-wave phase shifts. The uncertainties are those quoted in [22], whereas solid-black, green-dashed, red-dotted and blue dot-dashed lines stand for for the central results for \(r_c=\,0.9,\, 1.2,\, 1.5,\, 1.8\) fm, respectively

7 Coarse graining

The basic notion of coarse graining in the scattering problem was outlined by Avilés as early as 1972 in an insightful but forgotten paper [93] within the context of NN interactions. In this article the potential was effectively represented as a sum of delta-shell potentials. This form has important simplifications and a recent comprehensive mathematical study of this specific case has been carried out [94]. Here we extend the approach to account for the \(\chi \)PT potential tail at long distances (see also Refs. [33,34,35] for a parallel treatment of the NN case).

7.1 Short distance potential

The basic idea is as follows. If we want to describe the two-particle CM wave-functions limited to the range \(\Delta p\), only gross information can be determined in an interval \(\Delta r\), with \(\Delta r\Delta p\sim 1\). Thus, for a maximal \(s_\mathrm{{\max }}=2\) GeV\(^2\), we have \(\Delta p \sim \sqrt{s_\mathrm{{max}}/4-m_\pi ^2}\sim 0.70\) GeV and one obtains \(\Delta r \sim 0.3\) fm. This uncertainty suggests that for a limited energy range the potential only needs to be known in a limited number of points. With this in mind, we consider for \(r< r_c\) the \(\pi \pi \) potential as a sum of a subsequent number of \(\delta \) functions separated by about \(\Delta r\). Thus, up to \(r_c\), it requires to introduce \(N_\delta \) delta shellsFootnote 14

$$\begin{aligned} U_\mathrm{Short} (r)= U_{\Delta r}(r) \equiv \sum _{n}^{N_\delta }U (r_n) \Delta r \delta {(r-r_n)}\, , \end{aligned}$$
(62)

with \(r_n = n \Delta r\). Thus, we have two relevant scales in our setup: the separation distance \(r_c\) and the coarse graining scale \(\Delta r\). While one expects the results and main properties of the potential to be independent of the particular choices of \(r_c\) and \(\Delta r\), this can only happen when accurate information on the interaction is available at all energies. In our case \(\Delta r\) acts as an UV regulator whereas \(r_c\) works as an IR regulator.

We will analyze this problem for four different values of \(r_c\). For convenience, we will set \(r_c=0.9\), 1.2, 1.5 and 1.8 fm. The solution of (25) for the delta-shell potential in (62) is straightforward. One has

$$\begin{aligned} u_{l,n} (r) = {\hat{j}}_l (pr) - \tan \delta _{l,n}\,{\hat{y}}_l (pr)\, , \qquad r_n< r < r_{n+1}\,,\nonumber \\ \end{aligned}$$
(63)

where \(\hat{j}_l(x) = x\,j_l(x)\) and \(\hat{y}_l(x) = x\,y_l(x)\) are the reduced Bessel functions of first and second kind, respectively, and \(\delta _{l,n}\) is the accumulated phase shift. The discontinuity in the logarithmic derivative at \(r=r_n\)

$$\begin{aligned} \frac{u_l' (r_{n}^+)}{u_l(r_{n}^+)} -\frac{u_l' (r_n^-)}{u_l(r_n^-)} =U(r_n)\,\Delta r\,, \end{aligned}$$
(64)

with \(r_n^+ \equiv r_n +0^+\) and \(r_n^- \equiv r_n +0^-\) leads after using the unit Wronskian condition \(\hat{j}_l'(x) \hat{y}_l (x)- \hat{y}_l'(x) \hat{j}_l (x)=-1 \) to the bilinear recursion relation for \(\tan \delta _{l,n}\)

$$\begin{aligned} \tan \delta _{l,n+1}= \frac{A_{l,n} (p) + B_{l,n} (p) \tan \delta _{l,n} }{ C_{l,n} (p) +D_{l,n} (p) \tan \delta _{l,n}}\,, \end{aligned}$$
(65)

with

$$\begin{aligned} A_{l,n} (p)= & {} \Delta r \, U (r_n) j_l (p r_n)^2\,, \nonumber \\ B_{l,n} (p)= & {} \Delta r \, U (r_n) j_l (p r_n) y_l (p r_n)+k\,,\nonumber \\ C_{l,n} (p)= & {} -\Delta r \, U (r_n) j_l (p r_n) y_l (p r_n)+k\,, \nonumber \\ D_{l,n} (p)= & {} \Delta r \,U (r_n) y_l (p r_n)^2\,, \end{aligned}$$
(66)

with the initial and final conditions

$$\begin{aligned} \delta _{l,0}(k)=0\,, \qquad \delta _{l} (k) = \delta _{l,N} (k)\,. \end{aligned}$$
(67)

As shown in [95], these equations can be interpreted as the discrete version of the variable phase equation of Calogero [96], corresponding to the limit \(\Delta r \rightarrow 0\). Of course, the previous equations define an integration method in this limit. We stress that the idea of coarse graining is that if U(r) is determined from data with a maximum CM momentum \(p_\mathrm{max}\), the natural resolution of the problem is \(\Delta r \sim 1/p_\mathrm{max}\), and \(U(r_n)\) are the natural fitting parameters.Footnote 15

7.2 Fitting procedures

In this section we present our numerical results based on standard \(\chi ^2\)-fits. Precise \(\pi \pi \)-scattering phase shifts have been obtained in [22, 24] using Roy equations up to \(\sqrt{s_\mathrm{{max}}}=1.42\) GeV and we will use their tabulated values to determine our fitting parameters. While the standard strategy in \(\pi \pi \) scattering studies has been to use the physical region \(\sqrt{s}> 2 m_\pi \), for reasons to be justified in Sect. 7.6 we also include in the S0 fit the subthreshold region, \(0 < \sqrt{s} \le 2m_\pi \) , also deduced in those Roy-equation analyses [22, 24].

Table 1 Energy independent fit values of the \(\lambda _i \equiv U(r_i) \Delta r\) coefficients in GeV units defined in (68) for each partial wave and value of \(r_c\). For \(r_c=\,0.9,\,1.2,\,1.5\) and 1.8 fm, the corresponding number of delta-shells is 3, 4, 5 and 6, respectively. The maximum fitted energy \(\sqrt{s}\) (GeV) is chosen as the maximum one for which one finds a \(\chi ^2/N_p=1\pm \sqrt{2/N_p}\), with \(N_p\) the number of data points

In order to exemplify this formalism, we will focus just on the lowest \(\pi \pi \) partial-waves, namely the isoscalar and isotensor S waves S0 and S2, respectively, and the P wave. For any isospin channel we will take the potential

$$\begin{aligned} U (r)= & {} \left[ \sum _{n=0}^{N} \lambda _n \delta {(r-r_n)} \right] \theta ( r_c- r) \nonumber \\&\quad +\, U^{(4)}_\chi (r)\theta (r - r_c)\,, \end{aligned}$$
(68)

where \(\lambda _n = U(r_n) \Delta r \). For such a maximal energy this means \(\Delta r =0.3 \) fm. Thus, up to \(r_c=\{0.9,1.2,1.5,1.8\}\) fm, it requires to introduce \(N_\delta =\{3,4,5,6\}\) delta shells.

Anticipating the result, we will carry out the study, first without including the chiral tail, since as we will see in Sect. 7.7, it plays a minor role in the resulting fitting parameters \(\lambda _n\).

7.3 Energy independent coarse graining

To start with, we will assume on purpose sufficiently large separation distances so that the long-distance field theoretical contribution \({\mathscr {O}} (e^{-2 m_\pi r})\) can be safely neglected.

One of the main difficulties of the fitting procedure is that we do not have a priori any idea on what a reasonable values of the parameters \(\lambda _j\) are. To overcome this difficulty and in order to constrain as much as possible the delta-shell coefficient values \(\lambda _j\), we will introduce each delta-shell adiabatically, one by one. In addition, the starting fit values for the most internal delta-shell parameters, \(\lambda _0\) and \(\lambda _1\), will be extracted by fitting the \(\pi \pi \) scattering lengths and slope parameters when only one and two delta shells are considered, respectively. These values will change once the physical energy region is fitted. Nevertheless, this procedure will help us to asses the expected size of the inner delta-shell parameter values. The \(\pi \pi \) S0-, S2- and P-wave phase shift obtained in this way for the different values of \(r_c\) previously chosen are plotted in Fig. 6, whereas the value of the delta-shell parameters are given in Table 1. The statistical uncertainties are smaller than the quoted numbers and hence they are not included in Table 1. The maximum fitted energy \(\sqrt{s}\) for each partial waves is chosen as the largest one for which one finds a \(\chi ^2/d.o.f=1\pm \sqrt{2/d.o.f}\). This provides us an indication that we are not underfitting the data and hence that it defines an admissible boundary region. Note that for the S2 wave this condition is already satisfied for the minimum number of parameters and the maximum number of of data points. In addition, as we will discuss in Sect. 7.6, we also include for the S0 fit the subthreshold region \(m_\pi<s<2m_\pi \).

We can see from these results that with a few free parameters, from 3 to 6 depending on the particular choice of \(r_c\), this formalism already allows one to describe the \(\pi \pi \) scattering in the elastic region. In the case of the S2 wave, where the inelasticities are very small at low energies, it is possible to obtain a perfect description up to the maximum energy at \(s_\mathrm{{max}}=2\) GeV\(^2\). For the P wave one obtains a good description, up to energies around the \(\bar{K}K\) threshold at \(\sqrt{s}=1\) GeV, whereas for the S0-wave the description is limited to the region around 0.9 GeV, where the phase presents a huge rise due to the effect of the inelastic \(f_0(980)\) resonance.

Our results reproduce qualitatively features found in [74] using inverse scattering methods, where by construction the potential is a local and continuous function. In particular, for the S waves we find a short distance barrier in the isoscalar and a repulsive core in the isotensor. Of course, we only get a coarse grained version of those potentials, which befits our idea of a finite resolution sampling.

7.4 Threshold parameters

Our results for the threshold parameters defined by (8) are presented in Table 2. As expected from the quality of the fits, they agree within uncertainties with the results obtained in previous analyses [14,15,16,17, 19,20,21,22], which for completness are also given in Table 3.

Table 2 Scattering lengths and slope parameters for the different fits without TPE with an increasing number of delta shells separated by \(\Delta r=0.3\) fm
Table 3 Scattering length and slope parameter values for the \(IJ=00,20,11\) from [16] and [22]

7.5 Resonance poles

Resonances are determined by looking for poles in the second Riemann sheet, which is defined by the complex wavenumber \(k_R= k_r + i k_i \) with \(k_r \equiv \mathrm{Re}\,k_R > 0\) and \(k_i \equiv \mathrm{Im}\,k_R <0\). This can be achieved directly by substituting \(k \rightarrow k_R\) in the Schrödinger Eq. (19), or its coarse-grained delta-shell implementation (65) and (66).

Our numerical results are presented in Table 4. The numerical values are slightly different from those quoted in benchmarking studies based on Roy equations [16, 43, 44], which for completeness are also given in Table 5. This is not fully surprising since the analytic structure of our scattering amplitude are not the same as in the analyses based on the Roy equations, but only to leading order in the chiral expansion and hence the extension to the complex plane is not determined from the phase-shift analysis only. Nonetheless, we find that the result is encouraging and we expect to return in the future in order to implement higher orders to analyze the effect.

Table 4 \(f_0(500)\) and \(\rho (770)\) resonance poles obtained from the different fits without TPE with an increasing number of delta shells separated by \(\Delta r=0.3\) fm
Table 5 \(f_0(500)\) and \(\rho (770)\) resonance poles obtained in [16, 43, 44] using Roy equations

7.6 Crossing and comparison with Roy equations

As we have mentioned in the introduction, there are no direct \(\pi \pi \) scattering data. The closest thing to it are possibly the outcome of Roy equations, which incorporate by construction crossing, analyticity, unitarity and Regge behavior [16, 17, 22]. Moreover, a proper identification of the LECs requires by definition the satisfaction of crossing, particularly in the sub-threshold region. For definiteness, we show in Fig. 7 the real part of the amplitude starting at the edge of the left hand cut \(s=0\), covering the subthreshold region \(0 < \sqrt{s} \le 2 m_\pi \) (where the amplitude is purely real) and the physical elastic region, \(2 m_\pi<\sqrt{s}<2 m_K\), where the imaginary part of the amplitude is determined from unitarity and the corresponding real part. Note that we are neglecting possible inelasticities coming from multi-pion channels. In particular, one can clearly see the Adler zeros of the S0 and S2-wave amplitudes, a characteristic features of the subthreshold region and a distinct fingerprint of chiral symmetry [11].

Fig. 7
figure 7

Real \(\pi \pi \) partial-wave amplitudes as a function of the CM energy, starting from the subthreshold value \(\sqrt{s}=0\) for \(I=0\) (top), \(I=1\) (middle) and \(I=2\) (bottom). We compare the different delta-shell results with the results of the Roy-equation analysis (gray band)

The comparison in Fig. 7 supports the view that fits to phase shifts in the physical region, even to relatively high energies, do not constrain the subthreshold region. This has a large effect on the location of the \(\sigma \)-resonance and the implications will be discussed in more detail elsewhere. In addition, our fits have not imposed crossing correlations and the consideration or not of the subthreshold region in the S0 channel can be grasped by comparing the location of the corresponding resonance. If the subthreshold region \(m_\pi<s<2m_\pi \) is not imposed in the fit, one gets substantially smaller values for the real part. For instance, for \(r_c=0.9\) fm and \(N_\delta =3\) delta-shells one gets \(\sqrt{s_\sigma }= (324 - i\,194 )\) MeV, whereas the pole moves to \(\sqrt{s_\sigma }= (453 - i\,212 )\) MeV when the subthreshold region is also fitted. This large influence of the subthreshold region in the S0 wave is not surprising, since the pole is rather far from the real axis. In contrast, the location of the \(\rho \)-resonance pole proves rather insensitive to the subthreshold region in the P-wave. Although we are not fitting the subthreshold region in this channel, one gets for \(r_c=0.9\) fm and \(N_\delta =3\) delta-shells \(\sqrt{s_\rho }= (762 - i\,70)\) MeV, even when the subthreshold extrapolation of the fit, middle panel of Fig. 7, shows a significant discrepancy with the outcome of the Roy equations. Finally, we also report a discrepancy in the isotensor channel, where the phase shifts are compatible up to \(s_\mathrm{max}\). Actually, in all cases the amplitude has a zero at \(s=0\), due to the factor \(\sqrt{s}\) in the definition of \(t_{IJ}(s)\). Thus, the failure to satisfy the subthreshold behavior in the S2 wave is intrinsic.Footnote 16 We have checked that this does not improve by incorporating the TPE tail of the chiral potential (see also next subsection). The fact that the subthreshold region is quantitatively as relevant as the physical region, is another manifestation of the relevance of crossing in \(\pi \pi \) scattering, and calls for improvement when all relevant partial waves are considered.

7.7 Full potential: Inclusion of two pion exchange

Once the chiral tail has been defined, we can finally construct a full potential. At short distances (\(r< r_c\)), the non-perturbative regime of QCD is encoded by a sum of delta-shells separated by the minimum De Broglie wave-length considered in the analysis. At large distances (\(r> r_c\)), the chiral potential constructed in (46) with the correct analytical properties and left-hand cut contribution is included. Higher order corrections in the chiral expansion are \({\mathscr {O}} (f^{-6})\) and the determination of the corresponding potential requires going in the quantum mechanical picture to second order in the Born approximation. Their analysis is left for future research.

The resulting fits are reported in Table 6. As one can see, and with the exception of the S2 wave with three delta shells, we do not find a substantial improvement due to the explicit incorporation of the TPE potential. A similar trend is found for the corresponding low energy threshold and resonance parameters, see Tables 7 and 8. In any case, chiral effects due to explicit TPE are generally found to play a minor role.

In our analysis we have restricted to the lowest S- and P-waves. The implementation of higher partial waves is straightforward and requires introducing further short range contributions. Due to the centrifugal barrier term, we expect that the number of grid points will be reduced as the angular momentum increases.

Table 6 Energy independent fit including the chiral tail for \(r>r_c\). Values of the \(\lambda _i\) coefficients in GeV units defined in (68) for each partial wave and value of \(r_c\). For \(r_c=\,0.9,\,1.2\,,1.5\) and 1.8 fm, the corresponding number of delta-shells is 3, 4, 5 and 6, respectively. The maximum fitted energy \(\sqrt{s}\) (GeV) is chosen as the maximum one for which the \(\chi ^2/N_p\) is around 1, with \(N_p\) the number of data points

8 Analytic properties and the N/D method

Traditionally, much of the discussion of \(\pi \pi \) scattering has been marked by analyticity and dispersion relations, which at the partial-wave level and in the unsubtracted case read

$$\begin{aligned} t_{IJ}(s)= & {} \frac{1}{\pi } \int \limits _{-\infty }^0 \text {d}s' \frac{\mathrm{Im}\,t_{IJ}(s')}{s'-s-i 0^+}+ \sum _{n} \frac{ g_n}{s_n-s} \nonumber \\&\quad +\,\frac{1}{\pi } \int \limits _{4 m_\pi ^2}^\infty \text {d}s' \frac{\mathrm{Im}\,t_{IJ}(s')}{s'-s-i 0^+}, \end{aligned}$$
(69)
Table 7 Scattering lengths and slope parameters for the different fits including the TPE chiral tail with an increasing number of delta shells separated by r = 0.3 fm
Table 8 \(f_0(500)\) and \(\rho (770)\) resonance poles obtained from the different fits including the TPE chiral tail with an increasing number of delta shells separated by \(\Delta r=0.3\) fm

where \(0< s_n < 4 m_\pi ^2\) are the possible bound states, which we keep for generality, and \(\mathrm{Im}\,t_{IJ}(s) = \rho (s) |t_{IJ}(s)|^2\) for elastic scattering. In this section, we will discuss this issue within the context of our coordinate space framework and for the specific chiral potential derived in Sect. 4, but the results are general. Actually, we will see that the subtractions are not explicitly needed, although they are somewhat encoded into the short distance component of the potential. While many of the issues discussed here have been known for potential scattering since many years [97, 98] (see also [9, 10]), we feel it is necessary to review the main aspects for completeness.

The analytical properties of the scattering amplitude in the complex energy plane are determined by the long-distance behavior of the interaction. This is the basis for dispersion relations, which own their popularity to their link to axiomatic field theory and its straightforward implementation in terms of leading singularities. A frequent method, which has been used in these regard to implement known analytical properties, is the so-called N/D approach. In the N/D approach the partial wave amplitude is written in the form [99]

$$\begin{aligned} t_{IJ}(s) = \frac{N_{IJ}(s)}{D_{IJ}(s)}, \end{aligned}$$
(70)

with \(N_{IJ}(s)\) having only left-hand cut singularities and \(D_{IJ}(s)\) having only right-hand cut singularities. This method has often been invoked in \(\pi \pi \) scattering and implemented in several approximations (see e.g. [100]).

8.1 Jost functions

The way of realizing the N/D representation in potential scattering is well-known. We will describe here the formalism within the coarse grained approach for completeness as well as to provide the \(N_l(k)\) and \(D_l(k)\) functions explicitly, (we drop the isospin index for simplicity). The discussion is naturally carried in terms of the quantum mechanical scattering amplitude defined in (7), \(f_l(k)= 2 t_{l} (s)/\sqrt{s}\), as a function of the CM momentum, k, with \(s=4 (k^2+m_\pi ^2)\). Note that s is invariant under \(k \rightarrow -k \) and hence two-valued. The discussion below entitles to take \(\mathrm{Im} k >0\) for the first Riemann sheet and \(\mathrm{Im} k <0\) for the second.

For a regular potential, the Jost functions, \(F_l(k)\), are defined by the regular solutions of the wave equation at short distances, i.e. \(u_l(r) \rightarrow \hat{j}_l ( kr)\) and subjected to the asymptotic condition at \(r \rightarrow \infty \) [72]

$$\begin{aligned} u_{k,l}(r) \rightarrow \frac{1}{2} \left[ F_l(k)\,\hat{h}_l^{(2)} (kr) + F_l(-k)\,\hat{h}_l^{(1)} (kr) \right] , \end{aligned}$$
(71)

where \(\hat{h}_l^{(1,2)} (x) = \hat{j}_l(x) \pm i \hat{y}_l(x) \) are the reduced Hankel functions fulfilling \(\hat{h}_l^{(1)} (x)^* = \hat{h}_l^{(2)} (x) \) and \(\hat{h}_l^{(1)} (-x) = (-1)^{l+1} \hat{h}_l^{(2)}(x) \). Thus, the S-matrix is then defined as

$$\begin{aligned} S_l(k) = \frac{F_l(-k)}{F_l(k)} = e^{2 i \delta _l(k)}, \end{aligned}$$
(72)

so that the scattering amplitude becomes

$$\begin{aligned} f_l(k) = \frac{F_l(-k)-F_l(k)}{2 i k F_l(k)}. \end{aligned}$$
(73)

The Jost functions fulfill the reflection conditions

$$\begin{aligned} F_{l}(-k^*) = F_{l}(k)^* \end{aligned}$$
(74)

for complex k. Furthermore, due to (74), it can be shown [9, 10] that the functions

$$\begin{aligned} n_l(k) = \frac{F_l(-k)-F_l(k)}{2 i k}\quad \mathrm{with}\quad d_l(k) = F_l(k) \end{aligned}$$
(75)

fulfill the relations

$$\begin{aligned} n_l(k)^* = n_l(-k^*), \qquad n_l(k) = n_l(-k), \end{aligned}$$
(76)

which means that \(n_l(k)\) is purely real for \(\text {Im}\,\,k=0\). Moreover, \(\mathrm{Re}\, n_l(k)=\mathrm{Re}\,n_l(-k^*)\) and \( \mathrm{Im}\, n_l(k) = -\mathrm{Im}\, n_l(-k^*)\), so that \(n_l(k)\) is purely imaginary for \(\mathrm{Re}\,k=0\). In addition,

$$\begin{aligned} \lim _{k \rightarrow \infty } n_l (k) = 0,\qquad \lim _{k \rightarrow \infty } d_l (k) = 1. \end{aligned}$$
(77)

Thus, they have the desired properties for a potential constructed as a superposition of Yukawa potentials. A straightforward consequence of these properties is Levinson’s theorem. While the proofs of these statements have been known for a long time [97, 98], to our knowledge they have not been considered within the present context. We will review them adapted to our complete potential, which can be decomposed into two pieces: a cut-off potential and a Yukawa superposition. While they are discussed separately in the literature, our case at hand involves both cases at the same time. We will show next the pertinent steps.

8.2 Analytic properties

The scattering amplitude obtained perturbatively from \(\chi \)PT enjoys analytical properties deduced from the corresponding Feynman diagrams, i.e. particle exchange generated at the partial-wave level a left-hand cut, whose discontinuity has been used to reconstruct the chiral potential. A relevant question is whether our resulting full amplitudes obtained by solving the Schrödinger equation, share these properties beyond first order in perturbation theory, i.e. whether they satisfy, as expected, dispersion relations.

In order to check the analytical properties of the quantum mechanical amplitude, we note that our chiral potential is indeed a superposition of Yukawa potentials with the exception of the explicit s-dependence in the spectral function. It is convenient to write the potential as a superposition of exponentials in the form

$$\begin{aligned} U_I(r,s) = \int \limits _{2 m_\pi }^\infty \text {d}\mu \sigma _I (s,\mu ) e^{-\mu r}, \end{aligned}$$
(78)

where

$$\begin{aligned} \sigma _I(s,\mu ) = \int \limits _{2 m_\pi }^\mu \text {d}\mu ' \rho _I (s,\mu ') \end{aligned}$$
(79)

and \(\rho _I\) was defined in (44). This is indeed our case for the TPE potential in \(\pi \pi \) scattering, as one can see comparing (43) and (78). Being a Laplace transformation, it corresponds to the so-called analytical potential in the complex-r plane for \(\mathrm{Re}\, r >0\) [101], which fulfills the relation \(\lim _{\rho \rightarrow 0}\,U_I( \rho e^{i \theta },s)=0\) for \(-\pi /2< \theta < \pi /2\), with \(r=\rho e^{i \theta }\).

As we will see below, while the exponential falloff of the potential at long distances suffices to prove the analyticity in the strip \( | \mathrm{Im}\,k | < m_\pi \), with \(k=\sqrt{s/4-m_\pi ^2}\), the spectral decomposition is indeed needed to establish the cut along the line \(\mathrm{Re}\,k=0\) and \( m_\pi < \mathrm{Im}\,k \) of the S-matrix. There are several versions of the proof. On the one hand, it can be proved by estimating a bound for the Jost function using directly the spectral representation. On the other hand, by profiting from the analytical character of the potential and deforming the integration in r into the complex plane so that \(k r >0 \). For completeness we will review next and in a sketchy fashion the second method, as it does not require a bounded spectral function or the use of an spectral regularization (see Sect. 6.3).

The determination of the analyticity domain of the quantum mechanical problem is based on the equivalent integral equation for the Jost Functions as follows

$$\begin{aligned} u_{k,l} (r) = \hat{j}_l(kr) + \int \limits _0^r \text {d}r' {\mathscr {K}}_{k,l}(r,r') U(r',s) u_{k,l} (r'), \end{aligned}$$
(80)

which is a Volterra type integral equation and the kernel is given by

$$\begin{aligned} {\mathscr {K}}_{k,l}(r,r') = \frac{i}{2k} \left[ h_l^{(1)} (kr') h_l^{(2)} (kr) - h_l^{(1)} (kr) h_l^{(2)} (kr') \right] \,.\nonumber \\ \end{aligned}$$
(81)

Taking the large-r limit and comparing with (71), we get

$$\begin{aligned} F_l(k) = 1 + \frac{i}{k} \int \limits _0^\infty \text {d}r\,\hat{h}^{(1)}_l (kr) U(r,s) u_{k,l} (r). \end{aligned}$$
(82)

This equation is the basis for the analytical continuation to the complex-k plane. Actually, in the limit \(r \rightarrow \infty \) one has from (71) that \( |h^{(1)}_l (kr) | = {\mathscr {O}} (e^{-\mathrm{Im}\,k\,r})\) for \(\mathrm{Im}\,k> 0\). Thus, even when \(|u_{k,l} (r)| ={\mathscr {O}} (e^{\mathrm{Im }\,k\,r})\), one has a finite integral provided that the potential goes to zero. As a consequence, \(F_l(k)\) is analytical for \(\mathrm{Im }\,k> 0\). On the contrary, for \(- m_\pi< \mathrm{Im}\,k< 0\) one has \(|h^{(1)}_l (k\,r)\,U(r,s)\,u_{k,l} (r) | = {\mathscr {O}} (e^{-2\mathrm{Im}\,k\, r- 2 m_\pi r} )\) which is convergent for \(\mathrm{Im k}> -m_\pi \). Therefore, \(F_l(k)\) is analytical for \(\mathrm{Im k}> -m_\pi \) and hence \(F_l(-k)\) is analytical for \(\mathrm{Im }\,k < m_\pi \). In conclusion, \(S_l(k)= F_l(-k)/F_l(k)\) is analytical in the strip \( | \mathrm{Im}\,k | < m_\pi \).

For \(k=|k| e^{i\theta }\), due to the analytical character of U(rs), we can deform the contour to \(r=\rho e^{-i\theta }\), so that, for any term in the spectral integral over \(\mu \), one has a pole at \(k = -i \mu /2\). In the case of \(F_l(k) \) and \(F_l(-k) \), this pole becomes a cut after \(\mu \)-integration along the lines \(\mathrm{Re}\,k=0\) and \(-\infty< \mathrm{Im}\,k < -m_\pi \) and \( m_\pi< \mathrm{Im}\,k < \infty \), respectively. The fact that \(F_l(k)\) is analytical for \(\mathrm{Im}\,k>0\) and \(\lim _{k\rightarrow \infty } F_l(k) \rightarrow 1\) allows one to write a dispersion relation in the upper half-circle

$$\begin{aligned} F_l(k)=1 + \frac{1}{\pi } \int \limits _0^\infty \text {d}k' \frac{k'}{k'^2-k^2} \mathrm{Im} F_l(k'), \end{aligned}$$
(83)

where the antisymmetry of \(\mathrm{Im} F_l(-k)=-\mathrm{Im} F_l(-k)\) from (74) has been used. Thus, passing to the variable \(s=4(k^2+m_\pi ^2)\), which is one-valued in \(\mathrm{Im} k >0\), one can define the functionFootnote 17

$$\begin{aligned} t_{l} (s) = \frac{N_l(s)}{D_l(s)}, \end{aligned}$$
(84)

so that we can identify \(D_l(s)=F_l(k)\) for \(\mathrm{Im} k > 0\)

$$\begin{aligned} D_l(s) \equiv 1 + \frac{1}{\pi } \int \limits _{4 m_\pi ^2}^\infty \text {d}s' \frac{\mathrm{Im} D_l(s')}{s'-s-i 0^+}, \end{aligned}$$
(85)

where \( \mathrm{Im}\, D_l(s') = \mathrm{Im}\,F_l(s') \) has a right-hand cut in \(4 m_\pi ^2<s<\infty \) and is real for \(s<4 m_\pi ^2\). Furthermore, in the elastic approximation, using (6) we get

$$\begin{aligned} \mathrm{Im} D_l(s)= - \sigma (s) N_l(s) \, , \qquad s > 4 m_\pi ^2 \, , \end{aligned}$$
(86)

where

$$\begin{aligned} N_{l} (s) = \frac{F_l(-k)-F_l(k)}{2 i \sigma (s)} \, , \end{aligned}$$
(87)

is real for \(s> 4 m_\pi ^2\). For real \(s <0\) we get that \(s+i0^+ \leftrightarrow 0^+ + i \kappa \) with \(\kappa > m_\pi \) and hence

$$\begin{aligned} N_{l} (s+i0^+)= & {} \frac{F_l(-0^+ - i \kappa )-F_l(0^+ + i \kappa )}{2 i \sigma (s)} \, , \end{aligned}$$
(88)
$$\begin{aligned} N_l(s-i0^+)= & {} \frac{F_l(-0^+ + i \kappa )-F_l(0^+ - i \kappa )}{2 i \sigma (s)} \, . \end{aligned}$$
(89)

Thus, the discontinuity is

$$\begin{aligned} \mathrm{Disc}\,N_{l} (s) = 2 i\,\mathrm{Im}\,N_l(s) = \frac{\mathrm{Im}\,F_l(-0^+ - i \kappa )}{i\,\sigma (s)}, \end{aligned}$$
(90)

where we have used that for \(\mathrm{Im} k>0\), \(F_l(k)\) is analytical and hence \(F_l(0^+ + i \kappa )=F_l(-0^+ + i \kappa )\). Finally, note that the s-dependence appearing in the spectral function does not spoil these analytic properties. This completes the proof that the scattering amplitude obtained by solving the Schrödinger equation for the potential in (60) has the correct analytical properties.

8.3 Coarse graining

In order to calculate the Jost functions, in practice we coarse grain the interaction using the delta-shells and we proceed as before. In the discretized version we define the accumulated Jost functions \(F_{l,n}(k) \) and \(F_{l,n}(-k)\) as

$$\begin{aligned} u_{l,n}(r)= & {} \frac{1}{2} \left[ F_{l,n}(k) \hat{h}_l^{(2)} (kr) + F_{l,n}(-k) \hat{h}_l^{(1)} (kr) \right] , \nonumber \\&r_n< r < r_{n+1} \, . \end{aligned}$$
(91)

However, we keep track of both the continuity of the wave functions and the discontinuity of the derivative separately,

$$\begin{aligned} u_l (r_{n}^+) -u_l (r_{n}^-)= & {} 0, \nonumber \\ u_l' (r_{n}^+)- u_l' (r_n^-)= & {} U(r_n) u_l (r_{n}) \,\Delta r\,, \end{aligned}$$
(92)

with \(r_n^+ \equiv r_n +0^+\) and \(r_n^- \equiv r_n +0^-\). Again we use the fact that the Wronskian \(2i ( \hat{h}_l^{(1)} (x) \hat{h}_l^{(2)\,\prime }(x)- \hat{h}_l^{(1)\,\prime } (x) \hat{h}_l^{(2)} (x))= 1\) so that

$$\begin{aligned} F_{l,n+1}(k)= & {} A_{l,n}(k) F_{l,n} (k) + B_{l,n}(k) F_{l,n} (-k), \nonumber \\ F_{l,n+1}(-k)= & {} C_{l,n}(k) F_{l,n} (k) + D_{l,n}(k) F_{l,n} (-k), \end{aligned}$$
(93)

where we have introduced the coefficients,

$$\begin{aligned} A_{l,n}(k)= & {} 1+\frac{i \Delta r \, U(r_n) h_l^{(1)}(k r_n) h_l^{(2)}(k r_n)}{2 k}\,, \nonumber \\ B_{l,n}(k)= & {} \frac{i \Delta r \, U(r_n) h_l^{(1)}(k r_n){}^2}{2 k}\,, \nonumber \\ C_{l,n}(k)= & {} -\frac{i \Delta r \, U(r_n) h_l^{(2)}(k r_n){}^2}{2 k}\,, \nonumber \\ D_{l,n}(k)= & {} 1-\frac{i \Delta r \, U(r_n) h_l^{(1)}(k r_n) h_l^{(2)}(k r_n)}{2 k} \end{aligned}$$
(94)

and the initial conditions are

$$\begin{aligned} F_{l,0}(k) = 1, \quad F_{l,0}(-k) = 1, \end{aligned}$$
(95)

which corresponds to take the normalization \(u_{l,0}(r) = \hat{j}_l(k,r)\). The final values are

$$\begin{aligned} F_{l}(k) \equiv F_{l,N}(k) \, , \qquad F_{l}(-k) \equiv F_{l,N}(-k). \end{aligned}$$
(96)

For illustration purposes, we plot in Fig. 8 the functions \(N_{IJ}(s)\) (which is real) and \(D_{IJ}(s)\) (which is complex) separately above the threshold for the five delta-shell case. Of course, these functions reproduce the phase shifts presented before. As expected, \(\mathrm{Im}\,D_{IJ} (4m_\pi ^2)=0\). In this representation the Breit-Wigner position of the resonance is given by computing the zeros of \(\mathrm{Re}\,D_{IJ}(s_R)=0\). While it is not shown in the pictures, the Jost functions display some oscillatory behavior at higher energies due to the explicit delta-shells. Nevertheless, they still go to the expected values \(D(\infty )=1\) and \(N(\infty )=0\).

Fig. 8
figure 8

\(N_{JI}\) and \(D_{JI}\) (dimensionless) functions appearing in the N/D method (see main text) as a function of the CM energy for \(I=0\) (top), \(I=1\) (middle) and \(I=2\) (bottom)

9 Coarse graining inelasticities

9.1 Energy dependent coarse graining

As mentioned in the introduction, at sufficiently high energies particles are produced and elastic scattering happens in the presence of absorption, which one may view as a leak or hole in the probability. The inelastic region is characterized by the recombination of the two-pion internal configuration, which, of course, demands a complex energy-dependent potential. The complexification of the potential can be understood in terms of the loss of probability of the elastic channel. In order to have an idea of the size of the inelastic hole \(a_\mathrm{inel}\), we show in Fig. 9 the inelastic \(\pi \pi \) total cross sections in different isospin channels, defined as

$$\begin{aligned} \sigma ^\mathrm{inel}_{I}(s)= \frac{4 \pi }{k^2}\sum _{J} (2 J+1) \left[ 1- \eta _{IJ}^2(s) \right] \, . \end{aligned}$$
(97)

If we take \(\sigma ^\mathrm{inel}_I = 4 \pi a_\mathrm{inel}^2\), the largest inelastic cross section in Fig. 9 is compatible with an inelastic hole of less than half a fm, i.e. \(a_\mathrm{inel} \lesssim 0.5\) fm. Of course, this corresponds according to (97) to the contribution of all partial waves.

Fig. 9
figure 9

Inelastic cross-section in the different isospin channels as a function of the CM energy below \(\sqrt{s}=1.42\) GeV

As we have already seen, this loss of probability at the partial-wave level is parameterized in terms of a momentum-dependent inelasticity \(\eta _I^J(s)\), see (4). The imaginary part of the potential plays exactly the same role. At this point, we will assume that this complexification of the potential can be implemented just in the most inner layer of the potential. This assumption can be justified by analyzing the phenomenological inelasticities as a function of the impact parameter given by the relation \(b = (l+1/2)/p\), with l the angular momentum quantum number. In order to provide the range of the inelasticity, we show this dependence in Fig. 10. For the S waves it is found that only for the smallest value, \(b_\mathrm{min} = 1/(2 p_\mathrm{max}) \lesssim \Delta r \), the inelasticity is \(\eta \ll 1\). As we see, if we take \(\Delta r \sim 0.3 \) we may parameterize the inelasticity by one single energy dependent and complex parameter.Footnote 18

Fig. 10
figure 10

Inelastic profile as a function of the impact parameter \(b= (l+1/2)/p \) for the S0 (full, red), P (dotted, blue) and S2 (dashed, black) waves when the maximum CM energy becomes \(\sqrt{s}=1.42\) GeV

Hence, starting from our elastic previous description, we will assume \(\lambda _0\) in (62) to be a complex unknown function of the momentum and we will fit again the pseudo data given in [22] for the phase shift and inelasticity in [22] up to the maximum energy value provided \(\sqrt{s_\mathrm{max}}=1.42\) GeV. This procedure allows us to describe each partial wave (phase shift and inelasticities) exactly at each energy point. The results for the four \(r_c\)-values considered are plotted in Fig. 11, whereas the value of the real and imaginary part of the inner delta-shell layer is depicted in Fig. 12.

Fig. 11
figure 11

S0-, P- and S2-wave phase shifts (left panels) and inelasticities (right panels). The uncertainties are those quoted in [22], whereas solid-black, green-dashed, red-dotted and blue dot-dashed lines stand for the central results for \(r_c=0.9,\) 1.2, 1.5 and 1.8 fm, respectively. Nevertheless, the procedure described in the main text allows one to describe the input exactly at each energy point, so the four lines coincide exactly

Fig. 12
figure 12

Real and imaginary part of the inner delta-shell potential for the S0, P and S2 partial waves whereas solid-black, green-dashed, red-dotted and blue dot-dashed lines correspond to the results with \(r_c=0.9\), 1.2, 1.5 and 1.8 fm, respectively. The error bands have been computed using at each energy point a bootstrap with a uniformly distributed sample of 1000 points and taking the 68% of their distribution as the standard deviation

9.2 Traces of analyticity in the inelastic case

An interesting feature which can be appreciated in Fig 12 is the close resemblance of the real and imaginary parts of the inner delta-shell coefficient with the expected behavior from dispersion relations around a pole or a inelastic threshold. These two effects reflect in the S0 and P channels higher resonances or inelastic channels not explicitly included in the present optical potential analysis. Interestingly, Cornwall and Ruderman [7] found the fixed-r dispersion relation for the optical potential V(rs) given in (20). The implementation of such a dispersion relation in our analysis would reduce the number of fitting parameters in the inelastic region but it would also require a clear understanding of the high energy behavior. We leave this interesting investigation for future research.

The results for the real and imaginary part of the energy-dependent inner delta-shell when the chiral tail is included differ from those without, plotted in Fig. 12, only at low energies. The comparison is depicted in Fig. 13 showing again the rather small effect introduced by chiral corrections in the inelasticity parameter.

Fig. 13
figure 13

Real part of the inner delta-shell potential for the S0, P and S2 partial waves when the chiral tail is included or excluded. The only differ at low energies. Gray (light-gray), green (light-green), red (orange) and blue (cyan) bands correspond to the results with (without) the chiral tail for \(r_c=0.9\), 1.2, 1.5 and 1.8 fm, respectively

Another interesting possibility which deserves some further investigation is the generalization to the coupled channel case, to account explicitly for the opening of the \(K \bar{K}\) and \(\eta \eta \) thresholds, while keeping multi-pion channels in the inelasticity factor.

10 Subtractions, low energy constants and the number of parameters

In our construction of the chiral \(\pi \pi \) potential the low energy constants (LECs) have been discarded as they do not contribute for \(r \ne 0 \). In fact, the spectral representation does provide the same coordinate dependent potential regardless on the number of subtractions.

This raises the problem of where is this counterterm information gone within the present approach. In this last section we want to address the relation among subtraction constants in dispersion relations, low energy constants and the number of independent parameters within our coarse graining approach. We warn the reader that we have not succeeded in finding a unambiguous relation, which may ultimately be traced to two aspects. Firstly, it has to do with known ambiguities in mapping different regularization methods, namely the one used in \(\chi \)PT and the coordinate space regularization used here. Secondly, there is a difficulty in separating the short distance parameters in perturbation theory, particularly if we use a non-perturbative resummation scheme to fit the parameters.

In order to elaborate on this issue and appreciate the difficulties, we proceed in perturbation theory and re-write the problem in terms of a Fredholm integral equation,

$$\begin{aligned} u_{k,l}(r) = {\hat{j}}_{l} (kr) + \int \limits _0^ \infty \text {d}r'\, G_{k,l}(r,r') U_{l}(r')\, u_{k,l}(r') \, , \end{aligned}$$
(98)

where \({\hat{j}}_{l}(x)=x j_{l}(x)\) is a reduced spherical Bessel function of the first kind and \(G_{k,l}(r,r')\) is the Green’s function satisfying

$$\begin{aligned} \left[ \frac{\partial ^2}{\partial r^2} - \left( \frac{l(l+1)}{r^2} - k^2 \right) \right] G_{k,l}(r,r') = \delta (r-r') \, . \end{aligned}$$
(99)

An analytic expression for the Green’s function \(G_{k,l}(r,r')\) can be written in terms of two function ansatz u(r) and v(r) as

$$\begin{aligned} G_{k, l}(r,r') = u(r) v(r') \theta (r-r') + u(r') v(r) \theta (r'-r) \, . \nonumber \\ \end{aligned}$$
(100)

Inserting this in (99), it follows that u(r) and v(r) are solutions of the homogeneous equation with a unit Wronskian, i.e.,

$$\begin{aligned} \left[ \frac{\partial ^2}{\partial r^2} -\left( \frac{l(l+1)}{r^2} - k^2 \right) \right] u(r)= & {} 0, \nonumber \\ \left[ \frac{\partial ^2}{\partial r^2} - \left( \frac{l(l+1)}{r^2} - k^2 \right) \right] v(r)= & {} 0, \nonumber \\ u'(r) v(r) - u(r) v'(r)= & {} 1 \, . \end{aligned}$$
(101)

We choose one of the two solutions to be proportional to the regular one, \({\hat{j}}_{l} (kr)\). Then the other linearly independent solution with the desired Wronskian has to be proportional to \({\hat{y}}_{l}(kr)=kr\; y_{l}(kr)\), i.e. the reduced spherical Bessel function of the second kind. Therefore the Green’s function of the ordinary differential equation with the proper normalization can be written as

$$\begin{aligned} G_{k,l}(r,r') = \frac{1}{k}\,{\hat{j}}_{l} (kr_<)\, {\hat{y}}_{l} (kr_>) \, , \end{aligned}$$
(102)

where \(r_<= \mathrm{min} \{ r,r'\}\) and \(r_>= \mathrm{max} \{ r,r'\}\).

For the normalization of (98) the scattering amplitude can be written as

$$\begin{aligned} t_{IJ} (s) = -\frac{\sqrt{s}}{p^2} \int \limits _0^{\infty } \, \text {dr}\, \hat{j}_J (pr) U^I (r) u_{k,l}(r), \end{aligned}$$
(103)

which using the perturbative expansion inferred from reiteration of the integral equation gives

$$\begin{aligned} t_{IJ}^{(2)} (s)= & {} -\frac{\sqrt{s}}{p^2} \int \limits _0^{\infty } \, \text {dr}\, [\hat{j}_J (pr)]^2 U^{(2)} (r), \end{aligned}$$
(104)
$$\begin{aligned} t_{IJ}^{(4)} (s)= & {} -\frac{\sqrt{s}}{p^2} \int \limits _0^{\infty } \, \text {dr}\, \text {dr}' \, \hat{j}_J (pr)\times \nonumber \\&\quad U^{(2)} (r) G_J (r,r') U^{(2)} (r') \hat{j}_J (pr') \nonumber \\&\quad -\,\frac{\sqrt{s}}{p^2} \int \limits _0^{\infty } \, \text {dr}\, [\hat{j}_J (pr)]^2 U^{(4)} (r). \end{aligned}$$
(105)

Actually, the singularity structure in coordinate space suggests introducing a short distance cut-off \(r_c\) and hence a short distance potential \(U_\mathrm{Short}(r)\).Footnote 19 That means that in order to identify numerically the counterterms in perturbation theory we may split the integrals as

$$\begin{aligned} \int \limits _0^\infty = \int \limits _0^{r_c} + \int \limits _{r_c}^\infty \end{aligned}$$
(106)

so that we get

$$\begin{aligned} \int \limits _0^{r_c} r^2 \, \text {dr}\, [j_l (pr)]^2 U_\mathrm{Short} (r) + \int \limits _{r_c}^\infty r^2 \, \text {dr}\, [j_l (pr)]^2 U_\mathrm{Long} (r)\, \end{aligned}$$

and apparently the matching to the one-loop result could be undertaken in a straightforward fashion. This is actually not so, since this implies disentangling the fitting parameters in a chiral expansion, namely

$$\begin{aligned} \lambda _n \equiv U(r_n) \Delta r = \lambda _n^{(2)} + \lambda _n^{(4)} + \cdots \end{aligned}$$
(107)

where we can arbitrarily shift \(\lambda _n^{(2)} \) and \( \lambda _n^{(4)} \) by equal but opposite constants keeping \(\lambda _n\) constant, say the values of Tables 1 or 6.

This situation is not exclusive to the present approach and in fact is common to all unitarization schemes. For instance in the IAM method the fitted LEC’s are different than those of \(\chi \)PT or the predicted unitarized amplitudes from \(\chi \)PT develop huge uncertainties [68].

The fact that this perturbative matching can ultimately provide a successful description and still fulfilling the condition \(\lambda _n^{(2)} \gg \lambda _n^{(4)} \) is expected. We note, however, that the small changes between the \(\lambda _n\) parameters corresponding to the case without TPE potential listed in 1 and the case with TPE listed in 6 suggest this hypothesis.

11 Conclusions

The optical potential in \(\pi \pi \) scattering is a meaningful object under the most common and general assumption of the validity of the Mandelstam double spectral representation. Therefore, it plays a relevant role in the analysis of such an interaction within a invariant mass formulation of the relativistic two-body problem. Contrary to the more employed Bethe–Salpeter equation, such an approach is free from well-documented spurious singularities, which are triggered by incomplete calculations embodying subsets of Feynman diagrams with particle exchange. In addition, it is a much simpler approach in the CM frame, as it effectively reduces to a Schrödinger equation for equal mass particles.

Within such a framework, in the present paper we have analyzed \(\pi \pi \) scattering from a coarse-grained point of view at distances smaller than the elementarity radius of the pion. This means sampling the interaction in coordinate space at a resolution of the order of the shortest de Broglie wavelength. In our case, where we choose a maximal CM energy of \(s_\mathrm{max}=2\) GeV\(^2\), the resolution turns out to be \(\Delta r \sim 0.3\) fm. As a result, we obtain successful fits in the S0, P and S2 partial waves with the expected number of parameters. We have also analyzed the role of inelasticities by an energy dependent coarse grained interaction. The implications from chiral symmetry have also been analyzed in terms of a long-distance potential featuring the two-pion-exchange mechanism. This potential has been determined for the first time.

A non-perturbative renormalization of the amplitude, based on boundary conditions in coordinate space, is precluded by the energy dependence of the chiral potential. For a finite short distance cut-off about 1.2–1.5 fm this energy dependence becomes irrelevant.

A somewhat surprising result of our analysis is that explicit chiral corrections play a minor role, since at the distances above 1.2–1.5 fm chiral potentials are almost negligible. A rewarding consequence in this regard concerns the extraction of phase-shifts from energy shifts calculated in a finite box in the relative distance by means of the Lüscher formula [47]; its applicability requires the interaction to vanish above a given size which provides a lower limit to the size of the box, L. While the mere \({\mathscr {O}} (e^{-m_\pi L})\) nominal estimate suggests \(L \gtrsim 2 \mathrm{fm}\), our analysis is compatible with taking \(L \sim 1-1.5 \mathrm{fm}\), a much smaller value.

We have also shown that the quantum mechanical re-interpretation of the problem does not spoil the proper analytical properties of the partial wave scattering amplitude, which are explicitly fulfilled by the Feynman diagrams. Thus, the conventional dispersion relations with both a left-hand cut due to particle exchange in the crossed channel and the right hand cut due to unitarity are fulfilled. The standard N/D decomposition of the amplitude is explicitly realized, albeit without subtractions. We expect subtractions to be encoded in the short distance components of the potential, but the explicit determination and identification of the subtraction constants in terms of the short-distance parameters remains to be accomplished.

An important advantage of the coarse graining perspective is that the number of independent parameters is determined a priori by the shortest wavelength and the available crossing constraints. These constraints become increasingly large as we increase the angular momentum of the partial wave. This point deserves further investigation and the determination of all partial waves with this minimal number of parameters is left for future research.

A traditional objection to the successful data-driven unitarity methods is based on their lack of a power counting scheme, reflecting field dependence and off-shell ambiguities absent in the conventional bona fide EFT framework. The present coarse graining approach to \(\pi \pi \) scattering makes no further assumptions than those usually made. Namely, it implements unitarity in the elastic regime and it matches \(\chi \)PT in perturbation theory above a given separation distance, which is estimated to be about 1.2–1.5 fm. In contrast, it does have the advantage that we can estimate the number of fitting parameters a priori. More demanding fits due to an increase in the number of (consistent) data should not require more fitting parameters, but rather determining short distance parameters more accurately. Besides, the method is quite simple as it parameterizes the unknown short distance behavior regardless of any power counting. This also allows a discussion a posteriori of the chiral contributions which turn out to be minor.