Unitarization Technics in Hadron Physics with Historical Remarks

We review a series of unitarization techniques that have been used during the last decades, many of them in connection with the advent and development of current algebra and later of Chiral Perturbation Theory. Several methods are discussed like the generalized effective-range expansion, K-matrix approach, Inverse Amplitude Method, Pad\'e approximants and the N/D method. More details are given for the latter though. We also consider how to implement them in order to correct by final-state interactions. In connection with this some other methods are also introduced like the expansion of the inverse of the form factor, the Omn\'es solution, generalization to coupled channels and the Khuri-Treiman formalism, among others.


Introduction
The use of perturbative calculations within Chiral Perturbation Theory (ChPT) as input for non-perturbative S-matrix based methods is a general procedure several decades old. Due to the fact that ChPT results are perturbative, given in terms of an expansion organized in increasing powers of the external four-momenta and light pseudoscalar masses, unitarity is only satisfied in the perturbative sense, similarly as in a standard Born series. A well-known example in this regard is the calculations in Quantum Electrodynamics, where the expansion is done in powers of α (the fine structure constant), so that if the leading calculation is O(α) then unitarity contributions start at O(α 2 ), and they have the form (leading amplitude) 2 ×(phase space)= O(α 2 ). However, the fulfillment of unitarity implies to square the calculated amplitudes, and not only their calculation up to a lower order.
It is somewhat astonishing that already in 1970 one can read about motivations for unitarizing phenomenological chiral Lagrangians, introduced to construct realizations of the current algebra approach [1]. Rephrasing the original remarks by Schnitzer [2,3], the ideas he put forward are still the main reasons to advocate the unitarization of ChPT amplitudes: 1. The tree approximation to the scattering amplitudes violate badly unitarity. This could also be said for perturbative unitarity, at least in some partial waves. 2. The Lagrangians are nonlinear and nonrenormalizable, which makes difficult to compute higher-order corrections. Nowadays, we would better say that there is a rapid proliferation of counterterms as the order of the calculation increases in ChPT, with the status of the art at the two-loop level in ChPT. It is typically simpler and much more predictive to implement lower-order calculations of ChPT within non-perturbative methods. 3. Even if such corrections could be computed, the resultant renormalized perturbation series would probably diverge, since the perturbation parameter has the strength characteristic of strong interactions. This is clear from phenomenology because hadronic interactions are characterized by plenty of resonances and a rapid saturation of unitarity in many partial-wave amplitudes (PWAs).
Although the interest in the present writing is on ChPT and the associated chiral expansion, among those early papers of Schnitzer we also quote Ref. [3], which builds an approach to extract the consequences of chiral symmetry based on constructing a particular realization of the current-algebra, such that their associated Ward-Takahashi identities are satisfied, and two-body unitarity is implemented by means of an effective-range-type parameterization (a method discussed in Sec. 3.1).
One possibility to improve the agreement with data of the perturbative calculations within ChPT is to apply the chiral series expansion to an interaction element of the amplitude, which is afterwards implemented within non-perturbative techniques. This is one of the basic ideas behind unitarization methods for the chiral series of scattering amplitudes. The first works along these lines considered the application of an effective-range-type parameterization to unitarize ππ scattering [3,4], once the ππ scattering lengths were calculated at leading order in ChPT by the application of the current-algebra techniques and the partial conservation of the axial-vector current (PCAC) [5]. A similar unitarization method was later applied to the first calculation at next-to-leading order (NLO) in the chiral counting of the ππ PWAs in the chiral limit (m π = 0). The calculation of the latter ones, as well as their unitarization by applying a generalized effective-range expansion (ERE) [6] were undertaken in Ref. [7], as discussed in more detail in Sec. 3. This calculation explicitly shows that the ππ scattering in the chiral limit is finite.
Special mention deserve the pioneering works by Truong and collaborators [8][9][10], in which the role of the isoscalar S-wave ππ final-state interactions (FSI) are stressed, having significant effects on several physical processes. In a first instance [8], the authors correct the current-algebra result for η → 3π by S-wave ππ rescattering, a reaction which is also discussed by the application of the Khuri-Treiman (KT) [11] formalism in Sec. 4.3. For that Ref. [8] multiplies the current-algebra transition amplitude by an Omnès function [12,13], in which the isoscalar scalar ππ phase shifts were, however, taken from experiment. As a result, the Watson final-state theorem is fulfilled [14]. In another work of the saga [9], the input phase shifts were generated consistently by the theoretical scheme followed after taking the one-loop ChPT result for the scalar and vector pion form factors and imposing the fulfillment of unitarity, as discussed in Sec. 4.2. It is also stressed that in this form a resummation of the ChPT series is achieved that may also give rise to resonant effects.
The 90's of the past century experienced a boost in the interest of applying chiral effective field theories (EFTs) for the study of nuclear interactions. To large extent this was triggered by the seminal articles of Weinberg [15], in which the systematic application of ChPT order by order to calculate the nuclear potentials V is established. As the chiral order increases however extra derivatives with respect to r act on the potential, so that it becomes more singular for r → 0. Because of this complication the application of ChPT for the calculation of the low-energy NN PWAs by implementing the chiral potentials in quantum-scattering integral equations has not reached a complete satisfactory status yet. These ideas also triggered the application of ChPT to the study of the non-perturbativeKN scattering in coupled channels, particularly in connection with the Λ(1405) [16][17][18][19].
The non-perturbative character of the NN interactions is due to two basic aspects. (i) One of them is a quantum effect of kinematical origin within the typical scales of the problem. The typical distance of propagation of two nucleons as virtual particles is l NN ∼ 1/E kinetic ∼ m/p 2 ∼ (m/m π )b π , where m and m π are the nucleon and pion masses, respectively. The range of the NN interactions is given by the Compton wavelength of the pion b π = m −1 π (in our unitsh = c = 1). As m/m π 1 this travel distance for virtual particles is large enough for having several repetitive collisions between the propagating two nucleons. The same conclusion is reached if one focuses on the propagation of real nucleons. For a typical three-momentum m π they have a velocity of order m π /m. Thus, the time for crossing a distance b π is ∼ m/m 2 π >> b π . (ii) Nonetheless, if the coupling between two nucleons were small enough the scattering would be perturbative despite (i). This does not happen since the coupling due to one-pion exchange between two nucleons is of the order g 2 A m 3 π / f 2 π , where f π 93 MeV is the pion weak decay constant and g A 1.26 is the axial coupling of the nucleon. This factor times l NN implies the dimensionless number l NN 16π which is about 0.5. Therefore, the NN interactions should be treated non-perturbatively as a general rule. In this equation the phase-space factor 1/16π is included, which accounts for the two-nucleon propagation in all directions. We also distinguish in Eq. (1.1) the scale [20,21] which has a striking small size despite it is not proportional to m π . This is another reflection of the the non-perturbative character of the NN interactions.
The unnaturally large size of the NN S-wave scattering lengths (a s ), so that they are much bigger in absolute value than the Compton wavelength of the pion, |a s | b π , introduces a new scale at low energies. For instance, the scattering length for the isovector 1 S 0 NN PWA is a s −25 fm. As a result, the dimensionless number in Eq. (1.1) becomes even larger by a factor |a s |m π 1. Therefore, when the center of mass (CM) three-momentum is smaller than m π , in which case the ERE applies, the NN interactions are manifestly non-perturbative and the NN potential has to be iterated. Precisely, in this energy region one finds the bound state of the Deuteron in the PWA 3 S 1 and an antibound state for the 1 S 0 .
It is also remarkable the confirmation by unitarization methods of the existence of the σ resonance in pion-pion interactions at low energies. This resonance is nowadays called f 0 (500) in the PDG [22] and its pole position is given there at (400 − 550) − i(200 − 350) MeV. The standard view of ChPT, based on the spontaneous symmetry breaking of SU(2) × SU(2) chiral symmetry [23], considered as highly unlikely that such a low-mass resonance could happen in ππ scattering, where the small expansion parameter is claimed to be m 2 π /m 2 ρ 0.03, with m ρ the mass of the ρ(770) meson. However, for the isoscalar scalar ππ scattering the unitarity corrections are affected by a large numerical factor that could make that the actual expansion parameter in the momentum squared dependence of this PWAs is actually much larger. This was explicitly shown in Ref. [24] by performing the exercise of determining the value of the renormalization scale µ needed in order to generate the ρ(770) pole by unitarizing the leading-order (LO) ChPT amplitude. It was obtained that a huge unnatural value for the realm of QCD was needed, with µ 1 TeV, while the same value for generating the σ resonance had the natural value in QCD of µ 1 GeV.
It is instructive to also show the main equations for the completion of this exercise. For the isovector vector ππ interactions, where the ρ(770) resonance appears, one has the unitarized expression of the LO ChPT PWA T 11 (s), which reads [24] (as also discussed in Sec. 5) Here, T I J is the PWA of the two-pion system with isospin I, J is the angular momentum and the LO ChPT amplitude is (s − 4m 2 π )/6 f 2 π . The function g(s) corresponds to the two-pion unitarity loop function, given by In turn, the unitarized expression for the I = J = 0 PWA, which contains the f 0 (500) pole, is [24,25] T 00 (s) = f 2 π s − m 2 π /2 where the LO ChPT PWA is (s − m 2 π /2)/ f 2 π . The main difference between Eqs. (1.3) and (1.5) is the factor 6 dividing the LO ChPT T 11 (s) compared to T 00 (s), because s m 2 π in the region where the σ or ρ poles lie. Indeed, in order to get a resonance of mass m ρ in the I = J = 0 PWA one needs a µ of around 1.8 GeV, in comparison with µ around 1 TeV that is needed in the I = J = 1 case. The reason for this dramatic change in the needed values of µ is because g(s) only depends logarithmically on µ. This facts reflects that the unitarity corrections for the scalar isoscalar sector are numerically enhanced, and this enhancement is enough to generate resonant effects that strongly impact the phenomenology in the physical region and make fallacious to think in the possibility to reach accuracy by a straightforward application of ChPT for many reactions. As a result, the infinite set of unitarity bubble diagrams should be resummed in order to account for this numerical enhancement.
This phenomenon is also seen in the strong corrections affecting the I = J = 0 ππ scattering length originally calculated by Weinberg at LO [5] with current algebra methods. The expressions for the a I J scattering lengths up to NLO or O(p 4 ) in ChPT from Ref. [23] are affected by chiral loops which are the dominant NLO contribution in the limit M π → 0. They read: It follows then that a 00 has the largest NLO contribution in the limit m π → 0. In order to appreciate better the relatively large size of this correction, it is worth comparing it with the pion mass calculated in ChPT up to NLO [23], withm 2 π the bare mass squared. The NLO term here is a factor 9 smaller in absolute value than that for the a 00 . Indeed, this was one of the reason for developing a non-perturbative dispersive approach that could provide an improvement in the prediction of the ππ scattering lengths. The idea is to make use of the Roy equations [26] and to match with ChPT in the subthreshold region, where the ChPT expansion is better behaved, since it is away from the threshold cusps [27,28]. In this way, the two subtractions constants needed for solving the Roy equations can be predicted by ChPT, applied at different orders. The resulting convergence properties of the prediction for the ππ scattering lengths is much improved and a reliable estimate of the uncertainties can be also provided. Another more recent advance was the development of a new set of Roy-like equations in Refs. [29,30], the so-called GKPY equations. The difference is that these equations have only one subtraction instead of two and, e.g., they have given rise to an accurate determination of the f 0 (500) pole from experimental data, without relaying on the ChPT expansion.
In this work we have always followed the order of first discussing scattering, mostly in PWAs, and then FSI. After a brief review on unitarity in Sec. 2 we then move to discussing several unitarization methods in Sec. 3, establishing links between them. The generalized ERE, the K-matrix approach, the Inverse-Amplitude Method and the Padé resummation are then considered. The Sec. 4 is dedicated to the implementation of re-scattering effects in probes and several methods are presented, with some of them clearly related to the already presented ones in the section dedicated to scattering. Other ones are introduced that could be applied to any given set of PWAs. We dedicate Sec. 5 to discuss the N/D method for PWAs and FSI. This section ends with a brief account of the exact N/D method recently developed for non-relativistic scattering. The last section is a brief ending discussion.

