Solving integral equations in $\eta\to 3\pi$

A dispersive analysis of $\eta\to 3\pi$ decays has been performed in the past by many authors. The numerical analysis of the pertinent integral equations is hampered by two technical difficulties: i) The angular averages of the amplitudes need to be performed along a complicated path in the complex plane. ii) The averaged amplitudes develop singularities along the path of integration in the dispersive representation of the full amplitudes. It is a delicate affair to handle these singularities properly, and independent checks of the obtained solutions are demanding and time consuming. In the present article, we propose a solution method that avoids these difficulties. It is based on a simple deformation of the path of integration in the dispersive representation (not in the angular average). Numerical solutions are then obtained rather straightforwardly. We expect that the method also works for $\omega\to 3\pi$.


Introduction
The study of the η → 3π decay process is interesting, first and foremost, in the context of the determination of the quark mass ratio In order to extract the value of Q to high precision, it is very important to have a robust control on the final-state interactions in this decay, which lead to a strong effect in the width. To this end, one would like to have a non-perturbative framework, allowing the resummation of a certain class of the final-state interactions to all orders. Dispersion relations are the ideal tool for this. Incorporating 2-particle unitarity and crossing symmetry then leads to a system of coupled integral equations for the 3 isospin amplitudes in this decay. It took quite some time until these equations were written down in their final form. The development started with the pioneering work of Khuri and Treiman [1]. The mathematical structure of this type of equation was investigated in the following decade [2][3][4][5][6][7][8][9][10], see also the monograph [11]. Later, interest in the dispersive method waned. A revival occurred in the nineties, when it was demonstrated [12][13][14][15] how dispersion relations allow one to incorporate final state interactions of S-and P-waves in a reliable and calculable manner. Refs. [13,14] contain a detailed discussion of the role of subtractions in the dispersive representation, while the uniqueness of the solutions is investigated in [14]. Further, current algebra [16] and chiral perturbation theory results [17] were used to relate the parameter Q to this framework (normalization of the dispersive amplitude), and to get a handle on the subtraction constants. For a review of the early developments until 1990 we refer the reader to Ref. [12], see also the lecture notes [18]. An improved representation of the ππ phase shifts [19,20], the evaluation of electromagnetic corrections [21,22] as well as new experimental information on this decay [23][24][25][26][27][28][29][30] triggered new investigations of η → 3π in the dispersive framework [31][32][33][34][35][36][37][38][39][40][41]. A very comprehensive analysis has recently been presented in Ref. [42]. For applications of the dispersive framework in other three-body decays, see Refs. [34,35,[43][44][45][46][47]. Further attempts to incorporate a class of final state interactions in the η → 3π amplitudes may be found in Refs. [17,[48][49][50][51][52]. The present article is devoted to a discussion of numerical methods to solve the integral equations that occur in the formulation of Ref. [42]. In this connection, the investigations [2][3][4][5][6][7][8][9][10] revealed the following two important aspects: i) In the evaluation of the angular averages of the amplitudes, the integration path in z = cos θ must be deformed into the complex plane, in order not to destroy the holomorphic properties of the amplitudes. The deformation is fixed by providing the (eta mass) 2 with an infinitesimal positive imaginary part: M 2 η → M 2 = M 2 η + iδ , with δ → 0 + at the end of the calculation.
ii) Singularities emerge in the angular averaged amplitudes, at the pseudothreshold s = (M η − M π ) 2 . Some of these singularities are of a non-integrable type in the Lebesgue sense. The (positive) imaginary part of the eta mass, however, acts as a regulator, that can be removed after the dispersive integral has been performed. Handling these singularities properly is a delicate affair, as is illustrated e.g. by the discussions in Refs. [13,15,31,35,41,42,45,52].
We shortly mention several previous methods to obtain solutions of these equations. In Refs. [13, 15, 31-33, 42, 43], the equations were solved directly through iterations, with a careful treatment of the mentioned singularities. In Ref. [35,41], integral kernels are introduced, which make the numerical procedure faster, and which identify the mentioned singularities at the pseudothreshold in a clear fashion. Through the interchange of the order of integrations, the so-called Pasquier inversion [9,10], it is possible to carry out one integration and to obtain integral equations in one variable that are further solved by iterations (see Refs. [36,39,53,54] for more details). Finally, it is possible to obtain integral equations for the angular averaged amplitudes instead of the amplitudes themselves. This technique has recently been invoked in Refs. [44,45]. Last but not least, the convergence of the iterative procedure to solve the equations is not a priori guaranteed. The iterative procedure converges very well (typically, after 3-4 iterations) in the η → 3π decays, but the convergence becomes an issue for decaying particles with heavier masses, or for equations with more subtractions [45]. For this reason, e.g., in Ref. [44], a matrix inversion method was used to find a solution beyond the iterative procedure. These enterprises are numerically very demanding, in particular for the reasons spelled out above.
In the following, we propose a numerical framework that avoids the difficulties i) and ii) altogether. We show that one may deform the integration path in the dispersion integral (not in the angular average) into the complex plane. After choosing a properly deformed path, the integrand becomes regular, even at δ = 0 (except at threshold, where a mild, integrable square root singularity persists). In addition, the angular integration can be carried out in the original interval −1 ≤ z ≤ 1. Further, discretizing the integrals through a Gauss-Legendre quadrature, the integral equations are transformed into a set of linear matrix equations, whose solution is easily found by iterations. The procedure is pretty straightforward, takes very little CPU time and does not lead to any spurious irregularities in the solutions 1 .
The idea to avoid singularities in integration through path deformation is not new. An early reference is the work of Hadamard [56], whose method was rediscovered 50 years later by investigation of analytic properties of perturbation theory [57,58], and then extended and extensively used in S-matrix theory [58]. The analogue of the procedure for the quantum-mechanical three-body problem is well known since decades [59,60]. As already said, Aitchison and Pasquier mentioned this possibility for Khuri-Treiman-type equations in Ref. [55]. For the numerical evaluation of loop graphs in Quantum Field Theory, path deformations in momentum space [61] as well as in Feynman parameter space [62] are used. All these techniques, including the present one, may be summarized under the heading "Applications of Cauchy's integral theorem".
The layout of the article is as follows. In section 2, we display the integral equations for the η → 3π amplitudes in the form worked out recently in Ref. [42]. We describe in some details the technical difficulties that one encounters while solving the equations in section 3, while in section 4, we define the deformed path in the dispersion integral and show that in this manner, the singularities in the integral equations disappear. In section 5 we describe the numerical procedure for solving the equations, whereas a summary and conclusions are given in section 6. In appendix A we comment on the phase shifts used and collect some notation. The holomorphic continuation of the integrand is discussed in appendix B, path deformations in general are investigated in appendix C, and the case ω → 3π is shortly discussed in appendix D.

