Data-driven dispersive analysis of the $\pi \pi$ and $\pi K$ scattering

We present a data-driven analysis of the resonant S-wave $\pi\pi \to \pi\pi$ and $\pi K \to \pi K$ reactions using the partial-wave dispersion relation. The contributions from the left-hand cuts are accounted for using the Taylor expansion in a suitably constructed conformal variable. The fits are performed to experimental and lattice data as well as Roy analyses. For the $\pi\pi$ scattering we present both a single- and coupled-channel analysis by including additionally the $K\bar{K}$ channel. For the latter the central result is the Omn\`es matrix, which is consistent with the most recent Roy and Roy-Steiner results on $\pi\pi \to \pi\pi$ and $\pi\pi \to K\bar{K}$, respectively. By the analytic continuation to the complex plane, we found poles associated with the lightest scalar resonances $\sigma/f_0(500)$, $f_0(980)$, and $\kappa/K_0^*(700)$ for the physical pion mass value and in the case of $\sigma/f_0(500)$, $\kappa/K_0^*(700)$ also for unphysical pion mass values.


Introduction
There is a renewed interest in the hadron spectroscopy, motivated by recent discoveries of unexpected exotic hadron resonances [1,2,3]. Currently, LHCb, BESIII, and COMPASS collected data with unprecedented statistics, BELLE-II and GlueX just started to operate, and more facilities are planned in the near future, such as PANDA and EIC. Besides, lattice QCD has been applied to a broad range of hadron processes and recently was able to calculate the lowest excitation spectrum with the masses of the light quarks near their physical values [4].
To correctly identify resonance parameters one has to search for poles in the complex plane. This is particularly important when there is an interplay between several inelastic channels or when the pole is lying very deep in the complex plane. In these cases, the structure of the resonance is quite different from a typical Breit-Wigner behaviour. In order to determine the pole position of the resonance, one has to analytically continue the amplitude to the unphysical Riemann sheets. At this stage, the right theoretical framework has to be applied. The latter should satisfy the main principles of the S-matrix theory, namely unitarity, analyticity, and crossing symmetry. These constraints were successfully incorporated in the set of Roy or Roy-Steiner equations [5]. In a practical application, however, the rigorous implementation of these equations is almost impossible, since it requires experimental knowledge of all partial waves in the direct channel and all channels related by crossing. Therefore, the current precision studies of ππ [6,7,8,9] and πK [10,11] scattering are based on a finite truncation, which in turn limits the results to a given kinematic region, and require a large experimental data basis. Furthermore, applying Roy-like equations for coupled-channel cases is quite complicated and has not been achieved in the literature so far. Because of the above-mentioned difficulties, in the experimental analyses, it is a common practice to ignore the S-matrix constraints and rely on simple parameterizations. The most used ones are a super-position of Breit-Wigner resonances or the K-matrix approach. The latter implements unitarity, but ignores the existence of the left-hand cut and often leads to spurious poles in the complex plane.
A good alternative to the K-matrix approach and a complementary method to Roy analysis is the so-called N/D technique [12], which is based on the partial-wave dispersion relations. In this method, the dominant constraints of resonance scattering, such as unitarity and analyticity are implemented exactly. Since the time it was introduced by Chew and Mandelstam [12], the N/D method has been extensively studied for different processes [13,14,15,16]. The required input to solve the N/D equation are the discontinuities along the left-hand cut, which are typically approximated one way or another using chiral perturbation theory (χPT). In the present paper, we extend the ideas of [15], where the left-hand cut contributions were approximated using an expansion in powers of a suitably chosen conformal variable. In contrast to [16], however, we follow here a data-driven approach and adjust the unknown coefficients in the expansion scheme to empirical data directly. In this way, the model dependence is avoided, and the method can also be applied to the reactions which do not include Goldstone bosons, like for instance the J/ψJ/ψ scattering [1].
• Even though the ππ → ππ (and to a lesser extent πK → πK and ππ → KK) amplitudes are known very well from the Roy (Roy-Steiner) analyses [7,8,9,10,11,30], in the practical dispersive applications the final state interactions (FSI) are implemented with the help of the so-called Omnès function, which does not have left-hand cuts. Indeed, the left-hand cuts are different for each production/decay mechanism, while the unitarity makes a connection between the production/decay and the scattering amplitudes only on the right-hand cut. In the N/D ansatz, the Omnès functions come out naturally, as the inverse of the D-functions.
• Recently, it has become possible to calculate ππ and πK scattering using lattice QCD with almost physical masses [31,32,33,34,35,36,37]. Since, both the σ/ f 0 (500) and κ/K * 0 (700) states lie deep in the complex plane, the reliable extraction of their properties requires the use of the formalism that goes beyond the simple K-matrix parametrization and incorporates in addition the analyticity constraint.
The paper is organized as follows. In the next section, we focus on the formalism that we adopt in this paper. We start with the review of the N/D method in Sec 2.1. We then discuss the left-hand cut contributions in Sec 2.2. In Sec 2.3 we make the connection to the Omnès functions. The numerical results are presented in Sec. 3. We start with I = 0, ππ singlechannel analysis of both experimental and lattice data, which is followed by the coupled-channel {ππ, KK} analysis of the experimental data. These results are then used to determine the two-photon coupling of σ/ f 0 (500) and f 0 (980). At the very end, we focus on πK, I = 1/2 scattering of both experimental and lattice data. A summary and outlook is presented in Sec. 4.

N/D method
The s-channel partial-wave decomposition for 2 → 2 process is given by where θ is the c.m. scattering angle and ab are the coupledchannel indices with a and b standing for the initial and final state, respectively. For the following discussion, we focus only on the S-wave (J = 0) and therefore will suppress the label (J). The different normalization factors (N ππππ = 2, N ππKK = √ 2 and N KKKK = N πKπK = 1) are needed to ensure that the unitarity condition for identical and non-identical two-particle states are the same and can be written in the matrix form as where the sum goes over all intermediate states. The phase space factor ρ c (s) in Eq. (2) is given by with p c (s) and s th being the center-of-mass three momentum and threshold of the corresponding two-meson system. Within the maximal analyticity assumption [38], the partial-wave amplitudes satisfy the dispersive representation where s L is the position of the closest left-hand cut singularity and the discontinuity along the right-hand cut is given by (2). For unequal masses, as in πK scattering, the left-hand singularities of the partial-wave amplitude do not all lie on the real axis and the integration in the first term in Eq. (4) goes partly along the circle. We note, that the separation into left and righthand cuts given in (4) is only possible for the systems where no anomalous thresholds are present [39,40]. The unitarity condition (2) guarantees that the partial-wave amplitudes at infinity approach at most constants. In accordance with that, we can make one subtraction in Eq. (4) to suppress the high-energy contribution under the dispersive integrals. Thus we rewrite Eq. (4) as where we combined the subtraction constant together with the left-hand cut contributions into the function U ab (s). The choice of the subtraction point s M will be discussed later. The solution to (5) can be written using the N/D ansatz [12] t ab (s) = where the contributions of left-and right-hand cuts are separated into N(s) and D(s) functions, respectively. The discontinuity relation along the right-hand cut Disc D ab (s) = −N ab (s)ρ b (s) allows us to write a dispersive representation for the D-function, which up to a Castillejo-Dalitz-Dyson (CDD) ambiguity [41] 1 is given by Due to the non-uniqueness of the N/D ansatz, we have normalized the D-function in Eq. (7) such that D ab (s M ) = δ ab . Since D ab (s) is a complex matrix above the threshold, the position of s M has to be chosen such that all of its elements are real at this point, i.e. s M ≤ s th . To arrive at an integral equation for the N(s) function, one can write a once-subtracted dispersion relation for c D ac (s) (t(s) − U(s)) cb and fix its subtraction constant by requiring that which follows from Eq. (5). As a result, it yields [43] N ab (s) = U ab (s)+ (9) The above integral equation can be solved numerically given the input of U ab (s). Knowing the N ab (s) function on the righthand cut, the D ab (s) function is calculated by (7) and finally the partial-wave amplitude is produced with Eq. (6). In other words, if the discontinuities across all the left-hand cuts were known 2 the exact solution can be obtained by N/D method. An important property of Eq. (9) is that the input of U(s) is only needed on the right-hand cut. In the case of many channels, both the diagonal and off-diagonal t-matrix elements have a right-hand cut starting at the lowest threshold s th . However, only the input of the off-diagonal U ab (s) is required outside the physical region, while in order to solve (9), the input of the diagonal U aa (s) is needed in the physical region due to the phase space factor. It has a direct relevance for the {ππ, KK} case, where in the KK → KK channel the overlap of left-and right-hand cuts happens, but only in the non-physical region, 4m 2 π < s < 4(m 2 K − m 2 π ), and therefore does not require any modifications of the dispersion integrals. We also emphasize that by means of Eq. (6), the scattering amplitude can be rigorously continued into the complex plane, where one can determine pole parameters of the resonances. In our convention the scattering amplitude in the vicinity of the poles on the unphysical Riemann sheets (or physical Riemann sheet in the case of the bound state) is given by, where the normalization factor N ab comes from Eq. (1) and g pi denotes the coupling of the pole at s = s p to the channel i = a, b. We wish to comment on the case when there is a bound state in the system, since it happens for the relatively large unphysical pion masses. To find the binding energy s B , one searches for a zero of the determinant of the D ab matrix for energies below threshold, det(D ab (s B )) = 0, s B < s th .
In this case, the solution obtained using the set of N/D equations (6) with input from (13) satisfies the dispersion relation (5) combined with the bound state term, At the same time, it is straightforward to show that including such a bound state term into the definition of U ab (s) does not change the solution of (6) or the integral equation (9), provided that the residues g Ba g Bb are dialed properly using the det(D ab (s B )) = 0 condition.

Left-hand cuts
In a general scattering problem, little is known about the lefthand cuts, except their analytic structure in the complex plane. 2 in that case the subtraction constant is probably unnecessary to introduce. The progress has been made in [15], by considering an analytic continuation of U ab (s) to the physical region, which is needed as input to Eq. (9), by means of an expansion in a suitably contracted conformal mapping variable ξ(s), which is chosen such that it maps the left-hand cut plane onto the unit circle [44]. The form of ξ(s) depends on the cut structure of the reaction (i.e. {ab}) and specified by the position of the closest left-hand cut branching point (s L ) and an expansion point (s E ) around which the series is expanded, ξ(s E ) = 0. Since for the {ππ, KK} system all the left-hand cuts lie on the real axis, −∞ < s < s L , one can use a simple function where s L (ππ → ππ) = s L (ππ → KK) = 0 and s L (KK → KK) = 4 (m 2 K − m 2 π ). For the case of πK → πK, the left-hand cut structure is a bit more complicated (see Fig. 1). In addition to the left-hand cut lying on the real axis −∞ < s < (m K − m π ) 2 , there is a circular cut at |s| = m 2 K − m 2 π . The conformal map that meets these requirements is defined as where s L (πK → πK) = m 2 K − m 2 π . We note that, given the forms of ξ(s) in Eqs. (14) and (15), the series (13) truncated at any finite order is bounded asymptotically. This is consistent with the assigned asymptotic behavior of U(s) in the once-subtracted dispersion relation (5).
For reactions involving Goldstone bosons, in principle, χPT allows to calculate the amplitude over a finite portion of the closest left-hand cut and can be used to estimate C n in (13) as it has been done for other processes in [15,16]. However, it is not clear at which point χPT calculated to a given order still represents a good approximation. In addition to that, in order to merge the conformal expansion with the chiral expansion, the expansion point s E should lie within the region where χPT can be computed safely. For instance, for the elastic ππ → ππ scattering the natural choice would be to identify s E with the two-pion threshold. However in that case, the last data point, which can be described with the elastic unitarity, corresponds to ξ(s 1/2 max = 0.7 GeV) 0.45. On the other side, the faster convergence of the sum in Eq. (13) can be achieved for the choice of s E in between the threshold and s max , i.e. in the regions where χPT is at the limit of its applicability. Besides, for the coupledchannel case, one needs to rely on SU(3) χPT, which converges slower than SU(2) version of it.
In our paper, we determine the unknown C n in Eq. (13) and the optimal positions of s E directly from the data and use χPT results only as constraints for the scattering lengths, slope parameters, and Adler zero values. We note that the latter brings a stringent constraint on the scattering amplitude, since for both ππ and πK scattering the Adler zero is located very close to the left-hand cut (see Fig.1), and cannot be determined precisely from the fit to the data. However, once the Adler zero is imposed as a constraint, it improves drastically the convergence of (13) in the threshold region.

Relation to the Omnès function
The unitarity connects the partial-wave amplitudes in production (or decay) and scattering processes. Therefore, the reactions like γp → ππp, γγ → ππ, J/ψ → ππγ, η → 3π, etc. are very sensitive to the FSI. In a dispersive formalism, FSI are typically implemented with the help of the so-called Omnès function [45], Ω ab (s), that fulfills the following unitarity relation on the right-hand cut and analytic everywhere else in the complex plane, i.e. it satisfies a once-subtracted dispersion relation Therefore, for the case of no bound states or CDD poles, the D ab (s) function obtained in (7) can be easily related to the Omnès function as For the single-channel case, the Omnès function can be expressed in the analytic form in terms of the phase shift δ(s), with the convention that δ(s th ) = 0. Therefore, in singlechannel approximations, the Omnès function is frequently computed directly from the existing parametrizations of the phaseshift data and various assumptions about its asymptotic behavior at infinity. The latter constrains the asymptotic behavior of the Omnès function: for δ(∞) → απ one obtains Ω(∞) → 1/s α . In our approach, the phase shift curves are obtained from fits to the data using the N/D method. The high-energy asymptotic of the phase shift is coming from the approximation of the lefthand cut by conformal expansion and subsequent solution of the once-subtracted dispersion relation. As a result, in this scheme, the obtained Omnès function (or its inverse) is always asymptotically bounded, if there is no bound state or CDD pole in the system. When there is a bound state in the system, the relation between the Omnès function and the D(s) function given in Eq. (19) changes, where the extra factor (s − s B )/(s M − s B ) removes the zero of D(s). Due to this extra factor, the obtained Omnès function grows linearly at infinity and satisfies the twice-subtracted version of the dispersion relation given in Eq. (17). This can also be seen from the Levinson's theorem, which relates the contribution from the number of bound states n B to the phase shift at infinity as δ(∞) → −n B π (using the convention δ(s th ) = 0). For the multi-channel case, the Muskhelishvili-Omnès equations (17) do not have analytic solutions [46,47], and one needs to find a numerical solution, by employing for instance a Gauss-Legendre procedure [47]. In order to achieve that, however, one needs to know the off-diagonal scattering amplitude in the unphysical region and again make the assumption about the highenergy asymptotics. On the other side, with the N/D method, both the scattering amplitude and the Omnès function are obtained simultaneously from the fit to the available data. Additional information about the off-diagonal scattering amplitude in the unphysical region can be used as a constraint and not as a necessary requirement to obtain the Omnès matrix. Also, as discussed above, in most of the cases the obtained Omnès function (or its inverse) is asymptotically bounded. Therefore, this approach is useful in many practical applications.
As a check of our numerical calculations, we verified that the Omnès functions obtained using Eqs. (7) and (18) satisfy Eq. (17).

Numerical results
In this paper, we study the resonant ππ and πK scattering in the S-wave. These are the channels where σ/ f 0 (500), f 0 (980), and κ/K * 0 (700) resonances reside. Both ππ and πK channels have been measured experimentally [49,50,51,52,53]. However, throughout the whole energy range there are large differences between different data-sets and a careful choice of the data is required to achieve a controllable data-driven description of the phase shifts and inelasticity. For the ππ scattering, the situation is a bit better than for πK scattering, since there is very precise low-energy data coming from K l4 decays [51] and, in general, SU(2) χPT is a much more accurate theory than the √ s E , MeV  (9) 0.220 (5) 0.276 (6) 90 (9) 0.220 (5)   SU(3) version of it. In order to be consistent with χPT in the threshold region, we employ the effective range expansion where a is the scattering length and b is the slope parameter. For the ππ and πK scattering both a and b have been calculated at NNLO in χPT [9,48]. As expected, for the πK scattering, the chiral convergence is a bit worse than for the ππ scattering [48], however the results for the scattering length and slope parameter do not show large discrepancies with the Roy-Steiner results [10,11]. As for the Adler zero, we have checked that its position does not acquire large higher order corrections, and for simplicity one can take the LO result. In all numerical fits, however, we take the NLO result [54] as a central value, with the uncertainties from the omitted higher orders as |NLO − LO|, which should provide a conservative estimate. The NLO values for the low-energy constants are taken from [55]. For the case of non-physical pion masses with m π = 236 MeV and m π = 239 MeV, we only use Adler zero positions as a constraint, while for m π = 391 MeV, where σ/ f 0 (500) shows up as a bound state, no constraints are imposed. The free parameters in our approach are the conformal coefficients in (13), which determine the form of the left-hand cut contribution U ab (s) in Eq. (5). Apart from the standard χ 2 criteria, the number of parameters is chosen in a way to ensure that the series (13) converges. The uncertainties are propagated using a bootstrap approach. In several cases, however, we will be fitting Roy (Roy-Steiner) solutions, which are smooth functions and their errors are fully correlated from one point to another. In these cases, χ 2 /d.o. f loses its statistical meaning and can be < 1. In our fits, this scenario will simply indicate that we obtained the N/D solution which is consistent with the Roy (Roy Steiner) solutions, and we just make sure that the obtained uncertainty is consistent with that from Roy analyses.
Before entering the discussion of the results of the fits, we would like to briefly comment on the freedom of the choice of the subtraction point s M in the dispersion relation (4). The common choice in the application of the Omnès functions is s M = 0, due to its relation to scalar form factors and matching to χPT. On the other side, one can fix s M at the threshold, s M = s th , and then relate n max n=0 C n ξ n (s th ) to the scattering length. Similarly, one can fix s M at the Adler zero 3 , s M = s A , which would imply that n max n=0 C n ξ n (s A ) = 0. The last two choices can therefore reduce the number of fitted parameters by one. Eventually different choices of s M redefine the fitted coefficients C n in the U ab (s) function and the results of the N/D method are immune to that (after computing the D-function, it can be re-normalized to any other point below threshold). Since not in all the fits we impose threshold or Adler zero constraints, we decided to make Our results

Roy-like analyses
Exp., SC 458 (7) Table 3: Poles and couplings of the σ/ f 0 (500), f 0 (980), and κ/K * 0 (700) resonances calculated in data-driven N/D approach compared with the results of Roy-like analyses. SC and CC stand for single-channel and coupled-channel analyses, respectively. For the f 0 (980) or κ/K * 0 (700) poles we take a conservative dispersive average between [7] and [56] or [10] and [11], similar as it was done for σ/ f 0 (500) in [6]. In our results, the first error is the statistical one, while the second one comes from a variation of s E and has a systematic nature. the choice s M = 0 (22) in all the cases for simplicity. As for the expansion point s E , we choose it in the middle between the threshold and the energy of the last data point that is fitted, Note, that in the coupled channel case, s th in Eq. (23) denotes the physical threshold for the diagonal terms U ab (s), while for the off-diagonal terms it is the lowest threshold. We emphasize, that this particular choice guarantees a fast convergence of the conformal expansion (13) in the region where the scattering amplitude is fitted to the data and also where it is needed as input to Eq. (9). Unlike the physical region, where the reaction models are typically fitted to data, the pole extraction may carry significant systematic uncertainties, especially if the pole lies deep in the complex plane [59,60]. To assess these, we vary the parameter s E around its central value fixed to (23). We allow for a conservative variation by 25% of the difference √ s max − √ s th , in order to have a compromise between √ s th and √ s max . Note, that the extreme choice of 50% would correspond to s E = s th or s E = s max , which we clearly want to avoid, since it would bias the fit towards one or the other region. In the following results, the first error will indicate the statistical uncertainty (i.e. reflect the errors of the data and χPT input), while the second one will be associated with a variation of s E . All results presented below have been checked to fulfill the partial-wave dispersion relation given in Eq. (5) or Eq. (12) in the case when there is a physical bound state in the system. In addition we checked that there are no spurious poles or bound states in the considered cases 4 .
3.1. Single channel ππ → ππ analysis of the experimental and lattice data As a first step, we consider only the elastic ππ scattering, which should be enough to get a realistic estimate of the resonance position of σ/ f 0 (500), which is known to be connected almost exclusively to the pion sector. The reason for that is twofold. In many practical applications it is convenient to remove the KK (or f 0 (980)) effects, which do not influence much the σ/ f 0 (500) pole parameters, but at the same time require a proper coupled-channel treatment. Additionally, the current lattice QCD result for m π = 236 MeV covers only the elastic region [32]. Therefore, as a necessary prerequisite of a meaningful σ/ f 0 (500) pole extraction for unphysical pion masses, one has to test the N/D formalism first for physical quark mass values, where the position of σ/ f 0 (500) has already been obtained from the sophisticated Roy analyses [6,7,8,9]. The inclusion of the KK channel (or f 0 (980) resonance) will allow for a slightly more precise evaluation of σ/ f 0 (500) parameters and will be given in the next subsection.
Relying only on the available data up to where a strong influence of the KK threshold is not yet expected, we obtain a decent fit even without imposing chiral constraints. The pole occurs at √ s σ = 463(8) +6 −7 − i 217(6) +8 −9 MeV. The scattering length and slope parameters turn out to be compatible with those of χPT due to the presence of K l4 data. As we discussed above, this is not the case for the Adler zero, which is located too close to the left hand cut, i.e. where the series (13) simply converges too slow. With the additional constraints for the scattering length, slope parameter and Adler zero, the best fit result contains four parameters and leads to √ s σ = 435(7) +6 −8 − i 250(5) +6 −8 MeV. This result is compatible with the value √ s σ = 446(5) +6 −9 − i 230(5) +7 −9 MeV, obtained by replacing the experimental data with the pseudo data from the Roy-like analysis [7]. As it is shown in Fig. 2 both N/D fits are consistent within the error. This provides a proof for our expectation, that even in the case where there is no available Roy analyses (like lattice QCD data), we can rely on the N/D approximation. For our final result of the singlechannel Omnès function with physical pion mass, we opt for fitting the result of the Roy analysis [7], as the best representation of the data. The values of the fitted parameters are collected in Table 1, which result in the fast convergence of the conformal expansion (13) as shown in the left panel of Fig. 2. Note, that in order to use these fit parameters as the starting values of the more complicated coupled-channel fit, we have chosen s E here to be the same as for the coupled-channel case, where we aim to describe the data up to √ s max = 1.2 GeV. Also, this choice slightly improves the obtained σ/ f 0 (500) pole positions, since it pushes s E further away from the threshold region, which is constrained accurately from χPT. In Table 2 we compare threshold parameters and Adler zeros to χPT values, while in Table 3 poles and couplings are collected. Overall we achieve a good description of the Roy analyses results. In Fig. 2 we also show phase shift and Omnès function. Note, that a similar result for the Omnès function can be obtained by using the phase shift from the single-channel modified Inverse Amplitude Method (mIAM) [61,62,63,58] and Eq. (19). In this method, the dispersion relation is written for the inverse amplitude, while the left-hand cut and subtraction constants are approximated by the chiral expansion. The result closest to the Roy analysis for the σ/ f 0 (500) pole is achieved by performing two-loop mIAM fit [64]. In elastic N/D and mIAM approaches the KK channel is separated naturally from the ππ channel, which is beneficial for the practical applications. Apart from the experimental data, the recent lattice analysis [32] provided the results for the energy levels for pion mass values of m π = 236 MeV and m π = 391 MeV. While the former case is much closer to the physical pion mass, the lattice result for the larger mass deserves special attention, since in that case σ/ f 0 (500) shows up as a bound state. In the lattice QCD analysis, the discrete energy spectrum in a finite volume is related to the infinite-volume scattering amplitude through the Lüscher formalism [65], which was extended in [66] to the case of moving frames. In the case of elastic scattering at low energies it gives a one-to-one relation to p cot δ. The lattice results for p cot δ with m π = 236 MeV and m π = 391 MeV were shown in [32]. To fit these data, we analytically continue p cot δ below threshold, such that it does not produce any cusp behaviour at the threshold, where ρ 0 is the same as ρ in Eq. (3), but without the Heaviside step function. For both m π = 236 MeV and m π = 391 MeV, we find that the three-parameter fit covers the data quite well (see central and bottom panels of Fig.2). Similar to the K-matrix fits performed in [32], we found σ/ f 0 (500) as a deep pole on the second Riemann sheet for m π = 236 MeV and as a bound state for m π = 391 MeV. In our approach, however, the obtained scattering amplitudes satisfy p.w. dispersion relations, which is a stringent constraint on the real part of the inverse of the amplitude. As a result, the pole position is determined much more precisely, see Table 3. We also checked that the obtained scattering length m π a = 0.98 (19) for m π = 236 MeV is consistent with the chiral extrapolation result m π a NNLO = 0.75 − 0.87 of [9] and therefore including such additional constraint in the fit barely affects the results of the σ/ f 0 (500) pole and coupling.
It is instructive to compare the obtained pole positions of σ/ f 0 (500) for non-physical pion masses with the predictions of unitarized chiral perturbation theory (UχPT). The most popular are two approaches: mIAM [64] and Bethe-Salpeter equation (BSE) [67]. Both observe the same qualitative behaviour of the σ/ f 0 (500) pole. With increasing pion mass values the imaginary part of the pole decreases, then σ/ f 0 (500) becomes a virtual bound state and as m π increases further, one of the virtual states moves towards threshold and jumps onto the first Riemann sheet and become a real bound state. For m π = 236 MeV, the extracted value from lattice data is consistent with UχPT predictions for the real part, but somewhat lower for the width,    [50,51] (dashed curve) and fit to the pseudo data from Roy analysis [6,7] (thick curve). Note, that for the sake of comparison with the coupled-channel case (see Fig.4), we adopted for this case s E based on √ s max = 1.2 GeV, as discussed in the text.
the results of UχPT are very sensitive to the chiral order. Both mIAM [62] and BSE [67] at one loop found σ/ f 0 (500) as a virtual bound state for m π = 391 MeV. However, including the higher-order corrections (two loop) in mIAM [64] predicted the conventional bound state very close to the lattice results confirming the proposed trajectory. However, as pointed out in [32], it would be useful to perform lattice calculation between 236 and 391 MeV, to see what really happens in the transition region between a resonance lying deep in the second Riemann sheet and the bound state.

Coupled-channel {ππ, KK} analysis of the experimental data
While the single-channel analysis allows us to reproduce the low-energy behavior of the phase shifts and gives very reasonable values of the σ/ f 0 (500) pole parameters, a comprehensive study of the region up to √ s = 1.2 GeV should account for the interplay between ππ and KK channels. In our normalization (see Eqs. (1-3)), the two-dimensional t-matrix, with channels denoted by 1 = ππ and 2 = KK, is given by   Under assumption of two-channel unitarity, the inelasticity is related to |t 12 (s)| as and due to Watson's theorem, In the physical region the t-matrix is fully described by experimental information on the ππ phase shift δ 1 (s) [49,50,51], the inelasticity η(s) (or |t 12 (s)| for s > 4m 2 K [68,69,70]) and the ππ → KK phase δ 12 (s) [69,68,71].
Similar to the single-channel analysis, we first fit the available experimental data supplemented with constraints for scattering length, slope parameter and Adler zero from χPT in the ππ → ππ channel. As for the ππ → KK channel, the complication stems from two facts. Firstly, the experimental data exist only in the physical region above KK threshold. Therefore, in order to stabilize the fits, we make sure that the obtained |t 12 (s)| stays small around 5 s = 0 as a manifestation of χPT. Secondly, the existing experimental data for both |t 12 (s)| and δ 12 (s) contains incompatible data sets and require to make some choice. Since the phase δ 12 (s) is fully defined below KK threshold by means of Watson's theorem, we discard the data from [69] as it suggests that ππ → KK phase goes much lower than it is forced by the presence of f 0 (980) resonance. Therefore, we fit the data from [68] and [71] which are consistent due to the large error bars of the latter set. As for |t 12 (s)|, the two data sets from [70] and [68,69] should in principle be treated separately. However, only the data from [70] is compatible with the ππ inelasticity around the KK threshold. In order to describe the data from [68,69], most likely one has to include the four-pion channel, which is beyond the scope of the present paper. The best fit with (4, 4, 3) parameters in (11,12,22) channels [70], provides σ/ f 0 (500) and f 0 (980) poles at √ s σ = 454(12) +6 −7 − 262(12) +8 −12 i MeV and √ s f 0 = 990(7) +2 −4 −17(7) +4 −1 i MeV. These results are remarkably close to the Roy (for ππ → ππ) and Roy-Steiner solutions for (ππ → KK) as shown in Figs. 3 and 4. The large error bars arise from scarce experimental data around KK threshold and almost unconstrained |t 12 | in the unphysical region.
On the other side, we have at our disposal very precise ππ → ππ Roy-like analyses from [7] and ππ → KK Roy-Steiner analyses from [10,11,30]. Unfortunately, they do not come from the coupled-channel Roy-Steiner analyses and may display some inconsistencies between each other. In particularly, the Roy results on the real and imaginary parts of the t 11 (s) amplitude can constrain δ 1 (s) and η(s). The latter, in the two-channel approximation, is related to |t 12 (s)| by Eq. (29) and turns out to be inconsistent with any available Roy-Steiner solution on ππ → KK [10,11,30]. Therefore in order to avoid possible conflict in fitting two independent analyses, we impose ππ → KK Roy-Steiner solution only as constraint on |t 12 (s)| in the unphysical region 4m 2 π < s < 4m K . Currently, there are three competing solutions: one from Büttiker et al. [10] and two (CFDc and CFDb) from Peláez et al. [11]. We let the fit decide which solution to choose. As for the δ 12 , we take advantage of experimental data of Cohen et al. [68] in the fit, which are quite precise. The good description of the data can be achieved with as low as (4, 2, 3) parameters in (11,12,22) channels, respectively. The results of the fit are collected in Tables 1,2 and 3 and shown in Fig. 4. As expected, the values for the fit parameters in the 11-channel do not deviate much from the single-channel analysis in Sec.3.1. In the coupled-channel analysis the σ/ f 0 (500) pole position comes a bit closer to the Roy analysis value, than in the single-channel study. Moreover, we are now in a position to calculate its coupling to the KK channel, which we include in Table 3. By inspecting Table 1, one can also see the striking similarity between the fit parameters in the 22 channel and the fit to lattice ππ → ππ data with m π = 391 MeV, for which there is a bound state. Similarly, f 0 (980) will be a bound state in the 22 channel, if we eliminate its connection to the 11 channel, i.e. by putting U 12 = 0. This feature is not new and has already been observed in UχPT calculations, see for instance [72]. As for the 12 channel, the fit clearly favours CFDc solution of [11]. This is also consistent with our previous "free" fit to to the experimental data, as shown by the dashed curves in Fig. 4. On the right panels of Fig. 4 we show the elements of the Omnès matrix calculated using Eq. (18). The previous version of them, with the fit to [10,30] has already been successfully applied for the dispersive coupled-channel study of γ ( * ) γ * → ππ(KK) [73] and e + e − → J/ψππ(KK) [28].
We leave the coupled-channel study of the existing lattice data on {ππ, KK} [78] with m π = 391 MeV for a future work. In our opinion, this channel has to be analysed together with {πη, KK} lattice data [79], to shed more light onto the differences between the light scalar resonances f 0 (980) and a 0 (980).
3.3. Two-photon couplings of σ/ f 0 (500) and f 0 (980) As an application of the obtained Omnès functions in the N/D approach, we would like to extract the two-photon cou-  plings of σ/ f 0 (500) and f 0 (980). In principle, the coupling to the external currents has the potential to infer the scalar meson composition. Furthermore, it characterizes the interaction strength of σ/ f 0 (500) and f 0 (980) in the two-photon channel. The latter is important for the light-by-light sum rule applications [80,81,82,83] and serves as a key input to estimate the isoscalar two-pion (kaon) contribution to the hadronic light-bylight scattering for (g − 2) of the muon [84]. The central result in this section will be obtained using a coupled-channel dispersive representation, however, for σ/ f 0 (500) we will employ as well the single-channel representation both for physical and non-physical pion masses.
The photon-fusion partial-wave amplitude γγ → ππ, which we denote by h (J) I,λ 1 ,λ 2 , is the off-diagonal element of the γγ, ππ, KK channels. Since the intermediate states with two photons are proportional to e 4 , they are suppressed, and one can reduce the (3 × 3) matrix dispersion relation down to the (2 × 1) form, which require the hadronic rescattering part, Ω(s), and the left-hand cuts as input [23,24,85,25,86,75]. For the lowenergies around σ/ f 0 (500) (and to lesser extent around f 0 (980)) the contribution from the left-hand cuts is dominated by the pion-pole contribution (Born term), which is exactly calculable. Therefore, in this approximation there is no need of modeling left-hand cuts in one way or another or introducing any subtractions. The photon-fusion p.w. amplitudes are readily obtained using the Muskhelishvili-Omnès representation. For more details, we refer to Ref. [73]. The two-photon couplings are extracted by calculating the residue of h (0) 0,++ (s) at the pole positions, s p . Following [87,86], in our convention it is given by   [74,75,24,56,76]. Right panel: total cross section of γγ → π 0 π 0 (| cos θ| < 0.8) from [73] with updated I = 0, J = 0 contribution. The data for the cross section are taken from [77].
where h (0) 0,++ (s) is evaluated on the first Riemann sheet for p = σ/ f 0 (500), f 0 (980). An intuitive way of re-expressing the twophoton couplings, shown in Table 3, is by using the formal definition of the corresponding two-photon decay widths Converted to (32), our results read where in square brackets the single-channel approximation is shown. As expected, its Γ σ→γγ is almost indistinguishable from the coupled-channel case. In Fig. 5 we compare our results with the recent dispersive estimates [74,75,24,56,76]. While the two-photon decay width of f 0 (980) is consistent with the coupled-channel amplitude analysis of [76] and the over-subtracted coupled-channel Muskhelishvili-Omnès analysis [56], the two-photon width of σ/ f 0 (500) is about 25% smaller than their values. On the other hand, the obtained twophoton width of σ/ f 0 (500) is consistent with the sophisticated Roy-Steiner analysis [24] and other dispersive analyses from [74,75]. Finally, we also predicted σ/ f 0 (500) two-photon coupling/width for the unphysical m π = 236 MeV, which would be interesting to confront with the direct lattice calculations. We note, that the errors quoted in Eq. (33) correspond solely to the uncertainties in the Omnès matrix. In principle, one can perform a more comprehensive study of the theoretical uncertainties, by the inclusion of more distant left-hand cuts in γγ → ππ(KK). This would require introducing subtraction constants which can be either fixed from the pion dipole polarizabilities or fitted directly to the cross-section data. Doing so would likely enlarge the error, but we do not expect a significant change of the central values, since the current parameter-free description of the cross-section data (see Fig.5) is quite impressive. The advantage of the approach that accounts only for pion pole left-hand contribution, is that in the absence of any singlevirtual data one can predict the behavior of the p.w. helicity amplitudes for finite virtualities [88,73], which are needed as input for (g − 2) µ [84].

I = 1/2 single-channel: data and lattice
For the πK → πK single channel analysis we begin by fitting the experimental data and imposing constraints from χPT for the scattering length, slope parameter, and Adler zero. The latter at LO is given by a simple relation, The most precise calculation of the scattering length and slope parameter in χPT has been performed at NNLO in [48]. While the result for the scattering length m π a = 0.22 is consistent with the recent Roy-Steiner predictions m π a = 0.223(9) [11], it seems that there is a small tension in the slope parameter value m 3 π b = 0.13 compared to m 3 π b = 0.108(8) from [11]. The calculation of uncertainties is a bit cumbersome at NNLO and has not been presented in [48]. Therefore in our fits we take NNLO χPT values as central results, but include the conservative error-bar, such that it covers the recent Roy-Steiner results [11]. As for the Adler zero, we take the NLO value, as explained at the beginning of Sec.3. The available experimental data for this process is scarce in the region close to the πK threshold, and often contains the discrepancies even within one dataset [52]. Since we consider only the single-channel approximation, we perform the fit till ηK threshold of the data from [52,53]. In this way we also exclude the influence of the K * 0 (1430) resonance. We observe a similar situation as for the ππ → ππ single-channel analysis, that fitting the experimental data [52,53] or Roy-Steiner analysis of [11] provides  MeV, respectively. In general, these results compare well with the Roy-Steiner pole position 653 +18 −12 − i 280(16) MeV which we take as a conservative average between [10] and [11]. The one-sigma difference in the resonance mass can be attributed to the fact, that we are fitting Roy-Steiner solution only in the elastic region. We also look forward to the results of the KLF Collaboration, which plans to study πK scattering using a secondary K L beam at Jefferson Lab [89]. It will further improve the position of the κ/K * 0 (700) resonance. For the unphysical pion mass, we again use recent lattice data from the Hadron Spectrum Collaboration [36]. We analyse the data for m π = 239 MeV, where an evidence of κ/K * 0 (700) was observed in the p cot δ distribution. Due to large uncertainties, the pole position was not determined by the lattice collaboration, calling for more sophisticated approaches that include in addition to unitarity also the analyticity constraint. By employing the data-driven N/D approach, the present data can be easily described with the two-parameter fit, leading to χ 2 /d.o. f = 0.4. In this case, however, the Adler zero of the amplitude is located relatively far from the χPT value, since the lattice data in the low p 2 region suffers from the large uncertainties. Also, as we discussed before, in the Adler zero region the conformal expansion (13) does not converge well by construction and one has to impose Adler zero as a constraint, which effectively calls for one additional parameter. In this way, the impact of two points with prominently small errors at p 2 ∼ 0.09 and ∼ 0.11 GeV 2 is balanced out. The results of the fit are collected in Tables 1,2  and 3.
Again we would like to compare our results for the pole position and coupling with predictions of mIAM. According to [53], at m π = 239 MeV, the imaginary part of the pole decreases by ∼ 17%, while the real part and coupling slowly increase by ∼ 4% and ∼ 8% respectively. Our values extracted from the lattice data show a similar behavior, with the decrease in the imaginary part of 7.0(7.7)%, increase in the real part and coupling of 6.4(5.8)% and 1.6(6.4)%, respectively.

Systematic uncertainties
In the end, we wish to comment on the size of systematic uncertainties of our results. As it can be seen in Table 3, under the change of √ s E by 25% of the difference √ s max − √ s th around the central value (23), the σ/ f 0 (500) and κ/K * 0 (700) poles acquire noticeable systematic errors which are of the size of statistical ones. We admit that s E variation only accounts for the dominant part of the systematic uncertainty and therefore only provides a lower bound on the systematic error. However, even if we go to the extreme case of 50%, which corresponds either to s E = s th or s E = s max , the statistical error will grow only by a factor of two, compared to the case of 25%. This is different from the K-matrix fits (see for instance [59]), which cannot extract accurately the pole parameters. We remind, that in our approach, as opposed to K-matrix models, the obtained amplitudes satisfy p.w. dispersion relations, which is an additional constraint on the amplitude both on the real axis and in the complex plane.

Conclusion and outlook
In this work, we presented a data-driven analysis of the resonant S-wave ππ → ππ and πK → πK reactions using the p.w. dispersion relation. In this approach unitarity and analyticity constraints are implemented exactly. We accounted for the contributions from the left-hand cuts using the Taylor expansion in a conformal variable, which maps the left-hand cut plane onto the unit circle. Then, the once subtracted p.w. dispersion relation was solved numerically by means of the N/D method.
Using existing experimental information and threshold constraints from χPT we tested the single-channel N/D formalism for the physical pion mass, where the positions of σ/ f 0 (500) and κ/K * 0 (700) have already been obtained from the sophisticated Roy and Roy-Steiner analyses. We demonstrated that the results for the pole parameters are stable and almost do not change if we replace the existing experimental data with the very precise pseudo data generated by Roy and Roy-Steiner solutions in the physical region. As a next step, we performed the fits to the lattice data of the Hadron Spectrum Collaboration for m π = 236, 391 MeV in the case of ππ → ππ and for m π = 239 MeV in the case of πK → πK. We provided an improved determination of the σ/ f 0 (500) and κ/K * 0 (700) pole parameters compared to the simplistic K-matrix approach and also compared them with UχPT predictions.
An important feature of the N/D method is that the Omnès function comes out naturally, as the inverse of the D-function. The knowledge of the Omnès function, in turn, allows employing the Muskhelishvili-Omnès representation for the vast majority of production/decay reactions involving two pions (or pion and kaon) in the final state. While for the single-channel case, the Omnès function can be obtained analytically from the parametrisation of the phase shift, this is not the case for the coupled-channel case. In order to cover the f 0 (980) region we extended our analysis for the coupled-channel {ππ, KK} case and extracted the corresponding Omnès matrix. In our construction it is asymptotically bounded (i.e. it satisfies oncesubtracted dispersion relation) and therefore useful in many dispersive applications. The unknown coefficients from the conformal expansion were adjusted to reproduce existing Roy and Roy-Steiner analyses. As a straightforward application of the Muskhelishvili-Omnès representation, we estimated the twophoton decay widths of the σ/ f 0 (500) and f 0 (980) resonances, which turned out to be consistent with the previous dispersive results. The obtained Omnès matrix serves as an important building block, which allows for the dispersive calculation of the isoscalar two pion/kaon contribution to the hadronic lightby-light part [90,88,91] of the anomalous magnetic moment of the muon (g − 2) µ [84]. In particularly, with the input from γ * γ * → ππ, KK [73] one can estimate dispersively the contribution from the f 0 (980) resonance, and compare it with narrow resonance results [82].
The proposed method is not only limited to the ππ and πK scattering. We considered these reactions in the present paper because they show up as building blocks in many hadronic reactions/decays and have been calculated recently using lattice QCD. In principle, the N/D method combined with the conformal expansion for the left-hand cuts can be applied to any hadronic reaction where there is data (experimental or lattice) which possesses a broad (or coupled-channel) resonance that does not have a genuine QCD nature. For the latter (like for instance ρ or K * resonances) one needs to extend the formalism to allow for CDD poles. Also, it has to be modified in the presence of anomalous thresholds.