Unitarity
The S-matrix operator S, gathering the transition matrix elements between in and out states, is unitary: The scattering operator T, also called the T matrix, is introduced such that the S matrix in terms of it reads The unitarity of the S matrix implies in turn that T fulfills the following relation, Multiplying this relation to the left by T −1 and to the right by T † −1 we have the interesting equation If we also take into account that the T matrix is symmetric in partial waves, the previous relation implies that the imaginary part of the matrix of partial-wave amplitudes (PWAs) is fixed by unitarity. If we write this matrix in brief as T L , then the right-hand side (rhs) of the previous equation, in the energy region in which it is saturated by two-boy intermediate states, can be written as In this equation, q is the diagonal matrix of the center of mass (CM) three-momentum for every two-body intermediate state and θ is also another diagonal matrix whose matrix elements are 1 for √ s larger than the threshold and 0 otherwise. Eq. (2.5) is equivalent to the probably more familiar unitarity equation for PWAs The phase-space diagonal matrix q(s)θ/(8π √ s) is sometimes denoted for short by ρ(s). Now, if we are considering simultaneously stronger and weaker interactions (where the weaker interactions could correspond e.g. to actually electromagnetic or weak probes, while the stronger ones typically refer to the strong interactions among hadrons), at leading order in the weaker interaction (they are supposed to be proportional to some small dimensionless parameter), the unitarity relation reads, In this equation F represents the matrix elements of the T matrix involving the weaker interactions, so that they vanish if these interactions are neglected altogether, while still the stronger ones would be acting. In the latter equation we have taken that the weaker interactions act in the initial state, otherwise write iF † T on the rhs of Eq. (2.7). In PWAs the unitarity relation of Eq. (2.7) gets its simplest form. In the physical region for the reactions to occur it reads where ρ j (s) corresponds to the phase space of the intermediate hadronic states (integrations could also be involved for multiparticle states) and s is the standard Mandelstam variable corresponding to the total center-of-mass (CM) energy squared. The opening of the threshold for the channel j, s th,j , is accounted for by a Heaviside function of the form θ(s − s th,j ) included as part of ρ j (s). For the one-channel case the sum on the right-hand side of the Eq. (2.8) collapses to just one term, where T (s) is the corresponding uncoupled PWA. Since the left-hand side of the equation is real then it follows that the phase of the form factor F(s) and the phase shift of T (s) are the same modulo π. This is the well-known Watson final-state theorem.

ERE, K-matrix, IAM and Padé approximants
Along this section we follow a transversal discussion relating different unitarization approaches, like the (generalized) effective-range expansion (ERE), K-matrix parameterizations, the Inverse Amplitude Method (IAM) and the Padé approximants.

ERE and K-matrix approaches
In the early days of PCAC, soft pions theorems and realizations based on chiral Lagrangians, it was customary to refer as (generalized) ERE to a unitarization method based on the identification of a remnant in the inverse of a PWA free of right-hand cut (RHC) which was expanded in powers of p 2 . The standard ERE was originally derived in Ref. [31] for NN interactions which, for an uncoupled PWA, has the form The remnant part is identified with p 2 +1 cot δ as it is well known, because of the relation between the T and S matrices in the normalization used typically for the ERE, which is the one used in Eq. (3.1). Namely, the steps are The NN scattering is non-relativistic (NR), with m 2 >> p 2 at low energies, so that the expansion of p 2 +1 cot δ is a Taylor series in p 2 . However, for pion-pion interactions, where p 2 ∼ m 2 π in the region of interest both theoretical and experimentally speaking, the series expansion in p 2 is a Laurent series for the S waves. The reason is because the Adler zeroes required by chiral symmetry in the S-wave PWAs [32], despite there is no centrifugal barrier. The latter is present for the higher partial waves, ≥ 1, which implies the standard zero at threshold so that T vanishes as p 2 for p → 0.
The phase space factor for relativistic systems changes in comparison with the NR expression of Eq. (3.1). The steps are the same as in Eqs. (3.1) and (3.2), but now instead of T one should use T / √ s so that In more recent times, the remnant of T −1 after discounting the factor −ip/ √ s, required by unitarity, cf. Eq. (2.5), is called the inverse of the K-matrix, K , instead of p cot δ . In this notation, T is written as Of course, Eqs. (3.2), (3.3) and (3.4) can be generalized straightforwardly to a matrix notation for coupled-channel scattering, with T and K replaced by the matrices T L and K L , respectively. The inverse of the later is usually referred as the M L matrix, M L = K −1 L [33]. We are surprised that in these first works, e.g. [2][3][4]6,7,34,35], it was common to refer to the (generalized) ERE without any mention at all to the K-matrix approach, a notion much more common in later times and, in particular, for more recent papers based on the unitarization of Chiral Perturbation Theory (ChPT). Probably this is related to the fact that the K-matrix parameterizations have been used in many instances in the literature over large energy intervals in order to fit experimental data. As a result, it does not really make sense to keep any memory of a particular threshold, as it is the case for the ERE. Indeed, in those earlier papers referred the basic object of study was ππ scattering or the π vector form factor.
Another fact worth stressing is that in those earlier references the expressions finally used for T −1 L had better analytical properties than the ones typically found later in papers using the K-matrix approach. The reason is because the later ones keep only the term −ip/ √ s in T −1 L while, in the first papers referred, the non-trivial analytical function h(s), which is 8πg(s), cf. Eq. (1.4), modulo a constant, was used by performing a dispersion relation (DR) along the RHC. Namely, The function g(s) is an analytical function of s in the cut complex s plane, having the RHC along the real s axis for s > 4m 2 π . As a trivial byproduct, the zero at s = 0 that occurs in the phase space factor −ip/ √ s in the simplest K-matrix parameterizations is absent when using the function g(s), which is the correct analytical extrapolation of the two-body unitarity requirement above threshold. Indeed, the removal of this spurious singularity at s = 0 was the argument used in Ref. [4] to construct the function h(s) without using any DR. This reference also notices the presence of the Adler zeros in the I = 0, 2 S-wave ππ PWAs and similar expressions to Eq. (1.5) are proposed for these PWAs. The main different, an important one indeed, between Eq. (1.5) and Ref. [4] is that the function g(s), contrary to h(s), contains a subtraction constant which is absent in the function h(s) of Brown and Gobble [4]. This is a crucial fact for the right reproduction of important features in low-energy ππ scattering, like the generation of the f 0 (500) resonance pole in good agreement with the latest and more sophisticated determinations [22]. As a matter of fact, the predicted I = J = 0 ππ phase shifts in Ref. [4] are around a factor 2 smaller than data for the energies in between 500 − 700 MeV, while the I = 2 S-wave ππ phase shifts are too large in modulus by the same factor. These deficiencies in the approach of Ref. [4] are cured once the subtraction constant of Eq. (3.6), with a natural value for µ 1 GeV, is taken into account [25].
For the I = J = 1 ππ PWA Ref. [4] performs a generalized ERE up to and including the effective range, The parameter a 1 is fixed from the current algebra prediction [5], a 1 = 1/12π f 2 π , while r 1 is fixed by the vanishing of the real part of T 11 (s) −1 at s = m 2 ρ . The resulting equation is therefore, Let us notice that r 1 /2 in Eq. (3.7) can be also considered as a subtraction constant of g(s). Attending to Eq. (3.5) the relation is with δr a correction of around a 20% of the term explicitly shown. This simple calculation illustrates the discussion at the Introduction regarding the huge unnatural value µ 1.7 TeV that results by the matching in Eq. (3.9), while the expected value is around 1 GeV. As a result of this analysis the authors of Ref. [4] predicted the width of the ρ(770) to be 130 MeV and the I = J = 1 phase shifts up to 1000 MeV, in good agreement with later experimental determinations. They also give an expression for the coupling of the ρ → ππ (g ρππ ) in terms of f π and m ρ , which drives to the KSFR relation [36], f 2 ρ = m 2 ρ /2 f 2 π , if one assumes vector-meson dominance (VMD) [37,38]. Here f ρ is the coupling of the ρ-photon transition which is equal to g ρππ within VMD [38]. The authors summarize their research by stating that the fulfillment of the low-energy current-algebra constraints together with the inclusion of extra energy dependence as required by general principles, such as it follows by implementing two-body unitarity and the correct analytical properties of PWAs, are able to provide good results in a large energy range, which is much larger than the one naively expected for current-algebra results. This is a conclusion that has been strengthened along the years, at the same time that the chiral calculations have been improved going to higher orders and the unitarization methods have become more sophisticated.