The equations
We start from the integral equations for the three isospin amplitudes M I (s) in the framework specified in [42], The various quantities are defined as follows. The Omnès functions are where δ I (x) denote the elastic ππ phase shifts (S-wave for I = 0, 2, P-wave for I = 1). The hat-functionsM I are defined in terms of angular averages, where Finally, the measures are The C m II ′ (x, κ) are polynomials in x and κ, and P I (s) denote subtraction polynomials, see [42]. The phase shifts δ I are needed as input to solve these equations. We relegate a discussion of them to the appendix A, where we also collect some of the notation used below.

Technical hurdles
During the numerical, iterative solution of the equations (2)-(6), two main problems occur.
i) The angular averages z m M I (x) amount to an integration along the path h(x, z) of the argument of the amplitudes M I . For fixed x, the argument runs along a straight line, from The Kacser function κ(M 2 , M 2 π , x) [5] is holomorphic in the complex x-plane, cut along the real axis from x = 0 to x = 4M 2 π and along a straight line from Fig. 1 2 . Therefore, fixing its value at some point in the complex x-plane renders it unique in the cut plane. In the following, we use the convention In Fig. 2, we display the endpoints h ± (x) as the integration variable x runs from the threshold to infinity. Because κ(M 2 , M 2 π , x) is holomorphic in the cut x plane, the It is seen that this line crosses the original path of integration in the dispersive representation (2). The same happens for x 1 < x < x 2 . See text for details. endpoints h ± move -for δ > 0 -along curves that are infinitely often differentiable with respect to the real variable x. The dash-dotted vertical line connects the endpoints h ± (x) with x 2 < x < x 3 . It is seen that a straight line between h − (x) and h + (x) crosses the integration path in the dispersive representation (2). This is also the case for x 1 < x < x 2 . The problem is discussed at many places in the literature -early references are [2,4]. Its solution amounts to deform the path in the z integration (or in the integration over the variable h after a change of variables z → h), such that this crossing is avoided. We do not discuss this point any further here.
ii) The second problem arises after the angular integration has been performed in one way or the other. The angular averages develop singularities of the typê where p = ( 1 2 , 3 2 , 1 2 ) for I = (0, 1, 2), and K denotes a constant. These singularities then show up in the dispersive integrals (2). It turns out that, after the integration over x has been performed, the limit M → M η exists. [See e.g. Refs. [13,42], where the dispersive integral has been worked out explicitly.] We note that the case p = 3 2 corresponds to a non-integrable singularity at δ = 0: one is not allowed to interchange the limit δ → 0 + with the integration over x. These singularities render the standard procedure to construct a numerical solution of the integral equations The two singularities of the integrand are indicated with filled circles. Since these are located in the upper complex plane, the path may be deformed into the lower half plane, e.g., into the polygonal line rather delicate and cumbersome. In the following, we present a method to solve these equations in an easy and straightforward manner, that avoids the problems i) and ii) altogether.
[Including higher partial waves leads to even stronger singularities [63]. Our method also covers these cases, see the following section.]

Avoiding singular integrals
There is a simple way to cope with the singularities of the angular averages. To illustrate, we consider the integral [13,42] G(s) = We display in Fig. 3 the singularities of the integrand with two black dots in the complex x-plane: We have also drawn a cut that emerges from the branch point at this second singularity, reaching out to x = ∞ + iδ (dotted line). The path of integration is indicated with a solid line, from x 1 to x 3 . We now observe that the integrand is holomorphic in the complex half-plane Im(x) < ǫ, δ. Therefore, we may deform the path of  Fig.(3).
integration into the polygonal line x 1 B 1 C 1 X 3 (dashed), without changing the value of the integral. There are then no singularities anymore on the integration path, and we may set ǫ = δ = 0 before performing the integration. This proves that the limit ǫ, δ → 0 + of the original integral (10) exists [13,42]. [It is obvious that this remains true if the exponent 3/2 in (10) is replaced by any p ∈ C.] On the other hand, approaching the real axis from below, one encounters a pinch singularity at s = (M − M π ) 2 , which results in a singular behaviour of the function G(s).
A numerical integration along the dashed path does not pose any problems, because the integrand is smooth as a function of x. To illustrate, we display in Fig. 4 the real part (solid line) and imaginary part (dashed line) of the function G(s) (analytic expression is given in [42]), in units of the pion mass. For comparison, the result of the numerical integration along the dashed polygonal line with 160 Gauss-Legendre points is shown with filled circles. It is seen that the agreement is perfect.
We now use the same method in the original equations (2). Suppose that the integrand in (2) is holomorphic in some region Im(x) < 0 (we discuss this point in appendix B). Then we can avoid the singularities generated by the zeros in the function κ(M 2 , M 2 π , x) by deforming the path as shown in Fig 1. Aside from avoiding the singularities inM I (x), this procedure has the following advantage. Consider the endpoints of the angular integration as they now occur on the deformed path, see Fig. 5. There, we display the two endpoints h ± (x). It is seen that the problem with leaving the holomorphy domain of M I (s) does not occur anymore. See also Fig. 6, where we display some of the paths h(x, z), − 1 ≤ z ≤ 1, that connect h ± (x).  (7), if the integration in the dispersive integral (2) is performed along the deformed path ABCDEF in Fig. 1. Compare with Fig. 2, which shows the situation in the standard case.
We conclude that one can solve the equations (2) in their original form, avoiding complicated path deformations and singular integrals, provided that the path of integration in (2) is deformed properly -e.g. according to Fig. 1.
The deformed path displayed in Fig. 1 is obviously not the only one with these properties -see appendix C for a discussion of this point.