ERE and IAM
Already at 1972 the calculation of the NLO ChPT amplitude was worked by Lehmann [7] in the chiral limit (m π → 0), much earlier than the seminal paper by Gasser and Leutwyler [23], which established the general framework for ChPT at O(p 4 ). The author did not need to work out the chiral Lagrangians at NLO order because he only used unitarity, crossing symmetry and analyticity to work out the chiral loops. The point is that because of unitarity a PWA satisfies Eq. (2.6). However, unitarity is only satisfied perturbatively in the chiral expansion, so that if we denote by T 4 (s) a one-loop ChPT PWA and T 2 (s) its LO, then perturbative unitarity requires that The PWA T 4 (s) has left-hand cut (LHC) and RHC. The discontinuity along the RHC is twice i T 4 (s), because of the Schwarz reflection principle. A DR that results by considering a closed circuit engulfing the RHC, implies the following contribution to T 4 (s), Three subtractions have been taken because T 2 (s) at most diverges like s in the limit s → ∞. By invoking crossing one can build up the one-loop contributions from the t-and u-channels for a given process. As usual the Mandelstam variables are indicated by s, t and u (s + t + u = 0 for massless pions).
In Cartesian coordinates for the pions and treating all of them on equal footing, so that all of them are e.g. incoming, one can write for the scattering amplitude Here crossing has also been used to properly exchange the arguments of the A(s, t, u) function. The previous expression is manifestly symmetric in the indices i 3 and i 4 which also implies that, because the pions are bosons, A(s, t, u) is symmetric under the exchange t ↔ u. Since the isospin coordinates run only from 1 to 3, two out of the four pions have the same coordinates necessarily.
In the calculation of Ref. [7] the resulting expression for A(s, t, u) has two parts. One of them corresponds to DR integrals of the type in Eq. (3.11), in all the s-, t-and u-channels, which can be evaluated in an algebraic close form. The other contribution is a second-order polynomial in the Mandelstam variables, whose general expression can be written as a + bs + cs 2 + c (t 2 + u 2 ), which can also be extra constrained. E.g. a = 0 because Goldstone particles do not interact in the limit in which masses and four-momenta vanish. The term bs is order p 2 and it is already accounted for in T 2 (s). As a result the one-loop calculation of Lehmann only involves two unknown parameters, nowadays typically called counterterms because they are associated to bare parameters appearing at the NLO chiral Lagrangian.
In terms of the A(s, t, u) amplitude one can calculate the different ππ isospin PWAs [39], T I J . An interesting point of Ref. [7] is the perturbative matching in the chiral expansion of the calculated PWAs at O(p 4 ) with the ERE expression for a PWA, cf. Eq. (3.2). The subtle point is that the former only satisfies unitarity in a perturbative way, as discussed above. Therefore, writing in the massless case that is not right. The correct procedure is to write a chiral expansion of 1/T I J up to NLO and from there to identify cot δ I J , Taking into account the perturbative unitarity satisfied by t 4 , one can extract from here the NLO expression for cot δ I J (with a numerical normalization factor properly chosen) as This is indeed the first example that we know of a paper in the literature deriving the expression of a PWA as This formula, generalized to any other two-body PWA and also to coupled channels, is the basic one for the so-called IAM [39,40]. It also illustrates the connection between these earlier treatments based on the ERE and this more modern method, which was named IAM after the general framework for the one-loop calculations in ChPT was established in Ref. [23]. The approach of Ref. [7] has the advantage over the previous ERE of Refs. [2][3][4]34,35] that chiral one-loop contributions in the crossed channels are also kept, so that the LHC is reproduced up to NLO in the inverse of the PWA.
The extension of Eq. (3.16) up to two-loop ChPT can be done straightforwardly by expanding the inverse of (t 2 + t 4 + t 6 ) −1 up to next-to-next-to-leading order (NNLO), or O(p 2 ). The result is, (3.17) Taking into account that perturbative unitarity requires that t 6 = 2t 2 ρ t 4 , it follows that T I J given by Eq. (3.17) fulfills exact unitarity, T I J = −ρ. The Eq. (3.17) is the IAM at the two-loop order [41].

IAM and Padé approximants
Another non-perturbative method used with the aim of improving the convergence of the Quantum Field Theory (QFT) calculations in perturbation theory is the Padé resummation technique [42]. It is also a unitarization method that was applied since the early days of current algebra calculations by Refs. [43,44], in which the linear σ model was considered too. An interesting qualitative agreement with data for the ππ S-, Pand D-waves was reported, despite the limitations of the theoretical input.
Given a function f (z) that is analytic at z = 0, its Taylor series expansion around this point converges within the circle of radius R, which is the distance to the nearest singularity. However, it is also known that the value of f (z) at a point z 1 within its domain of analyticity, but beyond the radius of convergence of the Taylor series around z = 0, is fixed by the coefficients in the later expansion. The idea of the Padé method is to provide a resummation of the Taylor series and build an approximation of f (z) beyond the radius of convergence of its Taylor series around z = 0.
The Padé approximant [n, m] is given by the ratio of two polynomial functions P n (z) and Q m (z) of degrees n and m, respectively, which has the same n + m first derivatives as f (z) at z = 0. Namely, Notice that in particular the approximant [n, 0] is identical up to O(z n ) with the Taylor series of f (z) at z = 0. It is also typically the case that the Padé approximations usually provide an acceleration in the rate of convergence of the Taylor series itself. For instance, one can write that By iteration it can be expressed as a continued fraction, which are particular cases of Padé approximants, etc. Let us compare this first four Padé approximants with the first four terms in the Taylor series, The formulas for the IAM at one-and two-loop ChPT, Eqs. (3.16) and (3.17), respectively, can also be obtained as Padé approximants, where a generic small parameter ε accounts for the chiral order. Formally, we then write t 2 → ε 2 t 2 , t 4 → ε 4 t 4 and t 6 → ε 6 t 6 . The one-loop IAM is a [1, 1] Padé approximant: To solve this type of equation, typically found in Padé approximants, it is convenient to rewrite Eq. (3.21) as By matching the different powers of ε 2 one has that From which it follows that The result of the matching is the same as in Eq. Therefore, as Eq. (3.17).

Final-State Interactions
As a canonical example of taking into account the final-state interactions (FSI) that correct the production processes due to weaker probes because of the rescattering by the stronger interactions, we start with the unitarization of the vector pion form factor, F V (s), within the ERE approach of Ref. [34]. We next move to the Omnès solution for a form factor and also consider the scalar pion form factor, F S (s), paying attention to a caveat in the use of an Omnès function that one should properly consider. Along the discussion we introduce the way FSI are treated in Ref. [9], as it is probably the first paper in which NLO ChPT is unitarized to account for FSI following the basic notions of unitarity, Watson final-state theorem and use of an Omnès function, which are the basic elements employed in the different modern approaches to resum FSI [12,45]. We end this section with a basic account of the Khuri-Treiman (KT) approach for η → 3π decays.