Solving the integral equations
Finally, we make use of the fact that now, the integrands are smooth along the path of integration for s < D, except at the threshold s = 4M 2 π , where an integrable singularity of the type const/ x − 4M 2 π occurs (note thatM I (x) → const. as x → 4M 2 π ). The singularity can be tamed with a variable transformation z = A + (B − A)τ 2 , 0 ≤ τ ≤ 1. Therefore, an integration using the Gauss-Legendre method [64] is adequate. As we now show, this method has the further advantage that the integral equations boil down to a matrix equation, with matrix elements that need to be evaluated only once, before the iteration, which then becomes trivial.
We use Gauss-Legendre quadrature both for integration over x and over z. The integration path over x is split into linear pieces: AB, BC, CD and DF (the point F corresponds to the upper limit of integration: the phase shifts δ I (x) are equal either to 0 or to π, if x moves right to F ). The variable x from each interval is mapped to the interval [0, 1], and the Gauss-Legendre mesh points x i and weights w x i are used to carry out the integration. The number of points on these intervals is N AB , N BC , N CD and N DF , respectively. Further, the integral over z for a given x always runs from −1 to 1. The mesh points and weights are denoted by z a and w z a , respectively, and the number of mesh points is chosen to be L AB , L BC , L CD and L DF for x belonging to the one of the above intervals.
It is useful to introduce a multi-index, Here, Combining the equations (2) and (4), we may write where Finally, letting the variable s run over the set h α , α = 1, . . . N and introducing the notations R βα II ′ = w α Ω I (h β )K II ′ (h β , h α ) and G β I = Ω I (h β )P I (h β ), we can write down Eq. (13) as a matrix equation This equation generates iteration series for the vectors M β I , I = 0, 1, 2, which are rapidly convergent for η → 3π (note that the quantities R βα II ′ need to be evaluated only once for a given set of phase shifts). In the interval s ∈ (−∞, D), the amplitudes can then be directly constructed from M β I by using the dispersive representation (2). For s > D, the angular averaged amplitude must be interpolated to perform the Cauchy integral. Moreover, if one manages to invert the large matrix R βα II ′ numerically, Eq. (15) can be solved directly, without iterations, even if the latter do not converge [44]. Because there was no need to do so in our case, we sticked to the iteration procedure.
Finally in Fig. 7, we display the isospin zero amplitude for one particular choice of the subtraction polynomials, and for a specific choice of the ππ phase shifts. The amplitude agrees with the one constructed by Lanz [65], except near the threshold, where the cusp structure is different, and near the pseudo-threshold in the I = 1 channel.
6 Summary and conclusions 1. We have considered numerical aspects of the integral equations Eqs.(2)- (6) that govern η → 3π decays [42]. The standard approach to solve these equations numerically is confronted with two main technical hurdles: angular averages along complicated paths in the complex plane, and singularities of the integrand near the integration path in the dispersive representation.
2. Holomorphicity of the phase shifts in the low energy region allows one to deform the original path of integration in the dispersive representation (not in the angular averages). Both problems disappear: The angular averages can be performed in their original form, and there are no nearby non-integrable singularities on the integration path.
3. As the integrands are smooth, a Gauss-Legendre integration becomes feasible. The integral equations turn into a matrix equation, whose iterative solution can be obtained very efficiently. 4. We have constructed the 6 fundamental solutions [42] of Eqs.(2)-(6) for several sets of phase shifts, in the region 0 < s < 20 M 2 π . The results generally agree with the solutions obtained by Lanz [65]. 5. We provide tables with the fundamental solutions for 8 different sets of phase shifts as ancillary files.
6. As we show in appendix C, the very same technique is expected to also work for ω → 3π [43]. It remains to be seen to what extent it can be applied to other three-body decays and to the calculation of form factors [34,35,[43][44][45][46][47]. on the manuscript and for information on the ππ phase shifts. We thank Stefan Lanz for providing us with his Fortran code to construct the numerical solutions using the standard approach, and for data files with his fundamental solutions. A.R. acknowledges the support from the DFG (CRC 110 "Symmetries and the Emergence of Structure in QCD"), from Volkswagenstiftung under contract no. 93562 and from Shota Rustaveli National Science Foundation (SRNSF), grant no. DI-2016-26. J.G. thanks the HISKP at the University of Bonn for warm hospitality. Part of this work was performed during his stay there.

A Notation and phase shifts
In the text and in figures, we use the symbols The new integration path is fixed by the vertices In section 4, we set In the case ω → 3π we use The elastic S-wave (P-wave) ππ phase shifts δ 0,2 (δ 1 ) are needed as an input to solve (2) - (6). We rely on the phase shifts used in Ref. [42]: In the low energy region x < E, we use a Schenk parameterization [66]. Above x = F , the phase shifts are set to a constant, For E ≤ x ≤ F , we use phase shifts from [42]. Tables with these phase shifts at discrete values of energies between E and F, as well as the Schenk parameters to describe the phase shifts below E are provided as ancillary files, together with the pertinent basis functions.

B Deforming the path of integration
Here, we wish to show that the original path of integration in (2) can be deformed into the lower complex x-plane. For this to achieve, we must prove that the integrand is holomorphic there. As an input, we need the phase shifts δ I (x) and a starting value of the hat-functionsM I (x). For the latter, we use the current algebra expressions, which are polynomials in x and thus even entire functions. Concerning the phase shifts, see appendix A. The Schenk parameterization used below 800 MeV has the virtue that the corresponding expression for the phase shift is holomorphic in a strip of the complex x-plane, cut for x < 4M 2 π and x > 700 MeV. This is so because the Schenk parameterization provides a holomorphic expression for the tangent of the phase shifts. The phase shift itself is then obtained via As long as tan δ I (x) = ±i, δ I (x) is holomorphic. We display in Fig. B.1 the complex quantities tan δ I (x) along the polygonal path ABCD displayed in Fig. 1, and along similar paths which are nearer to the real line. It is seen that one does not hit any singularity. We have checked that the parameterization in terms of conformal variables as worked out in Ref. [69] leads to the same conclusion. We conclude that sin δ I is holomorphic in that region as well. Consider now the denominator of the measure, |Ω I (x)|. At first, it sounds surprising that one can  [67,68]). It is seen that the singularity is never hit in the region spanned by the path ABCD with the choice (A.2) for the 4 vertices. Schenk parameterization [66] is used, with parameters from [19], see appendix A.
continue analytically the modulus of a complex function. The point is that this quantity stands for where − denotes a principal value integral. The relevant quantity is thus the integral It is useful to first consider the related ordinary integral The function H I (x) is holomorphic in the complex x-plane, cut along the real axis for x ≥ 4M 2 π . On the upper rim of the cut, G I and H I are closely related: Consider now the interval [A, D] in Fig. 1. We can continue here H I (x) across the cut, because δ I (x) are holomorphic there, and the integration path D A can be deformed into the lower half plane, e.g., into the polygonal line ABCD in Fig. 1.
The relation (B.5) shows that G I and thus |Ω(x)| can holomorphically be continued as well. The continuation is unique.
In practise, the continuation can be performed rather easily: One adds and subtracts the phase shift δ I (x) in the integrand of G I (x) and obtains This identity holds for x ∈ [4M 2 π , ∞). In the integral on the right hand side, one need not use the principal value prescription. The right-hand side is holomorphic in the region where the phase shift is holomorphic.
Finally, the Omnès functions Ω I (x) which enters the iterations in the region covered by the dashed blue lines in Fig. 6 may be evaluated in a manner analogous to |Ω I (x)| just discussed: Below the threshold, use the integral representation (3), as well as above threshold, for Im(x) > 0. In the case of Im(x) < 0, Re(x) > 4M 2 π , proceed as in the case of the principal value integral just discussed.
We conclude that we may indeed deform the path of integration as is indicated in the Fig. 1. [In order to render the proof watertight in view of the square-root singularities at the threshold x = 4M 2 π , one may slightly change the paths. Let A ′ = 4M 2 π + ǫ, ǫ > 0. Then ABCDF→ AA ′ BCDF. Numerically, this change is irrelevant, and we stick to ǫ = 0. ] We have checked that the results do not change under a change of the polygonal line ABCD.
The following two remarks are in order at this place.
1. We do not assume that the phase shifts are holomorphic at all energies. As mentioned, we glue the phase shifts continuously together, using different parameterizations. Only in the low-energy region do we rely on a holomorphic parameterization.
2. One might argue that we extrapolate data from the real line to the polygonal line ABCD, which would be an unstable procedure. Instead, we make use of a holomorphic parameterization to render the integrations easier to perform. For a given parameterization, the result is the same, whether performed in the standard manner, or using a deformed path in the dispersion relation. The latter method is, however, by far more efficient.