ERE, the Omnès solution and coupled channels
The application of the ERE for implementing the FSI of the pion vector form factor was pioneered in Ref. [34]. The main aim in that paper concerns the corrections because of the finite width of the ρ to the VMD dominance relation between Γ(ρ → e + e − ) and Γ(ρ → π + π − ), as well as to characterize the energy shape of Γ(e + e − → π + π − ).
Ref. [34] implemented the relationship between the I = J = 1 ππ PWA and the pion form factor F V (s) by writing F V (s) = T 1 (s)/t 2 (s), with t 2 (s) the LO ChPT amplitude. This relation is a consequence of the Omnès representation in the approximation in which: i) One assumes that the only zero in T 1 (s) in the region of interest is the one at threshold, s = 4m 2 π , because of the = 1 centrifugal barrier; ii) one also assumes the dominance of the ρ(770) exchange so that it is a good approximation to consider that T 1 (s) is dominated by s-channel dynamics. Under these assumptions T 1 (s) is approximately given by the Omnès function on the rhs of Eq. (4.24) times (s − 4m 2 π )/48π f 2 π . Thus, guaranteeing that F V (0) = 1 because of conservation of total charge. Next, Ref. [34] performs the same ERE of Ref. [4], which we have already discussed, cf. Eq. (3.7), which allows to finally write the form factor in a successful manner as The authors of Ref. [34] simplify further this expression by removing those terms involving the expansion of the real part of h(s) around s = m 2 ρ that are at least of O(s − m 2 ρ ). They finally write Again, one concludes that the extrapolation of the current-algebra results plus the extra energy dependence that arises by implementing the basic principles of two-body unitarity and analyticity allows one to reach much higher energies than expected, even above the 1 GeV frontier. The fact of writing a form factor proportional to a given PWA is usually employed in many cases in the literature. The basic reason is to write down an expression for the coupled form factors F i (s) that automatically satisfies the constrained imposed by the two-body unitarity, cf. Eq. (2.8). Following Ref. [33] one then writes where the sum is over the strongly-coupled channels. The functions α i are real and they are also expected to be smooth because all the RHC features in F i (s) are included in the PWAs T ij (s). As a result, the α i should not have nearby singularities, if any. They could involve crossed-channel cuts which could be mimicked typically by parameterizing these functions by low-degree polynomials. Nonetheless, in the case of the low-energy interactions of the lightest pseudoscalars, like pions, an extra feature is the presence of the Adler zeroes in the S waves. In particular, for I = J = 0 we have already discussed that this Adler zero is around s A = m 2 π /2, cf. Eq. (1.5). This Adler zero is a characteristic feature of the interaction of the Goldstone bosons, as said, but not necessarily for their production through external currents. To handle with such cases Ref. [33] proposes to explicitly remove the Adler zeroes in the T ij (s), when they are present, and any necessary zero in the production process is then explicitly included in the prefactors. Denoting by T ij (s) = T ij (s)/(s − s A ij ), with s A ij the Adler zero in T ij (s), the final expression proposed is For the case of only one coupled channel, the form factor can be expressed in terms of an Omnès function Ω(s). Due to the Watson final-state theorem the continuous phase of the form factor ϕ(s) is the same as the phase shift δ(s) for the PWA T(s). The Omnès function results by performing a DR for the logarithm of the function f (s) = F(s)Q(s)/P(s), where P(s) and Q(s) are the polynomials whose only roots are the possible zeros and poles of F(s), respectively, which are assumed to be finite in number. The discontinuity of log f (s) along the RHC is the discontinuity of its imaginary part, and it is given by 2iϕ(s). We can then write for the where we have taken n subtractions assuming that ϕ(s) does not diverge stronger than s n−1 when s → ∞. The Omnès function Ω(s) is defined in terms of ω(s) as One can always normalize the Omnès function such that Ω(0) = 1, which fixes a 0 = 1. In this manner we always take at least one subtraction. It is also clear that the ratio is a meromorphic function of s in the first RS of the cut complex s plane, being analytic in the whole plane if F(s) has no bound states. As it is well known, any analytical function in the whole complex plane is either a constant or it is unbounded, which is then the case for R(s) too under the stated assumptions. Therefore, diverges as much as or stronger than Ω(s) for s → ∞. The function ω(s) would have severe divergences for s → ∞ if its DR required for convergence more than one subtraction. The reason is that if ϕ(s)/s n−1 (n ≥ 2) has no zero limit for s → ∞, the DR for ω(s) would be affected by logarithmic divergences like s n−1 log s which could not be cancelled by the subtractive polynomial. In such circumstances it would be required that R(s) is a non-trivial analytical function in order to cancel such divergences and guarantee that F(s) can be represented as a DR. If the conditions are met for a DR of log F(s)Q(s)/P(s), cf. Eq. (4.6), then R(s) = Q(s)/P(s) is a rational function. Thus, from the previous analysis, we conclude that the DR of ω(s) in Eq. (4.6) involves only 1 For a more extensive discussion on the Muskhelishvili-Omnès problem the reader can consult Refs. [12,45]. one subtraction and it is then necessary that |ϕ(s)/s| < s −γ for some γ > 0 in the limit s → ∞. We can then write the following representation for F(s), ) The presence of P(s) makes clear that one can fix de normalization of the Omnès function, Ω(0) = 1, without any loss of generality. The asymptotic behavior of Ω(s) in the limit s → ∞ can be calculated as follows. Let us rewrite ω(s) in Eq. (4.12) as being the limit s → ∞ dominated by the logarithmic divergence, as the other two terms in this equation are constants. It follows from here the limit behavior This result, together with Eq. (4.10), implies that the asymptotic behavior for F(s) is where C Ω and C F are constants, and p and q are the number of zeros and poles of F(s), respectively (or equivalently, the degrees of P(s) and Q(s), in this order). Two interesting consequences follow from Eq. (4.16): i) If the asymptotic high-energy behavior of F(s) is known to be proportional to s ν , then ii) Under changes of the parameters when modeling strong interactions one should keep Eq. (4.17) unchanged.
As ν is a known constant, then It is worth stressing that by using Eq. (4.10) one can guarantee that Eq. (4.18) is fulfilled, while this is not the case for Ω(s). The use of this function without taking proper care of the rational function P(s)/Q(s), included in the expression for F(s) in Eq. (4.10), could drive to an unstable behavior under changes of the parameters, e.g. in a fit to data. This problem was originally discussed in Ref. [46] in connection with the scalar form factor of the pion F S (s), to which we refer for further details in the discussion that follows. This form factor is associated with the light-quark scalar source,ūu +dd, and is defined as where u and d are the up and down quarks,m is their masses, and s = (p + p ) 2 . Because of the quantum numbers of the non-strange scalar source, the FSI occur in the isoscalar scalar meson-meson scattering, introduced in Sec. 3. There, we discuss the Adler zero required by chiral symmetry and the pole of the σ or f 0 (500) resonance, being both of them related by unitarity, analyticity and chiral symmetry. At around the two-kaon threshold, √ s = 991.4 MeV, the KK channel makes a big impact. This energy almost coincides with the sharp emergence of the f 0 (980) resonance, which gives rise to a rapid increase of the ππ isoscalar scalar phase shifts δ 00 (s), since it is a relatively narrow resonance [22], cf. Fig. 2. The elasticity parameter η 00 also experiences a sharp reduction as soon as the KK channel open, since the f 0 (980) couples much more to KK than to ππ [47]. This phenomenon causes an active conversion of the pionic flux into the kaonic one.
The rapid rise of the isoscalar scalar ππ phase shifts, also implies the corresponding rise of the phase of the isoscalar scalar PWA T(s), ϕ(s), because they coincide below the KK threshold, i.e. for √ s < 2m K . However, above this energy the rise of ϕ(s) is interrupted abruptly if δ 00 (s K ) < π, with s K = 4m 2 K , while in the opposite case ϕ(s) keeps increasing. Quite interestingly, the two situations can be connected by tiny variations in the values of the parameters in the hadronic model, while keeping compatibility with the experimental phase shifts at around the f 0 (980) mass.
As a result, there is a jump in the limiting value of Ω(s) because ϕ(∞) changes by π. Thus, in order to keep constant Eq. (4.18) under an increase by π in ϕ(∞) for δ 00 (s K ) > π, it is necessary to increase p by one unit, so that a zero is necessary in F S (s) that is not present when δ 00 (s K ) < π. For completeness, we also mention that had we required the continuity from δ 00 (s K ) > π to δ 00 (s K ) < π then an extra pole (in the first RS) should be added. This latter scenario can be disregarded in ππ scattering because of the absence of bound states. 2 Let s 1 be the value of s at which the pion scalar form factor has a zero for δ(s K ) > π. Then, we can write an Omnès representation of the pion scalar form factor in terms of a modified Omnès function such that F S (s) = F S (0)Ω(s). From here it is clear that s 1 can be fixed by the requirement that F(s 1 ) = 0.
Because of the Watson final-state theorem in the elastic region we can write that F(s) = |F(s)| sin δ 00 (s)|/ρ(s) and it vanishes when δ 00 (s 1 ) = π, which allows to determine s 1 from the knowledge of δ 00 (s). The context clarifies whether the same symbol Ω(s) actually refers to Eq. (4.7) or Eq. (4.20). A clear lesson from the discussion here is that possible troubles could occur when using an Omnès function in fitting the free parameters because an unstable behavior could arise due to a jump in ϕ(∞). These regions of dramatic differences in exp ω(s) are separated by a discontinuity of ϕ(s) in the parametric space. As a consequence, it is important in the fitting process to satisfy the condition Eq. With respect to the difference between ϕ(s) and δ 00 (s), as indicated above, the f 0 (980) dominates the behavior of the isoscalar scalar meson-meson scattering around 1 GeV, and couples much more strongly to kaons than to pions. For instance, Ref. [47] calculates that its coupling to kaons is a factor 3 larger than that to pions. This makes that the mixing between the pion and kaon scalar form factors is suppressed, following each of them its own eigenchannel of the I = J = 0 PWAs.

The IAM for FSI
The first step of Ref. [9] is to write down twice subtracted DRs expressions for the scalar and vector pion form factors, F S (s) and F V (s), respectively, as Here, δ 00 (s) and δ 11 (s) are the J = 0 and 1 isoscalar and isovector ππ phase shifts, in this order. These DRs can be interpreted as singular integral equations (IEs) for the form factors F S (s) and F V (s) [50]. Let us remark, as in Ref. [9], that the solutions of the IEs of Eqs. (4.21) and (4.22) for F S (s) and F V (s), respectively, can be expressed in terms of the associated Omnès functions [51]. In the approximation of identifying the phases of the form factors with the phase shifts, strictly valid only for the elastic region, one has the approximate expressions where P S (s) and P V (s) are polynomials that take into account the zeros (if any) of the form factors in the first or physical RS. At the one-loop order in ChPT, or equivalently next-to-leading order NLO or O(p 4 ), we can replace inside the dispersive integrals of Eq. (4.21) the ππ scattering PWAs at leading order, f 0 (s) = sin δ 00 e iδ 00 = δ 00 (s) + O(p 4 ) = σ(s) 16π The function h(s) is defined in Eq. (3.5). By proceeding in an analogous way, a similar expression holds for the vector form factor at this level of accuracy, O(p 4 ), There is an important difference between the scalar and vector form factors. The unitarity corrections are enhanced by around a factor 6 for the former compared to the latter, because the leading order ChPT amplitude is around a factor 6 larger, compared Eqs. (4.25) and (4.26), as first noticed in Ref. [24] and already discussed above in detail. By invoking the Watson final-state theorem, one can calculate from the perturbative expressions of F S (s) and F V (s) in Eqs. (4.27) and (4.28) the ππ phase shifts for J = 0 and 1, respectively. Nonetheless, since the form factors are calculated perturbatively one should proceed consistently in order to extract from this perturbative information the corresponding phase shifts. In this way, denoting by F 2 (s) the LO form factors and by F 4 (s) = F r 4 (s) + iF i 4 (s) their NLO contributions, with the superscripts indicating the real (r) and imaginary (i) parts, we then have for the Watson final-state theorem: from where the phase φ can be extracted. We show in Fig. 1 by the solid lines labelled (a) the resulting phase shifts (with s given in units of m 2 π ) as compared with the experimental data available at the time of Ref. [9]. The left panel is for F S (s) and the right one for F V (s). It is clear that there is a strong departure between the calculated phase shifts from the NLO ChPT form factors and the experimental values even at low values of s. This is also clearly true for the moduli squared of F V (s), indicated by the dashed line with the label (c) as compared with the experimental data points (empty circles). It is particularly visible there the emergence of the resonance ρ(770), which dominates the phase shifts and |F V (s)| 2 , with tails extending up to threshold and affecting the low-energy results. This phenomenon can only be captured approximately in SU(2) ChPT by the large size of the counterterm¯ 6 , as estimated in Ref. [23]. For the vector case the cause of the large higher-order contributions is clearly associated with the prominent role played by the ρ(770) resonance. In turn, for the scalar sector the enhanced RHC is the one blamed for such effects. Indeed, these strong contributions from unitarity and analyticity even drive to the emergence of a pole in the complex s plane, the σ or f 0 (500) resonance, as already discussed, cf. Eq. (1.5).
Ref. [9] discusses that the application of the chiral series expansion should be performed on the inverse of the form factor rather than on the form factor itself. The main reason lies on sound and general grounds, as provided by unitarity and analyticity. Let us consider a DR representation of F −1 (s), analogous to Eq. (4.21). The point to be stressed is that the imaginary part of F −1 (s) is expected to be much smoother than the imaginary part of F(s) itself in the elastic region. The reason is that the imaginary part of the inverse of the form factor satisfies, because of unitarity in PWA, that As F(s) and T(s) share the same resonances, their propagators cancel in the ratio T(s)/F(s) that gives F −1 (s). Then, this ratio is expected to be smoother than F(s) = ρ(s)F(s) * T(s), where this cancellation does not occur but rather the resonance effects in F(s) and T(s) mutually enhance each other because of the product involved.
Then, let us write down a twice-subtracted DR for the inverses of the form factors F S (s) and F V (s). First, we neglect by now the possible presence of zeroes in the form factors in the 1st Riemann sheet (RS), which give rise to poles in the inverse of the form factors. We will discuss later the issue of a zero in F S (s) for certain types of T matrices, as discussed and first noticed in Ref. [46]. As a result we write, . . (4.36) One can see in Fig. 1 that the solid lines labelled by (b), calculated from Eqs. (4.35) and (4.36), are much closer to the experimental points than the phases deduced from the perturbative NLO ChPT calculation of the form factors in Eqs. (4.27) and (4.28), which are the solid lines indicated by (a). The same dramatic improvement also happens for the modulus squared of F V (s) calculated from Eq. (4.36) and shown by the dashed line with the label (d), as compared with the data points given by the empty circles. Notice that this improvement is achieved by employing the same perturbative input, namely the NLO ChPT results. It is a matter of properly reshuffling the chiral expansion in a way clearly motivated by unitarity and analyticity.