C Allowed paths
Here, we discuss paths in the dispersive representation that allow for an undistorted path in the angular average, and follow the discussion in [55]. For fixed x, the angular average amounts to a line integral between h ± (x). Starting from x = ∞ to lower values, it is clear from Fig. 2 that a critical value of x is reached as soon as this straight line touches the value 4M 2 π for some z ∈ [−1, 1], Deformed integration path used in the present work (displayed for s < 30M 2 π ). Any deformed path that does not cross the solid line allows for an undistorted integration in the angular average (4).
This condition may be brought to the form For a fixed value of z = ±1, the equation (C.2) has one real and 2 complex conjugate solutions for x. Varying z in the interval [−1, 1] traces out a curve in the complex x-plane. In Fig. C.1, we display the curve with Im(x) ≤ 0 (solid black line). If one uses a path in the dispersion integral that does not cross the hatched region, the angular average can be performed in the undeformed interval −1 ≤ z ≤ 1. The polygonal line ABCDEF (displayed for s < 30M 2 π ) is the version used in the present work. See appendix B for the actual values chosen for A, . . . , F and for x i .
The following two remarks are at order at this place. Dashed: A deformed integration path. Any deformed path that does not cross the solid line allows for an undistorted integration in the angular average (4). For notation see appendix A. Note that √ D 2 ≃ 0.75 GeV, i.e., the Schenk parameterization is still operative at D 2 and generates holomorphic phase shifts.
i) The reason to extend the polygonal line into the lower complex x-plane until x = D is the following. After having solved the integral equation for the angular averages, one obtains the physical amplitudes by performing the final dispersive integration with Cauchy kernel 1/(x − s − iǫ). Even at ǫ = 0, this kernel is then not singular on the integration path for real values of s with s ∈ [D, F ], except at the threshold s = 4M 2 π . There, the integrand in (2) develops an integrable square root singularity in case of S-waves. It can be tamed with a transformation of variables, z = A + (B − A)τ 2 , 0 ≤ τ ≤ 1. The integration can thus be performed at ǫ = 0, without further ado -an additional advantage of our procedure. For s ∈ [D, F ], due to the singularity of the Cauchy kernel, the angular average needs to be interpolated between the pertinent Gauss-Legendre points, before the integration can be done. This is the reason why we push D to higher values than required.
ii) In order to stay away from the threshold region x = 4M 2 π while performing the angular average, we have chosen the path BC with an imaginary part that is more negative than what is needed according to the condition (C.2), see also Figs. 5 and 6.

D The decay ω → 3π
We shortly consider the situation for ω → 3π decays [43]. In Fig. C.2, we display the analogue to the boundary in Fig. C.1, for M 2 η replaced by M 2 ω . A possibility for the integration path is displayed with dashed lines. It is seen that with Schenk parameterization [66], it is still possible to avoid the pseudohreshold singularity and to perform the dispersive integral without interpolation between Gauss-Legendre points in the interval s < D 2 , which contains the decay region. We therefore expect that the method also works in this case.