KT formalism
The KT formalism was originally developed by Ref. [11] to study the K → 3π decays and, up to including two-body intermediate states, it allows to implement unitarity and crossing symmetry. Later on this approach has been applied to study extensively the η → 3π decays, among others. These decays violate isospin because the G parity of the η is +1 and that of the pion in −1, so that it is proportional to m u − m d in pure QCD.
The application of ChPT to the η → 3π decays has been controversial, until accepting without any complex that FSI are so strong that a non-perturbative unitarization method is needed to be implemented in order to be able to confront well with experimental data [8]. The earliest calculations using current-algebra techniques obtained a value for the η → π + π − π 0 of around 65 eV [52], too small compared with the experimental result Γ(η → π + π + π 0 ) = (300 ± 12) eV [22]. Roiesnel and Truong [8] stressed that a non-perturbative calculation taking care of the isoscalar-scalar ππ FSI, by employing an Omnès function on top of the current-algebra result, increases the decay width up to 200 eV. A few years later, the NLO ChPT calculation [53] gives (160 ± 50) eV, which implies a large correction by a factor 2.4 over the LO calculation in the right direction, but still too small by around a factor of 2. In addition, the parameter α, typically employed in the parameterization of the Dalitz plot for the decay η → 3π 0 , is positive at NLO ChPT [53] while experimentally it is negative, α = −0.0318 ± 0.0015 [22]. The calculation at NNLO in ChPT of the η → 3π decays was performed in Ref. [54] but the proliferation of new counterterms prevented a sharp result. If resonance saturation is assumed to estimate the NNLO ChPT counterterms then the Dalitz plot parameters are not well reproduced. One then concludes that the η → 3π decays are sensitive to the detailed values of the O(p 6 ) counterterms, so that an accurate calculation requires a precise knowledge of their values. This controversial situation stimulated the interest in developing sophisticated calculations combining ChPT and non-perturbative methods, within unitarized ChPT [8,55,56] and the KT formalism [57][58][59][60][61].
We now describe the basic points of the one-channel KT formalism for η → 3π decays and refer the reader to Refs. [61,62] and the recent review [45] for further details. In particular, the generalization to coupled channels was worked out in Ref. [61], given in a more compact matrix notation in Ref. [45].
Let us consider the decay η(p 0 ) → π + (p 1 )π − (p 2 )π 0 (p 3 ), which is related by crossing symmetry to the scattering reactions η(p 0 )π 0 (−p 3 ) → π + (p 1 )π − (p 2 ) in the s-channel, η(p 0 )π − (−p 1 ) → π 0 (p 3 )π − (p 2 ) in the t-channel, and η(p 0 )π + (−p 2 ) → π + (p 2 )π 0 (p 3 ) in the u-channel. The Mandelstam variables s, t and u are given by The crossing-symmetry relations are These amplitudes in turn can be decomposed in scattering amplitudes with well defined isospin, M I (s, t, u), as The inversion of these relations gives us the M I (s, t, u), In the KT formalism the S and P waves are the ones that are subject to a non-perturbative treatment.
The PWAs have a RHC above the two-pion threshold s > 4m 2 π . Instead of writing the unitarity constraint as in Eq. (2.6), one should consider it as giving the discontinuity along the RHC because of the two on-shell intermediate pions. We then write 3 which is the relation finally used. A crucial feature of the KT formalism is to write down A(s, t, u) as the sum of three functions of only one Mandelstam variable, M 0 (s), M 1 (s) and M 2 (t) [60,61] A(s, t, u) which is invariant under the exchange t ↔ u, a feature that can be seen as a consequence of charge-conjugate invariance. This representation is valid up to O(p 8 ) in ChPT [54,60] because then the ππ D waves also contribute and higher polynomials in (s − t) and (s − u) would be required. The derivation of Eq. (4.44) can be understood by considering only J ≤ 1 PWAs in the s-channel and taking into account the isospin decomposition for the process ηπ 0 → π + π − and the crossed-channel ones, cf. Eq. (4.39). In this way, for the s-channel process there is no I = 1 contribution, which only happens in the crossed ones, cf. Eq. (4.39). As this is a P-wave we then write it as M 1 (t)(s − u) + M 1 (u)(s − t), that also keeps explicitly the symmetry under the exchange t ↔ u. The I = 0 contribution can only happen in the s-channel, because for the other channels the third component of isospin is not zero. This is the M 0 (s) contribution in Eq. (4.44). Finally, regarding the I = 2 it is clear from Eq. (4.39) that it appears in the combination −2M 2 (s)/3 + M 2 (t) + M 2 (u). Taking the expression for A(s, t, u) in the ones of M I (s, t, u), as given in Eq. (4.40), it follows that Due to the fact that in decay channel all the three pions are on-shell in the region (m η − m π ) 2 ≥ s ≥ 4m 2 π this is another source of imaginary part from the crossed-channel cuts that are also on-shell. For s = (m 2 η − m 2 π )/2 the branch point singularity at t, u = 4m 2 π happens for cos θ = ∓1. These crossed-channel cuts can be separated from the RHC one by giving a vanishing positive imaginary part to m 2 η and then proceed by analytical continuation in m 2 η [63].
Writing down the PWAs for I J = 00, 20 and 11 we have We have also introduced in Eq. (4.46) the angular averages The functionM I (s) have not discontinuity across the RHC so that the discontinuities of the PWAs M (I J) (s) can be expressed as, (20) (s) sin δ (20) Following the same steps as above in Eq.  where P (m) Requiring that A(s, t, u) diverges linearly at most at infinity in the Mandelstam variables [58], then M 1 (s) should be bounded by a constant and M 0 (s), M 2 (s) would diverge linearly at most in the limit s → ∞. Furthermore, we also know the asymptotic behavior in the same limit for the Omnès functions, cf. Eq. (4.15), with |Ω (I J) (s)| → s −δ (I J) (∞)/π . Depending on δ (I J) (∞) the value of m should be adjusted to the required asymptotic behavior of M I (s). For instance, Ref. [61] assumes that δ (00) (∞) = π, δ (11) (∞) = π and δ (20) = 0, so that m = 2 for I = 0 and m = 1 for I = 1, and 2.
The DRs in Eq. (4.54) constitute a set of coupled linear IEs because the angular averages M I n are also expressed in terms of the M I (s) functions. A standard way for solving these equations is by iteration. The subtraction constants can be determined by matching with the NLO ChPT calculation of A(s, t, u) and/or fitted to data, as done in Refs. [58,61]. A clear improvement is obtained in the calculated decay width for the η → π + π − π 0 in Ref. [58], where the value Γ η→π + π = 283 ± 28 eV was obtained. Other improvements concern the parameter α for characterizing the amplitude for η → 3π 0 in its Dalitz plot. NLO ChPT gives a value α = 0.0142 while the KT treatment of Ref. [61] gives α = −0.0337 (12), to be compared with the PDG average value of α = −0.0318(15).

The N/D method
In this section we elaborate on different aspects of the N/D method, first introduced in Ref. [6] to study uncoupled ππ PWAs. We first review on this method, discuss in more detail the limit in which the crossed-channel dynamics is neglected [24], and afterwards elaborate on how the latter can be treated perturbatively within the N/D method [18,64]. These results can also be used to take into account FSI in production processes [65,66]. For the case of NR scattering, thanks to recent developments [67], it is possible to know the exact discontinuity of a PWA along the LHC for a given potential. In this way, one can generate the same solutions as in the Lippmann-Schwinger (LS) equation, together with other ones that cannot be obtained in a LS equation when mimicking the short-distance interactions by contact terms in the potential [68].

Scattering
For the scattering of particles with equal masses there is only a LHC for s < s Left because of crossing. However, when the particles involved have different masses there are also other types of cuts in the complex s plane due to crossing. For instance, for the scattering of particles a + b → a + b, in addition to a LHC there is also a circular cut for |s| = m 2 b − m 2 a [69], where for definiteness we have considered that m b > m a . Nonetheless, when we refer in the following to the LHC we actually mean all the crossed-channel cuts. Indeed, had we taken instead the complex p 2 plane all the cuts would be linear and only a LHC would be present [69].
We introduce the N/D method following Ref. [24]. The uncoupled case is discussed first and afterwards we move to coupled-channel scattering. The discussion is restricted to two-body intermediate states. The discontinuity of the inverse of a PWA T (s) along the RH is 2i times its imaginary part, being the latter fixed by phase space because of unitarity, cf. Eq. (2.5).
In the N/D method T (s) is expressed as the quotient of two functions, where N (s) stands for the numerator function and D (s) for the denominator one. The former has only LHC and the later RHC.
To enforce the right kinematical threshold behavior of a PWA, vanishing as p 2 , Ref. [24] divides T (s) by p 2 , The N/D method is then applied to this function, It follows then that the discontinuities of N (s) and D (s) along the LHC and RHC, respectively, are Consistently with Eq. (5.7), the DR for N (s) can be written as . , and the last integral can indeed be performed algebraically. This is a linear IE for D (s) with s along the LHC. Once this solved one can calculate D (s) for s ∈ C and, in particular, along the physical region, s + i with s ∈ RHC. Other types of IEs could be deduced by taking more subtractions independently in D (s) and N (s). Fore more details on this respect the reader can consult [70]. The expression in Eq. (5.9) can be shortened and simplified for equal mass scattering with mass m by taking s 0 = 4m 2 , because then p(s ) 2 = (s − 4m 2 )/4. It follows that, The last integral in the previous expression can be written in terms of g(s), Eq. (1.4).
One of the subtraction constants can be fixed because we can freely choose the normalization of D (s) and N (s), since their ratio and analytical properties are invariant under a change in normalization. The standard choice is to fix D (0) = 1. However, given ∆ (s) along the LHC, the solution is not unique because of addition of extra subtraction constants in D (s) and N (s).
Historically, the possible addition of Castillejo-Dalitz-Dyson (CDD) poles [71] was the clear indication that extra solutions could be obtained even if ∆ (s) is assumed to be known along the LHC. They give rise to zeros of T (s) along the RHC and each zero comprises two real parameters, its residue and position. Phenomenologically the CDD poles correspond to the short-distance dynamics underneath the scattering process and might also be related to the addition of bare states [72]. Let us notice that T (s) −1 does not exist at a zero of T (s) and, therefore, Eq. (2.5) is not defined there. As in Ref. [71] let us introduce the auxiliary function λ(s) such that Denoting by s i the zeros of T (s) along the real axis above threshold, we can write λ(s) from Eq. (5.12) as where the λ(s i ) are a priori unknown. Thus, Eqs. (5.11) and (5.13) allow us to write

The last term in the previous equation can be rewritten as
The contribution a m s m and Eq. (5.14) can be rewritten as where a m , γ i and s i are constants not fixed by the knowledge of ∆ (s), and the CDD poles give rise to the last term.
Interesting results can be deduced under the approximation of neglecting the LHC, ∆ (s) → 0. Eq. (5.8) then becomes 17) and N (s) is just a polynomial, which can be reabsorbed in D (s) by dividing simultaneously both functions by N (s) itself. The expression for T (s) then becomes The number of real free parameters in the previous equation is + 1 + 2M , with M the number of CDD poles. A priori there is nothing to prevent the generalization of Eq. (5.18) such that some s i could also lie below threshold. We could adjust the position and residue of a CDD pole such that the real part of D (s) vanishes at the desired position. This would give rise to typical resonance behavior above threshold, or to a bound-state pole if this happens below threshold. This is why the parameters of the CDD poles are typically associated with the coupling constants and masses of the poles in the S matrix. In other instances, the CDD poles are needed because the presence of a zero cannot be related to ∆ (s), but they respond to fundamental constraints in the theory. This is the case of the Adler zeroes in QCD [32], which already occur at LO in the chiral expansion, while ∆ (s) = 0 only at NLO. It is therefore necessary to account for them by including CDD poles, such that the derivative of the PWA at the zero corresponds to the inverse of the residue of the CDD pole, γ i . For the ππ Adler zeroes the latter could be fixed in good approximation by the LO ChPT result. The other + 1 parameters emerge by having enforced the correct behavior of a PWA near threshold, which should vanish as p 2 . Let us stress that Eq. (5.18) gives the general form of an elastic PWA when the LHC contributions are neglected. Phenomenologically this assumption could be suited if the LHC is far away and/or if it is suppressed for some reason [73]. The free parameters in Eq. (5.18) can be fixed by fitting experimental data and/or by reproducing the Lattice QCD (LQCD) results at finite volume or when varying some of the QCD parameters, like N c or the quark masses [74][75][76][77]. Ref. [24] focuses on meson-meson scattering, whose basic theory is QCD. It studied the S-and P-wave two-body scattering between the lightest pseudoscalars (π, K and η), as well as the related spectroscopy. It was found that the full nonet of scalar resonances [78] f 0 (500), f 0 (980), a 0 (980) and K * 0 (800) arose from the self-interactions among the lightest pseudoscalars, while the more massive resonances f 0 (1370), f 0 (1500), a 0 (1450) and K * 0 (1430) stem from a nonet of bare resonances with a mass around 1.4 GeV. In addition, Ref. [24] included a bare scalar singlet with a mass around 1 GeV which gives also a contribution to the f 0 (980) [79]. Later on, Ref. [80] extended this model by including more channels and could determine a glueball state affecting mainly the f 0 (1700) with a reflection (because of the ηη threshold) on the f 0 (1500) as well. Of course, the same Eq. (5.18) can be applied to other interactions, e.g. Ref. [64] studied WW scattering in the Electroweak Symmetry Breaking Sector.
The generalization of Eq. (5.18) to coupled channels is rather straightforward by employing a matrix notation, where the T matrix in coupled channels is a matrix denoted by T L (s). As in Eq. (5.18) we take from the onset that crossed-channel dynamics can be neglected in a first approximation. Thus, the matrix element T L,ij (s) is proportional to p i i p j j , which gives rise for odd orbital angular momentum (unless i = j) to another cut between s th;i and s th;j due to the square roots in the expressions of p i and p j as a function of s. To avoid this cut we define the matrix T L , analogously to Eq. (5.2), as In this equation, the symbol p L corresponds to a diagonal matrix with matrix elements 20) and m 1i and m 2i are the masses of the two particles in the same channel i. The matrix unitarity relation along the RHC then reads where ρ(s) is another diagonal matrix whose elements are ρ i (s). The next step proceeds with the generalization to coupled channel of Eq. (5.1) by writing T L as 22) with N L and D L two matrices, the former only involves LHC and the later RHC, respectively. In our present case without LHC, the matrix elements of N L are polynomials functions. Multiplying to the left N L and D L in Eq. (5.22) by N L −1 we can always make that N L = I and write, with R(s) a matrix of rational functions which poles produce the CDD poles in D L . Let us notice that all the zeros in detT L correspond to CDD poles in the detD L (s). This is the generalization of the CDD poles for the coupled-channel case.
The resulting expression for T L (s) in Eq. (5.23) can be also recast as with g(s) the diagonal matrix with matrix elements g i (s) defined as where a i (s 0 ) is a subtraction constant and s 0 the subtraction point. The result of this integration can also be written as The parameter µ is a renormalization scale, such that a change in the value of µ can always be reabsorbed in a corresponding variation of a i (µ), while the combination a i (µ) − 2 log µ is independent of µ. The unitarity loop function g i (s) corresponds to the one-loop two-point function where ω j = m 2 ji + p 2 and the total four-momentum p 1 + p 2 is indicated by P. The integral in Eq. (5.27) diverges logarithmically, which is the reason why a subtraction has been taken in Eq. (5.25). The Eq. (5.26) also results by employing dimensional regularization and reabsorbing the diverging term in a i (µ).
Let us elaborate on the so-called natural value for the subtraction constant. The function g i (s) given by Eq. (5.26) has the value at threshold, This expression is compared with the one that results by evaluating g i (s) in terms of a three-momentum cutoff Λ. The resulting expression for the function g i (s), and denoted by g Λi (s), can be found in Ref. [39]. The natural size of a three-momentum cutoff in hadron physics is the inverse of the typical size of a compact hadron, which is generated by the strong dynamics binding quarks and gluons. Thus, according to this estimate we take Λ 1 GeV. For NR scattering g i (s) and g Λi (s) (m 1i , m 2i |p|) are given by the value at threshold of every function plus −ip/(8π(m 1 + m 2 )) + O(p 2 ) [70]. The value at threshold of g Λi (s th ) can be worked out explicitly with the result [74] g Λi (s th ) = − 1 By equating Eqs. (5.28) and (5.29) the following matching value for a i (µ) results, One should employ µ Λ 1 GeV in Eq. (5.30) to estimate the natural value for the subtraction constants, a procedure originally established in Ref. [18]. In this way, both the renormalization scale µ and the cut off Λ are used with values suitable to the transition from the low-energy effective field theory (EFT) to the shorter-range QCD degrees of freedom. As an example, let us take ππ scattering and Λ = 1 GeV. Then, from Eq. (5.30) The Eq. (5.24) is adequate for including perturbatively the LHC contributions in the T-matrix T L . This can be achieved by matching order by order with a calculation within an EFT. For instance, this has been used many times taking as input one-loop calculations in ChPT [18,47,64,74,75,[81][82][83][84][85]. The procedure is as follows. Let us take a meson-meson scattering amplitude calculated in ChPT up to one-loop or O(p 4 ), T L = T 2 + T 4 + O(p 6 ). Then the chiral expansion of Eq. (5.24), with V = V 2 + V 4 + O(p 6 ), g = O(p 0 ) [64], reads at LO, and at NLO, and similarly for higher orders. Thus, up to NLO the matching equations fix V 2 and V 4 to The LHC contributions arise because crossed-channel loops are calculated order by order in the ChPT results for T L . Figure 2. Results from Ref. [25] with only one free parameter for the S-wave meson-meson scattering with I = 0 and 1. From top to bottom and left to right, the isoscalar scalar ππ → ππ and KK → ππ phase shifts, the ππ inelastic cross-section with the same quantum numbers and a π 0 η event distribution around the isovector scalar a 0 (980) resonance are plotted. For more details and references of the experimental papers we refer to Ref. [25].
In order to appreciate the powerful of the method for some reactions we consider the LO matching, that is, with V = V 2 , applied in Ref. [25] to study the meson-meson S-waves with I = 0 and 1. This is a coupled-channel study with ππ and KK for I = 0 and πη and KK for I = 1. It is certainly remarkable that only one free parameter entered in the successful calculation of the PWAs from the ππ threshold up to around 1.2 GeV. This is shown in Fig. 2 by the ππ → ππ, KK → ππ phase shifts, inelastic reaction and π 0 η event distribution around the a 0 (980), from top to bottom and left to right, respectively. The resonances f 0 (500), f 0 (980) and a 0 (980), clearly visible in Fig. 2, are generated dynamically from the interactions between the pseudoscalars. The free parameter is the three-momentum cut-off with natural size Λ 1 GeV used in the evaluation of the unitarity-loop functions g Λi (s) employed in this study.
It is also the case in some instances [24,47,74,75,80,81] that the ChPT expansion is complemented with the exchange of bare resonances fields, so that the tree-level amplitude is crossing symmetric. One typically improves the convergence properties of the chiral expansion by including bare resonance fields because of the (partial) saturation of the chiral counterterms by the resonance exchanges [86]. Then, the matching process is undertaken up to O(hp 4 ,h 2 ), which means to neglect any two-loop contribution and any one-loop contribution beyond O(p 4 ). In this way, one could consider one-loop contributions involving higher orders because of the explicit inclusion of the resonance fields. The matching proceeds as in Eqs. The perturbative solution of the N/D equations with respect to the LHC contributions can also be organized as an iterative solution on increasing number of insertions of ∆ . The first-iterated N/D method consists on taking only one power of ∆ in the integrand of the DRs for D (s) and N (s). The approximation is obtained by settling D (s) = 1 into the integrand for the DR of N (s), Eq. (5.8), which is then denoted as N ;1st . Then, .
Since ∆ (s) is known the DR integral could in principle be calculated. This is usually a tree-level amplitude that can also be calculated in QFT, from which indeed ∆ (s) is actually derived. Therefore, we assume that in the first iterated N/D method N ;1st is also given. As a result, the calculation of D (s) in this approximation, denoted by D :1st (s), just reduces to perform the integration The first-iterated N/D method was used in Ref. [88] to discuss ππ scattering within linear realizations of chiral symmetry, taking into account the exchanges of a σ and ρ resonances. More recently, it has been employed to study ρρ scattering in Ref. [89] by taking the pure gauge-boson part of the non-linear chiral Lagrangian with hidden-local symmetry [90,91]. Its generalization to the SU(3)-related vector-vector scattering was undertaken in Ref. [92]. These studies were motivated by the earlier ones in Refs. [93], with still an on-going productive discussion in interpreting the results.

FSI
Let us consider the unitarity relation for a form factor, Eq. (2.9), with the expression of the T matrix in PWAs T L as given in Eq. (5.24). It then results that Since g(s) = −ρ(s) it is clear that g(s) + 2iρ(s)θ(s) = g(s) * , so that from Eq. (5.37) we have that along the RHC it is fulfilled that This formalism has been employed by Ref. [95] to study the γγ →meson-meson fusion reactions. Refs. [66,96] used it to study the scalar form factor of the pion (and of other pseudoscalar mesons) in connection with J/ψ and D decays, and Ref. [65] analyzed the vector form factor of the pion. This formalism was also very important to unveil the two-pole structure of the Λ(1405) in Ref. [18], because in previous studies the πΣ event distribution for this resonance was always taken to be proportional to the modulus squared of the πΣ → πΣ I = 0 S-wave.

The exact N/D method in NR scattering
For non-relativistic scattering one can calculate for a given potential the exact discontinuity of a PWA along the LHC. This has been a recent advance in S-matrix theory achieved by Ref. [67], to which we refer the reader for further details. The key point was to extrapolate analytically the LS equation to complex three-momenta for off-shell scattering. The solution of the LS equation for half-off-shell scattering is an analytical function in the off-shell three-momentum complex q plane with vertical cuts which extend along the lines (±)p ± iλ with |λ| ≥ µ 0 . Here the ± symbols are unrelated, µ 0 is the lightest particle exchanged, and p is the on-shell three-momentum (fixed by the energy E of the process, E = p 2 /2µ, with µ the reduced mass). E.g. for NN scattering the lightest particle exchange is the pion and µ 0 = m π . We denote in the following a PWA for half-off-shell scattering as T (q, p), where q is the off-shell three-momentum and p the on-shell one.
The discontinuity we are interested in, e.g. for its later application to the N/D method, is After some mathematical derivations that can be consulted in Ref. [67], this discontinuity can be obtained by solving an ordinary linear IE. This IE is written in terms of the discontinuity of the potential in momentum space v (q, p). Its writing gets simplified by usingv defined bŷ The discontinuity of the potential entering into the IE is with − < < + and + → 0 at the end of the calculation. After this preamble, the sought IE is (p = ik, k ≥ µ 0 and n = 2 + 2) In terms of f (ν) the discontinuity ∆ (p 2 ) is given by Thus, we need to solve the IE for ν ∈ [−k + µ 0 , k − µ 0 ], and the range of the integration in the IE for f (ν) is finite for a given p, contrary to the LS equation. This IE can be solved without ambiguity because ∆v(ν, ν 1 ) can be determined for a given potential and with it f (ν) by solving Eq. (5.46). For a general potential it is convenient to employ its spectral decomposition, v(q, p) = ∞ µ 0 dμ 2 η(μ 2 ) (q − p) 2 +μ 2 + . . . (5.48) where η(μ 2 ) is the spectral function, and the ellipsis indicate possible subtractions that due to its polynomial nature do not give contribution to the discontinuity of the potential. In terms of the spectral decomposition we can write that ∆v (ν, ν 1 ) = − 2 π ∞ µ 0 dμ 2 η(μ 2 )ρ(ν 2 , ν 2 1 ;μ 2 )θ(ν 1 − ν −μ) . The function ρ(ν 2 , ν 2 ;μ 2 ) is a polynomial in its argument and its fixed by the partial-wave projection involved in the case of interest. For brevity in the presentation offered here we have just referred to the uncoupled case, but the formalism can also be generalized easily to evaluate the LHC discontinuity for coupled PWAs [67]. A potential is said to be singular if for r → 0 it diverges stronger than 1/r 2 or as α/r 2 for α + ( + 1) < −1/4. In the opposite case the potential is said to be regular [67]. In the ChPT calculation of nuclear potentials the increase in the order of the calculation implies typically an increase in the degree of divergence of the potential for r → 0, because off-shell momentum factors give rise to spatial derivatives. This fact is the main reason why the original Weinberg's program for solving nuclear properties once the chiral potentials are calculated order by order has not been taken to full completion.
The resulting ∆ (p 2 ) obtained by solving the master Eq. (5.46) has a different qualitative behavior depending on whether the potential is regular, attractive singular or repulsive singular. General arguments, based on the scaling properties of the function ρ(ν 2 , ν 2 1 ;μ 2 ), were given in Ref. [67] to explain such differences in the behavior of ∆ (p 2 ). Explicit examples were also worked out in Ref. [67] corresponding to actual PWAs in NN scattering, with the chiral potential calculated at different chiral orders, from LO up to NNLO. The function ρ(ν 2 , ν 2 1 ;μ 2 ) is a polynomial in ν and ν 1 of degree m. Then, the argument of Ref. [67] follows by considering a re-scaling by a parameter τ of the variables k, ν and ν 1 in the limit k µ 0 . It follows from Eq. (5.46) that the n th iterated solution for f (ν) is subject to a re-scaling by τ (n+1)m−(2 +1)n = τ (m−2 −1)n+m . (5.50) The point is whether m − 2 − 1 is smaller or larger than zero. In the former case we have the behavior corresponding to a regular potential, so that each extra iteration implies at least an extra factor of 1/k and for k → ∞ the discontinuity ∆ (−k 2 ) tends to its Born approximation. However, when m − 2 − 1 > 0 each iteration increases the power of k in the asymptotic behavior of ∆ (−k 2 ), becoming more and more divergent as n increases. This is the situation for a singular potential. For the regular potentials ∆ (p 2 ) tends to its Born term contribution which vanishes at least as 1/p 2 for p 2 = −k 2 and k 2 → ∞. For such type of ∆ (p 2 ) it was shown in Ref. [70] that any N/D IE, irrespectively of the number of subtractions taken, has solution. However, for singular potentials the resulting |∆ (p 2 )| grows faster than any polynomial in the same limit. This is clearly shown in Ref. [67] by log-log plots in which the slop of |∆(−k 2 )| continuously grows with increasing k 2 . As a dramatic consequence of this result is that it is not possible to write down a DR representation for a NR PWA if the potential is singular. However, it is still possible to use the N/D method because what matters for the N/D IEs is the product ∆ (−k 2 )D (−k 2 ). The denominator function is known to behave asymptotically as s −δ(∞)/π , cf. Eq. (4.15), and δ(∞) = Nπ, with N the number of bound states, because of the Levinson theorem. It turns out that the number of such stats is infinite for attractive singular potentials [97] and, in this case, D (−k 2 ) vanishes also faster than any power law.
The exact N/D method is defined in Ref. [67] as the N/D method but using ∆ (p 2 ) stemming from Eqs. (5.46) and (5.47), which is the exact LHC discontinuity of the full PWA for a given potential. In this way, we showed in Ref. [67] that one reproduces exactly the LS-equation solutions for regular potentials. This is also true for the singular potentials when the potential is used in the whole range of integration in the LS equation, that is, for q ∈ [0, ∞] (the cut-off is sent to infinity). For the singular-potential case we refer to the standard kind of solutions, so that for a repulsive singular potential the solution has no free parameters and is determined, while for the attractive singular case the solution involves one free parameter that could be fixed e.g. by imposing a given value for the scattering length [97][98][99]. Several potentials were studied in Ref. [67], both for uncoupled and coupled PWAs. Within the latter group the 3 S 1 -3 D 1 coupled PWAs were studied and the 3 S 1 scattering length was taken as input. Needless to say, in all cases the LS equation with infinite cutoff and the N/D method agree perfectly in our numerical study.
The fact of having none or only one free parameter is a very constrained situation in practical applications, and it is the reason why it has not been possible to achieve yet a good agreement with data in NN scattering in terms of regulator-independent solutions (i.e. in which the three-momentum cut-off is taken to infinity). Notice that the number of free parameters in the solution of the LS equation for singular potentials is then not linked with the chiral order in the calculation of the chiral potential. However, in terms of the N/D method one can in principle add an arbitrary number of subtractions, which allows one to explore for extra solutions. We have already discussed this point in connection with the ambiguity associated with the CDD poles in Sec. 5.1. This possibility was explored in detail in Ref. [68] for the 1 S 0 NN PWA. The NLO and NNLO ChPT potentials for this PWA are actually attractive and singular. The standard solutions of the LS equation were reproduced, and a detailed numerical analysis was performed in order to show the agreement between the LS equation and the exact N/D method. But we also showed in this reference that one can also generate new solutions that cannot be achieved by the LS equation when the three-momentum cut-off is taken to infinity with contact interactions included in the potential to aim renormalization (in the form of polynomial counterterms in its momentum expression). In this way, a new solution was discussed that can reproduce the 1 S 0 scattering length, effective range and shape parameter v 2 . For this solution the DRs for N (s) and D (s) converge separately. It is also interesting to indicate that a solution within the exact N/D method for this PWA reproducing the scattering length and the effective range could not be found. Last but not least, a very attractive feature of the exact N/D method is that it allows to evaluate the PWAs in whole complex p 2 plane. Then, it is very convenient to look for resonance and (anti)bound states. In the case of the 1 S 0 PWA there is an antibound state which is found at p = −i0.066 MeV both at NLO and NNLO when all the first three ERE parameters are reproduced.

Concluding remarks
We have elaborated on several unitarization methods of perturbative calculations in ChPT that can be employed to study scattering and the re-scattering corrections to an external probe. Special attention has been given to the N/D method both for scattering and FSI implementation. The unitarization methods, since the earlier papers on current algebra techniques, have been able to extend to much larger energies the expected region of utility of ChPT calculations. This has been accomplished thanks to the extra energy and momentum dependence generated by using a non-perturbative theoretical framework which satisfies key properties of S-matrix theory, which stem from two-body unitarity and analyticity. Some of the most striking and important applications of the unitarization methods of input perturbative calculations has occurred in the field of spectroscopy. In this way it has been possible to study resonances and bound states, and even predict some of them, while unexpected properties have been unveiled too, as e.g. the two-pole nature of some resonances [100], as first shown for the Λ(1405) in Ref. [18].