On the amplitudes for the CP-conserving K±(KS) → π±(π0)ℓ+ℓ− rare decay modes

The amplitudes for the rare decay modes K± → π±ℓ+ℓ− and KS → π0ℓ+ℓ− are studied with the aim of obtaining predictions for them, such as to enable the possibility to search for violations of lepton-flavour universality in the kaon sector. The issue is first addressed from the perspective of the low-energy expansion, and a two-loop representation of the corresponding form factors is constructed, leaving as unknown quantities their values and slopes at vanishing momentum transfer. In a second step a phenomenological determination of the latter is proposed. It consists of the contribution of the resonant two-pion state in the P wave, and of the leading short-distance contribution determined by the operator-product expansion. The interpolation between the two energy regimes is described by an infinite tower of zero-width resonances matching the QCD short-distance behaviour. Finally, perspectives for future improvements in the theoretical understanding of these amplitudes are discussed.


Introduction
Rare kaon decays have been playing a crucial role in flavour physics, and more generally in particle physics. Since they proceed through strangeness-changing neutral currents, they are quite suppressed in the standard model (SM), and thus offer a window through which we might possibly catch a glimpse of new physics (NP). This research subject has been thoroughly reviewed in ref. [1], where a detailed list of references can be found. It is again becoming quite timely nowadays, thanks to experiments like NA62 at the CERN SPS, planning to collect, already in 2018, the data required in order to measure the decay rate JHEP02(2019)049 of K + → π + νν with a precision of 10 % of its SM prediction [2,3], and KOTO at J-PARC, aiming at measuring the decay rate of K L → π 0 νν, at around the SM sensitivity in a first stage [4][5][6].
The rare processes K ± → π ± + − and K S → π 0 + − are analogous to those mentioned in eq. (1.1), and appear thus as particularly suitable in order to uncover possible violations of lepton-flavour universality (LFUV) in the kaon sector [44]. The experimental programs of NA62 [45] (for charged kaon decays) and of LHCb [46,47] (for K S decays) offer quite interesting prospects in this regard. Trying, in parallel, to improve our theoretical understanding of these processes therefore constitutes a quite timely undertaking. It is thus not surprising that efforts in this directions have also become part of the agenda of the lattice-QCD community [48][49][50].
Taking only the high-precision data collected by the NA48/2 Collaboration [56,57] into account, one obtains instead the almost two times more accurate result R K ± [NA48/2] = 0.309 (43). In the neutral-kaon sector, the observed decay rates are [59,60] Br[K S → π 0 e + e − ] mee>0.165 GeV = [3.0 +1.5 −1.2 (stat) ± 0.2(syst)] × 10 −9 , Br[K S → π 0 µ + µ − ] = [2.9 +1.5 −1.2 (stat) ± 0.2(syst)] × 10 −9 . (1.4) These measurements have not yet reached the level of precision already available in the case of the charged kaon. At the theoretical level, the situation for the CP-conserving decays K ± (K S ) → π ± (π 0 ) + − is less favourable than for the K → πνν modes. Indeed, whereas the latter are dominated by short distances, the former are governed by the long-distance process K → πγ * → π + − [61][62][63], involving a weak form factor, W + (z) in the case of K ± → π ± γ * , and W S (z) in the case of K S → π 0 γ * . These form factors are given by the matrix elements of the electromagnetic current between a kaon state and a pion state, in the presence of the weak interactions, considered at first order in the Fermi constant. The corresponding momentum transfer squared, the di-lepton invariant mass squared s = zM 2 K , is small enough, it ranges from 4m 2 to (M K − M π ) 2 , where m is the mass of the charged lepton, so that these processes can be treated within the framework of chiral perturbation theory (ChPT) [64][65][66][67]. Conservation of the electromagnetic current implies a vanishing contribution at lowest order, O(E 2 ) in the chiral expansion. A non-vanishing form factor is generated at next-to-leading order, both by pion and kaon loops, as well as by order O(E 4 ) counterterms. The main feature of this one-loop representation of the form factor is the appearance of a unitarity cut on the positive real-s axis, due to the two-pion intermediate state (there is also a cut due to the opening of the KK channel, but the later is located far enough from the kinematic region of interest, so that its contribution can be, for all practical purposes, approximated by a polynomial). The corresponding absorptive part (discontinuity) is given by the product of the pion electromagnetic form factor F π V (s) with

JHEP02(2019)049
the P -wave projection f Kπ→π + π − 1 (s) of the weak Kπ → π + π − amplitude ((K, π) stands for either (K ± , π ∓ ) or (K S , π 0 )), both taken at tree level. These basic properties of the one-loop form factor already suggest some simple ways to improve upon this result and to collect some O(E 6 ) effects. Specifically, in ref. [63], the O(E 2 ) Kπ → π + π − vertex was replaced by the phenomenologically well-known Dalitz-plot expansion of the K → ππ + π − amplitude, and a slope was added to the pion form factor. Furthermore, slopes b +,S in the di-lepton invariant mass squared s, required by the data, but also by order O(E 6 ) counterterm contributions, were included in the weak form factors W +,S (z), in addition to the constant terms a +,S ∝ W +,S (0), already generated by the O(E 4 ) counterterms. There are, however, several reasons to go beyond the approximations considered by the authors of ref. [63]: • The parameters a +,S and b +,S appear merely as phenomenological constants, which have to be fixed from the available data. We will provide an update of the situation using the most recent data on the decay distributions. At this stage, one may already observe that the value (1.3) obtained for R K ± is rather different from the one quoted in ref. [63], R K ± = 0.167 (36), and based on the older data from refs. [52] and [51]. This quite substantial increase in the value of R K ± confirms the prediction R K ± > ∼ 0.23 made by the authors of ref. [63]. We will also test for potential effects of the data on the value of the curvature (quadratic slope) of the K ± → π ± π + π − Dalitz plots.
• Final-state rescattering effects are only partially taken into account by the parameterization of the form factors W +,S (z) proposed in ref. [63]. In particular, the one-loop P -wave projections of the weak Kπ → π + π − amplitudes also involve pion rescattering in the crossed channels. Since the two pions are in the P wave, these are not expected to be as large as, for instance, in the case of the K → ππ amplitudes [1], where the pions are in the S wave. Nevertheless, they could have an impact on the shape of the form factor, and hence on the determinations of b + and b S .
• The phenomenological values of a + and b + are comparable in size, whereas, from naive chiral counting, one would expect b + to be substantially smaller than a + . A theoretical explanation of |b + /a + | ∼ 1 is still lacking.
• Let us close this list with the most important point. The potential detection of any manifestation of LFUV hinges on the possibility to obtain independent information on a + and b + , such as to be able to make a prediction for the ratio R K ± . This definitely requires to go beyond the low-energy expansion itself.
In order to address these issues, we improve the existing theoretical description of the form factors W + (z) and W S (z) in several ways: • We provide complete order O(E 6 ) representations of these form factors. They are obtained as the result of a two-step recursive procedure. The first step reproduces the one-loop representations of ref. [61]. The second step includes, besides the effects already accounted for by the representation of ref. [63], all one-loop pion-pion rescattering effects, both in the pion form factor and in the K → 3π vertex.

JHEP02(2019)049
• We extend the previous representation of W + (z) beyond the low-energy domain upon using a simple parameterization of the pion form factor that describes the data over a large energy range including the region of the ρ(770) resonance. Since no data are available as far as the K + π − → π + π − scattering amplitude is concerned (existing data only concern the decay region, in the form of Dalitz-plot expansions), we proceed by applying a unitarization procedure to the O(E 4 ) P -wave projection of this amplitude. These two items allow us to construct the absorptive part of the form factor W + (z) arising from two-pion intermediate states. The form factor W + (z) itself is then recovered through an unsubtracted dispersion relation. Upon comparing the behaviour of this "exact" form factor at small values of z with its low-energy expansion, we obtain sum rules for the parameters a + and b + , which we can then compare to their direct determinations from data.
• We model contributions due to other intermediate states, which occur at higher thresholds, by an infinite sum of single-resonance states, with couplings tuned such as to correctly reproduce the known high-energy behaviour of the form factor in quantum chromo-dynamics (QCD).
For completeness, we should mention that there are other theoretical studies [68][69][70] of the form factors W (z) that go beyond the low-energy expansion. In ref. [68], a two-parameter representation for W + (z) is proposed, combining the lowest-order chiral expression with resonance exchanges. In ref. [69] the form factors are described within the large-N c treatment of weak hadronic matrix elements, see ref. [71] and references quoted therein, matching the quadratic cut-off dependence of the low-energy contribution with the logarithmic one from the short-distance part. Finally, the authors of [70] give a representation of the form factors in terms of meson form factors, but the issue of the matching with the short-distance part is not addressed. For a recent account on the theoretical situation of the K → π + − amplitudes, see ref. [72].
The content of this paper is consequently organized as follows. In section 2, we first discuss general features of the weak form factors W +,S (z) in QCD, and then address more specifically their short-distance behaviour. Long-distance properties of the form factors are described in section 3, where we recall the results of the existing one-loop calculations [61,73], as well as the beyond-one-loop representation of ref. [63]. Phenomenological aspects linked to the determination of a + and b + from the recent high-precision data are the subject of section 4. Starting from a discussion of the analyticity properties of the form factors, we next construct (section 5) a two-loop representation that accounts for all ππ rescattering effects at this order. We compare this representation of the form factors with the one used in ref. [63], and discuss the impact on the determination of the parameters a +,S and b +,S . Section 6 constitutes the main part of the paper as far as the issues raised above are concerned. There, we construct a dispersive model for the form factor W + (z). The corresponding absorptive part includes the contribution, now not restricted to low energies, of the ππ intermediate states, and an infinite sum over zero-width resonances, with couplings chosen such as to provide matching with the short-distance behaviour established previously in section 2. Assuming unsubtracted dispersion relations, we evaluate the pa-JHEP02(2019)049 rameters a + and b + through the corresponding sum rules. A final section summarizes our results and provides our conclusions, as well as critical comments and some perspectives for the future. Numerical values for various input parameters are collected in appendix A. Technical details related to the content of section 5 are displayed in appendix B.

Theory overview
In this section, we provide background material concerning the theoretical aspects of the amplitudes describing the weak Kπ → + − transitions in the standard model, which will also allow us to set up the notation to be used in this study. We first start with the description within the effective four-fermion theory below the electroweak scale, and turn next to the short-distance behaviour of the form factors in three-flavour QCD.

The structure of the form factors in three-flavour QCD
In the standard model, weak non-leptonic ∆S = 1 transitions of hadrons are described, at a low-energy scale, i.e. below the charm threshold, and at first order in the Fermi constant G F , by an effective lagrangian L ∆S=1 (x) given by [74][75][76][77][78][79] This expression involves the current-current four-quark operators Q 1 and Q 2 , as well as the QCD penguin operators Q 3 ,. . . Q 6 . At lowest order in both the fine-structure constant α and the Fermi constant G F , and from a low-energy (long distance) point of view, the two CP-conserving transitions K ± (k) → π ± (p) + (p +) − (p −) and K S → π 0 (p) + (p +) − (p −) proceed through the one-photon exchange process K → πγ * , so that their amplitudes will involve the weak form factors W +,S (z), z = s/M 2 K , s ≡ (k − p) 2 = (p + + p − ) 2 , defined as (we use the notation W (z) to stand for either W + (z) or W S (z), whenever the discussion applies to both channels) (2.2) In terms of these form factors the amplitudes read Here, j ρ (x) stands for the electromagnetic current corresponding to the three lightest quark flavours (e q denotes their charges in unit of the positron charge), Electroweak quark-penguin operators, as well as the mixing of Q 1 ,. . . Q 6 with them, give contributions of order O(α 2 G F ) to the amplitudes, and will not be considered here.

JHEP02(2019)049
The four-quark operators and their Wilson coefficients C I (ν) depend on the QCD renormalization scale ν, and their evolution with respect to this scale is given by the renormalization-group equations The pure QCD anomalous-dimension matrixγ(α s ) for three flavours, whose matrix elements γ I,J (α s ) occur in these equations, is known at leading-order (LO) and next-to-leading (NLO) accuracy,γ 6) and the corresponding coefficients γ I,J can be found in refs. [74][75][76]79] and [80][81][82][83][84], respectively. At this stage, let us make two remarks.
First, we should stress that the form factors W +;S (z; ν) defined by eq. (2.2) depend actually on the QCD renormalization scale ν. Indeed, although the two composite operators involved in this definition are separately finite, the electromagnetic current because it is conserved, the lagrangian L ∆S=1 for strangeness changing non-leptonic transitions by explicit renormalization of the four-quark operators and renormalization group evolution of the Wilson coefficients, as given in eq. (2.5), their time-ordered product is singular at short distances, and needs to be renormalized [85,86]. Indeed, in the short-distance analysis of the K → π + − transitions, two mixed quark-lepton four-fermion operators, having the factorized form of a quark current times a leptonic current, are also encountered, and provide an additional contribution to the effective lagrangian [77,[85][86][87]: The presence of the axial-current operator Q 7A reflects the contribution of the Z 0 and of W -box diagrams, which also contribute to Q 7V . More importantly, however, the operator Q 7V as well receives the contribution from the electromagnetic-penguin type of diagram, with heavy quarks in the loop, as discussed in [85,86]. This contribution will induce, already at order O(α 0 s ), i.e. even before QCD corrections are applied, a dependence on the renormalization scale ν in the Wilson coefficient C 7V , which can be expressed as (2.10)

JHEP02(2019)049
As revealed by their structures, the two operators Q 7V and Q 7A are finite, and do not depend on the renormalization scale ν, as long as only QCD corrections are considered, which will be the case here. At lowest order, the coefficients γ J,7 are given by [85][86][87][88] 6,7 = 0 (2.11) for three active flavours u, d, s. Their values at next-to-leading order are also available from the literature [89]. The Wilson coefficient C 7A does not depend on ν. Moreover, it is proportional to the very small quantity τ ≡ −V td V * ts /V ud V * us , |τ | ∼ 1.6 · 10 −3 , and the contribution of Q 7A can be neglected as long as one does not discuss issues related to the violation of CP or processes that are short-distance dominated. Keeping only the contribution from Q 7V , the amplitudes in eq. (2.3) actually read They involve the form factors f Kπ + (s), defined as (the plus sign applies for (K, π) = (K ± , π ∓ ), and the minus sign for (K, π) = (K S , π 0 ), in agreement with the phase convention chosen in eq. (2.2)) (2.13) and normalized to f Kπ + (0) = 1 in the limit where the up, down and strange quarks have equal masses. The contribution from Q 7A , which we have omitted, would also involve the form factors f Kπ − (s), but multiplied by the lepton mass and the leptonic pseudoscalar density. The amplitudes being observables, they should no longer depend on the scale ν. This means that the scale dependence coming from the short-distance singularity of the form factor W (z; ν) within three-flavour QCD has to cancel the ν-dependence of the Wilson coefficient C 7V (ν) generated at the electroweak scale, (2.14) We will come back to this issue in greater detail in section 2.2 below.
Second, notice that throughout we are working within the framework provided by pure QCD with three flavours of light quarks. Knowledge on the manner how three-flavour QCD is embedded into the full standard model is not required. In particular, the fourquark operators evolve, at all scales, according to the three-flavour matrix of anomalous dimensionsγ(α s ) (truncated, in practice, at NLO). Actually, from this point of view, the only input from the SM which is required, besides, of course, the structure of the effective JHEP02(2019)049 Figure 1. The diagram corresponding to the leading short-distance contribution to the operatorproduct expansion of the time-ordered product in eq. (2.15). The external lines correspond to the insertion of the current j µ (x) (wiggly line) and of the quark fields d ands. The circled vertex materializes the insertion of a four-quark operator Q I , I = 1, . . . , 6. lagrangian itself, as given by eqs. (2.1) and (2.9), are the "initial" values of the Wilson coefficient C I (ν 0 ), I = 1, . . . 6, and C 7V (ν 0 ) at some scale ν 0 slightly above 1 GeV, but in any case below the charm threshold.

Properties of the form factors at high momentum transfer
This subsection is devoted to the discussion of the behaviour of the form factors W +,S (z; ν) at high momentum transfer, z → −∞, within the framework of three-flavour QCD. For this purpose, we consider the short-distance properties (in what follows, q µ is an Euclidean four-vector, q 2 < 0, whose components become all simultaneously large) of the time-ordered product of the electromagnetic current with the various four-quark operators Q I that appear in L ∆S=1 , One can identify several contributions with the appropriate quantum numbers to the corresponding operator-product expansion (OPE). The leading contribution occurs at order O(q 2 ), and consists of the term It is shown in figure 1, and corresponds to a perturbative intermediate state, made up by a light-flavour quark-antiquark pair. In the absence of QCD corrections, only the treelevel four-quark operators are involved. Gluonic corrections also contribute to this class of short-distance behaviour. They will both renormalize the four quark operators and build up the Wilson coefficient for the O(q 2 ) term in the OPE. At order O(q), one encounters several possibilities, e.g. (D τ denotes the QCD covariant derivative) and so on. We will only consider the leading contribution, at order O(q 2 ), and without QCD corrections, although we briefly comment on the latter below.
In the case of the two current-current operators Q 1 and Q 2 , it is then enough to study the leading short-distance behaviour of (i, j, k, l denote colour indices)

JHEP02(2019)049
with T jl 1 ik = δ l i δ j k for Q 1 , and T jl 2 ik = δ j i δ l k for Q 2 . The evaluation of the loop diagram is straightforward. With dimensional regularization, one finds with (ν stands for the renormalization scale in the MS scheme, ν MS ≡ νe γ E /2 / √ 4π, and D = 4 − 2ε) (2.20) The discussion of the penguin operators Q 3 and Q 4 actually turns out to be quite similar to the one for the operators Q 1 and Q 2 . Indeed, from the expressions of these operators, one sees that there are two ways to perform the contraction with the current j µ shown in figure 1. The first one consists in contracting with the quark bilinears occurring in the sum over flavours. From the point of view of the colour structures, this contraction corresponds to T jl 1 ik for Q 3 , and to T jl 2 ik for Q 4 . However, since the quark masses are irrelevant for the leading short-distance behaviour, this leads to the sum of three identical contributions, weighted by the corresponding quark charges. But since e u + e d + e s = 0, this weighted sum vanishes. There only remains to consider the second possibility, where the contraction is done with the d (ors) quark from the term in front of the flavour sum, and with thed (or s) quark from the second (or third) term of this sum. But this amounts to the same computation as before for Q 1 and for Q 2 . The colour structure now corresponds to T jl 2 ik for Q 3 , and to T jl 1 ik for Q 4 . Furthermore, instead of multiplying by e u , one multiplies by e d + e s = −e u . Finally, for the remaining operators Q 5 and Q 6 one obtains a vanishing result. This is due to their (V − A) ⊗ (V + A) structure, which makes the factor arising from the Dirac matrices vanish, whether one uses the naive dimensional regularization (NDR) [90] or the 't Hooft-Veltman (HV) [91,92] scheme. The product of Dirac matrices in eq. (2.19) can also be simplified, but this time the result will depend on which scheme is being used. After minimal subtraction, one obtains where the general form of the Wilson coefficient reads The result of the OPE with the complete lagrangian (2.1) then writes as The preceding short-distance analysis tells us that the OPE of the time-ordered product of the (three-flavour) electromagnetic current with L ∆S=1 is dominated by the same axial-current operatorsγ ρ (1−γ 5 )d that also appears in the expression of Q 7V . Furthermore, it also exhibits a short-distance singularity, which is renormalized by minimal subtraction, leaving over a dependence on the MS renormalization scale ν. At the level of the dimensionally regularized form factor itself, this translates into (2.28) and implies that the scale-dependence of the form factors is given by (2.29) Turning now toward the condition (2.14), we find that it is indeed satisfied at order O(α 0 s ), where it reads

JHEP02(2019)049
since the comparison with eq. (2.11) shows that One may actually turn the preceding argument around, and, starting with the requirement that eq. (2.14) be satisfied, obtain information on the Wilson coefficients ξ I (α s ; ν 2 /s) on the right-hand sides of eqs. (2.27). Indeed, from The combination on the left-hand side will thus be scale independent at NLO provided that in addition to (2.31) the relations hold. Performing the calculation gives, for instance, where the number of active flavours is N f = 3. The coefficients ξ I p0 are not constrained by this type of argument. They can only be determined by an explicit calculation of QCD corrections to the diagram of figure 1, which we will however not attempt to perform here.
Before closing this section, let us briefly leave our three-flavour world, and consider how the preceding discussion is modified in the presence of a fourth quark flavour, which corresponds to the situation considered in lattice calculations [48][49][50]. For m c < ν < m b , the effective lagrangian reads

JHEP02(2019)049
The operators appearing in this expression are ≡ Q 2 and, moreover, that in the QCD penguin operators the sums over the quarks q, as in eq. (2.21), span the whole range of still active flavours, i.e. q = u, d, s, c in the present case. Likewise, the electromagnetic current is the one corresponding to four active flavours. Repeating the same exercise as before, one now obtains (recall that all quarks are considered to be massless) The scale dependence is now proportional to the small quantity τ , defined after eq. (2.11).
To the extent that CP-violating effects are not considered, these contributions can be safely omitted for all practical purposes. But strictly speaking, the absence of a short-distance singularity in the form factor holds only, in the case of four active flavours, within this approximation. This picture is consistent with the fact that, above the charm threshold, the Wilson coefficient C 7V of the operator Q 7V is also proportional to τ [89].

Low-energy expansion of the weak form factors
In this section, we recall the properties of the form factors W +,S (z) from the point of view of their low-energy expansion. We first give their expressions at one loop and discuss some of their properties. We consider next the regime where z 1 (s M 2 K ), so that the pion loops remain as the only sources of non-analyticity. In this regime, we then briefly discuss the description of the form factors proposed in ref. [63].

The form factors at one loop
Since the K ± → π ± γ * and K S → π 0 γ * transition form factors vanish at lowest order [61], the low-energy expansions of the form factors W +,S (z) start at order one loop. These one-loop expressions have been computed in ref. [61] as far as the octet component is concerned. The contribution of the 27-plet has only been worked out more recently, in JHEP02(2019)049 ref. [73]. In a notation slightly different from the one used in these references, the resulting expressions read The functionJ ππ is defined defined as (λ(x, y, z) denotes the Källen function) Explicit expressions for this function are given in refs. [65][66][67]. The first term in eq. (3.1) gathers the contributions from the O(E 4 ) counterterms, which can be expressed in terms of low-energy constants defined in ref. [93] and chiral logarithms (µ χ denotes the renormalization scale in the effective low-energy theory) as and These combinations of counterterms and chiral logarithms do not depend on µ χ : dw (8,27) +,S /dµ χ = 0. The second and third terms in the expression (3.1) come from the pion and kaon loops, respectively. α tree + (α tree S ) corresponds to the P -wave projection of the amplitude of the reaction K + π − → π + π − (K S π 0 → π + π − ) at order O(E 2 ) in the low-energy expansion, and coincides with the linear slope, evaluated at lowest order. of the Dalitz plot of the decay K ± → π ± π + π − (K S → π + π − π 0 ). The interpretation of the quantityα tree + (α tree S ) is similar, but in terms of the amplitude for the reaction In this case there is no decay region, andα tree + and α tree S rather correspond to subtheshold parameters in the expansions of the amplitudes for K + π − → K + K − and K S π 0 → K + K − around the centres of their respective Dalitz plots (s = M 2 K + M 2 π /3, t = u). In terms of the two constants g 8 and g 27 that describe the JHEP02(2019)049 K → ππ amplitudes at lowest order [1], one has (for the numerical values, see table 3 in appendix A) and As already mentioned in the introduction, the physical region of the decays In this range of values the function J KK (zM 2 K ) is well described by its Taylor expansion at z = 0. Performing this expansion in the contribution from the kaon loops, and keeping terms at most linear in z, allows to rewrite the one-loop form factors as Notice that a 1L +,S and b 1L +,S correspond to the values of the form factors and of their slope, respectively, at z = 0, The one-loop expressions of the constants a 1L +,S and b 1L +,S that result from this simple exercise read a 1L +,S = a CT-1L Predicting the values of a 1L +,S from their one-loop expressions requires some knowledge of the O(E 4 ) counterterms, in particular N 14 and N 15 , which occur in the dominant octet part. Several proposal have been made [93][94][95][96][97] in order to extend the determination of the low-energy constants through resonance saturation in the strong sector [98,99] to the weak sector. Unfortunately, these extensions involve unknown parameters (see also the discussion at the beginning of section 6.2 below), thus requiring additional assumptions. This JHEP02(2019)049 unfortunate situation introduces an uncontrolled model dependence in the predictions that can be made within this framework. The prospects for a phenomenological determination of some of these constants from (real or virtual) radiative kaon decays have been discussed in the recent account [100]. In this context, we also recall the "octet dominance hypothesis" discussed in refs. [61,68], which corresponds to the assumption that the (scale independent) . This brief description of the structure of the one-loop expressions of the form factors (we refer the reader to the original articles for further details) brings us to a few remarks: • As already noticed in ref. [61], there are substantial differences in the structures of the two form factors W +;1L (z) and W S;1L (z), and therefore also between a 1L + and a 1L S . In particular, pion loops are suppressed by the ∆I = 1/2 rule in W S;1L (z), but kaon loops are about three times as important as in W +;1L (z) (in absolute value).
• Besides providing tiny slopes b 1L +,S , the kaon loops also contribute to a 1L +,S , see eqs. These contributions are proportional toα tree + andα tree S , respectively. Higher-order corrections will moveα tree +,S to their phenomenological valuesα +,S , which are however not known, see the discussion after eq. (3.5). Assuming, for the sake of illustration, that this change is of about the same size as the change from α tree +,S to α +,S , i.e.α + ∼ 2.5α tree + andα S ∼ 1.3α tree S (see eqs. (3.17) and (3.18) below), we would conclude that the kaon loops could contribute to a + (a S ) at the level of ∼ +0.15 (∼ −0.25). Of course, this argument is at best indicative of the fact that contributions from kaon loops to a +,S , in contrast to b +,S , could be substantial. Besides, higher-order corrections will also produce other effects on the form factors than the simple replacement ofα tree +,S byα +,S , but their impact on a +,S are however more difficult to estimate without an explicit computation.
• It is also possible, and of interest for the sequel, to look at the one-loop form factors in terms of their analyticity properties, which follow from the structure of the Feynman diagrams shown in figure 2. These properties can be expressed in the form of a once-subtracted dispersion relation

JHEP02(2019)049
with the discontinuity along the positive real-z axis provided by one-loop unitarity.
Since at this order the only intermediate states are π + π − and K + K − , see figure 2, it reads ), and are the P -wave projections mentioned after eq. (3.5).
The expressions of Im W +,S;1L (s/M 2 K ) also show that one subtraction is necessary, but sufficient, in order to obtain a convergent dispersive integral. Performing this integration leads to Expanding the contribution from the kaon loop as before allows to cast this expression into the form given in eq. (3.8). Some information is lost within this dispersive approach, namely the way how a 1L +,S , which appears here as a mere subtraction constant, decomposes into its various components, as displayed in eq. (3.10). The issue raised in the preceding item of this list, having to do with local contributions from the kaon loops, can therefore not be addressed within this dispersive approach.
• On the other hand, extending the absorptive parts beyond their lowest-order expressions (3.12) provides a starting point for establishing the structure of the form factors at two loops (beyond two loops, also intermediate states with more than two mesons need to be considered). Furthermore, one may restrict the contributions to the absorptive parts to two-pion states from the outset, and include the other two-mesons states (KK, but also, for instance Kπ) into the subtraction polynomial, which becomes a first-order polynomial in z at two loops. The explicit construction of the two-loop form factors along these lines will be the subject of section 5.

The form factors beyond one loop
We now come to the expressions of the form factors considered in ref. [63]. They go beyond the one-loop expressions discussed in the preceding subsection, and read Their structure is easy to understand from the representation (3.8) of the one-loop form factors, as already described in the introduction. The pion form factor, which reduces to unity at lowest order in the low-energy expansion, is extended by the addition of a term linear in z, with a slope given by The same kind of modification is also implemented in the K + π − → π + π − (K S π 0 → π + π − ) vertex: the parameters α + and β + (α S and β S ) correspond to the linear slope and to one of the quadratic slopes (curvatures), respectively, describing the Dalitz plot of the decay K ± → π ± π + π − (K S →→ π + π − π 0 ). Expressed in terms of the parameters introduced in refs. [101][102][103], α + and β + read Using the values determined from data in ref. [104], one obtains (errors have been added quadratically, and the numerical values we use are collected in appendix A for the reader's convenience) α + = −20.84(74) · 10 −8 , β + = −2.88(1.08) · 10 −8 . These values are quite similar to the ones given in ref. [63] (the authors of this last reference use the notation of ref. [105]; ref. [106] gives the correspondence between the two sets of parameters) and used in order to analyze the data in refs. [55][56][57]. Notice that the linear slope α + is determined quite accurately, at less than 4%, whereas the curvature β + is only known with a relative precision of about 35%. Likewise, the parameters α S and β S read The numerical values again follow from ref. [104], see appendix A. They differ somewhat from those that were available to the authors of ref. [63], α S = −5.5(5) · 10 −8 and β S = +0.5(1.3) · 10 −8 . While the form factors W +,S;b1L (z) capture some contributions that arise at order two loops, they do account for all two-loop corrections. This issue will be discussed in detail in section 5 below. Notice also that the relations given in eq. (3.9) generalize in a straightforward manner to the representation (3.15), (3.19)  Figure 3. The left-hand plot shows the data for |W + (z)| 2 , in units of 10 −11 , in the electron mode from the NA48/2 experiment [56] (filled squares, in blue) and from the BNL-E865 experiment [55] (filled circles, in red). The right-hand plot shows the same data, but with the error bars of two points, one for each experiment, around z = 0.3, rescaled by a factor of 2.35. Throughout, only the statistical uncertainties are shown.
In the sequel, we will also consider other representations of W +,S (z) than the one-loop or beyond-one-loop ones. For each of these representations, we will define corresponding parameters a +,S and b +,S upon extending the relations (3.19) to the form factor under consideration.
4 Extraction of a +,S and b +,S from recent data As already mentioned in the introduction, the main recent feature as far as the data are concerned is the availability of quite precise information on the decay distributions for the charged-kaon channels K ± → π ± + − , see table 1. The relation between the decay distribution with respect to the di-lepton invariant mass squared s and the corresponding form factor is given by (recall that λ(x, y, z) denotes the Källen or triangle function, whereas M K stands either for M K ± or for M K S ) (4.1) Here, we will limit our attention to the results coming from the high-statistics experiments [55,56] for the electron mode, and [57] for the muon mode. For the electron channel, the experimental data for |W + (z)| are shown on figure 3. The NA48/2 data [56,57] are available on the HEPData repository [107,108]. 1 Data points for the E865 experiment [55] do not seem to be available on a public repository. For both experiments, the uncertainties given for the individual data points are statistical only. For the systematic uncertainties in the determinations of a + and b + , we refer the reader to the experimental articles quoted in table 1. These data have been confronted with various ansätze for the JHEP02(2019)049 form factor W + (z), from a simple constant plus linear (in z) type of parameterization, W lin (z) = G F M 2 K f 0 (1 + δz), with however little theoretical content, to theoretically more elaborate models [63,[68][69][70]. In the present study, we will only consider the "beyond one loop" form factor of ref. [63].
Fitting the expression (3.15) of the form factor, with the values (3.17) as input, to the data of ref. [55], we obtain (4.2) in reasonable agreement, although with somewhat larger uncertainties, with the values given in table 1 of ref. [55]. Repeating the exercise with the data of ref. [56] leads to the rather surprising outcome Thus, whereas the BNL-E865 data favour the "negative" solution [63], b + < ∼ a + < 0, the NA48/2 data rather point toward the "positive" one, b + > a + > 0, which, with b + about three times as large as a + , is more difficult to understand from a theoretical point of view, as already explained in the introduction. Looking however for a second minimum in the NA48/2 data, we indeed find one, corresponding to the negative solution, with only a slightly higher value of χ 2 ,  Up to now, we have analyzed the data with the form factor in eq. (3.15), but for fixed values, as given in eq. (3.17), of the "external parameters" α + and β + . As already observed, the value of α + is determined quite accurately, the precision on β + being much lower. Redoing the fits for various values of β + with decreasing values of |β + |, see figure 4, we notice that moving β + a few standard deviations away from the central value in eq. (3.17) somewhat improves the quality of the fits. The effect is quite mild in the case of the BNL-E865 data, but somewhat more pronounced in the case of the NA48/2 data. This leads us to consider the possibility of fitting simultaneously a + , b + , and β + . Performing this fit on each experiment separately, we obtain the results in table 2 (only the statistical errors are given).
In agreement with the trend shown by figure 4, the NA48/2 data yield values of β + that tend to be larger than the value quoted in eq. (3.17), obtained from a global fit to the data on the Dalitz-plot structure of the K → πππ decays. One also observes that the χ 2 function JHEP02(2019)049 is quite flat, in the vicinity of its minimum, in the β + direction, which leads to rather large ranges of values. Finally, the values of a + and b + produced by this three-parameter fit to the data come out quite similar to the ones obtained from the two-parameter fits, with however, as could be expected, somewhat larger statistical uncertainties (the effect is, however, more pronounced in the muon mode). If we fit the combined NA48/2 and E865 data (after the rescaling of the statistical uncertainties discussed above) in the electron mode, the 1σ interval of values allowed for β + becomes narrower, but confirms the trend toward larger values of β + than those extracted from the K → πππ data. We conclude this study with the following observations: • The NA48/2 data show a marginal preference for the positive solution, in contrast to the BNL-E865 data, which clearly favour the negative one.
• The data in the electron mode from the two experiments are compatible (up to two data points, one from each experiment, which require a rescaling of their error bars) and can be combined. The combined fit clearly favours the negative solution.
• All fits show a moderate improvement when somewhat larger values of β + (smaller values of |β + |) are considered. Taking the rather conservative attitude of letting β + increase by up to two standard deviations from its value in eq. (3.17), i.e. −4 < ∼ β + · 10 8 < ∼ −0.7, leads to the range of values (as |β + | decreases, |a + | also decreases, while |b + | increases) We finally close this section with a few words about the decay K S → π 0 + − . As shown in table 1, our experimental knowledge of these processes is much more scarce, and comes from the data collected by the NA48/1 collaboration in refs. [59] (electron mode) and [60] (muon mode), and analysed using the expression of the form factor from ref. [63], as given in Eq, (3.15). A combined analysis of the branching ratios of the two lepton modes using JHEP02(2019)049 the values of α S and of β S quoted in ref. [63] gives [60] (10.6). In both cases, due to the large uncertainties, the situation is not very conclusive, even as far as the signs of a S and b S are concerned.

The two-loop representation of the form factors
The discussion at the end of section 3.1 suggests a dispersive procedure for the construction of the form factors, based on the extension at next-to-next-to-leading order of the relation given in eq. (3.12). The successive steps allowing for such a construction, based on chiral counting, analyticity and unitarity, have been described in detail in refs. [109,110] in the case of the pion form factor, and our task will be to adapt them to the present situation.
As far as analyticity is concerned, the form factors W +,S (z) are analytic in the complexz plane with cuts along the positive real axis, starting at z = 4M 2 π /M 2 K . The discontinuity along these cuts is provided by unitarity. For our present purposes, we need only consider the contributions from two pion intermediate states, which is sufficient in order to account for the analyticity properties of the form factors in the range of values of z relevant for the processes K ± (K S ) → π ± (π 0 ) + − . As compared to the case of the treatment of the pion form factor in refs. [109,110], we also need to deal with the fact that the kaon becomes an unstable state in the presence of weak interactions. This has some general consequences as far as the analyticity properties are concerned. In particular, the absorptive and dispersive parts of W +,S (z) for z real are not real, and W +,S (z) do not satisfy the property of real analyticity [111], unlike, for instance, the electromagnetic form factor of the pion or of the kaon. As far as their analyticity properties are concerned, the form factors W +;S (z) are actually quite akin to the form factor f ηπ + (s) describing the isospin-violating τ → ηπν τ second-class transition, which is discussed in detail in ref. [112]. Many aspects related to the analyticity properties of f ηπ + (s) can be directly transcribed to W +;S (z). For instance, the discussion in ref. [112] on the absence of anomalous thresholds in f ηπ + (s) applies, mutatis mutandis, also to W +;S (z).
While the use of dispersive methods might appear as questionable in the presence of an unstable kaon, one should however recall that the computation of the form factors to two loops within chiral perturbation theory rests on the evaluation of Feynman diagrams, which have well-defined analyticity properties. A careful analysis [113,114] shows that these analyticity properties can be reproduced within a dispersive framework upon letting the kaon mass squared become a free variable M 2 K . One starts with a value of M 2 K which lies below the three-pion threshold (the fact that the kaon is also unstable due to its decay into two pions plays no role in the present discussion), so that the kaon becomes stable (with respect to its decay into three pions), and dispersion relations can be implemented. At the end, one moves M 2 K above the three-pion threshold through an analytic continuation, providing the kaon mass squared with a small positive imaginary part, M 2 K → M 2 K + iδ, δ → 0 + . It has in particular been shown [115,116] that this analytic continuation can be performed without encountering singularities in the case where the masses of the charged JHEP02(2019)049 and neutral pions are equal. We will assume this to be the case in the sequel, isospinbreaking effects being far from our present concern. The situation where the difference between the masses of the neutral and charged pions is taken into account would anyhow require a separate analysis.

Construction of the form factors to two loops
Having described the overall framework, let us now turn toward the explicit implementation of the procedure. The starting point of the construction is provided by unitarity. Since the partial-wave projections f Kπ→π + π − 1 (s) at next-to-leading order have a more complicated analytic structure than at the lowest order [117], [112], it is necessary to consider first the situation where M K < 3M π , as discussed previously. This will be understood to be the case from now on. Then the absorptive part of the form factor is real (it coincides with its imaginary part), and the unitarity condition at two loops, restricted to two-pion intermediate states, reads where we have written, for s > 4M 2 π , where σ π (s) = 1 − 4M 2 π s . The structure of the absorptive parts of the form factors W +,S;2L displayed in eq. (5.1) corresponds to the analyticity properties of the two-loop Feynman diagrams shown in figure 5. Notice that whereas the partial-wave projection f Kπ→π + π − 1 (s) is proportional to λ 1/2 Kπ (s) at lowest-order, this is no longer the case at next-to-leading order. Furthermore, the prescription of endowing M 2 K with an infinitesimal positive imaginary part then also specifies a suitable determination of λ 1/2 Kπ (s). The pion form factor at next-to-leading order can be written as [65,66,109,110] (as far as notation is concerned, we follow the last of these references) Re where ϕ +−;+− 1;ππ (s) is the lowest-order P -wave projection of the amplitude for π + π − → π + π − scattering, β being identified with the slope of this amplitude in its expansion around the center of its Dalitz plot, whereas r 2 π V denotes the mean square of the charge radius of the pion, and Then one has It is interesting to notice that the quantity λ

1/2
Kπ (s) does no longer occur in the right-hand side of this last formula. There remains then to compute where ψ +,S (s) describes the one-loop correction to Re f Kπ→π + π − 1 (s), see eq. (5.2). This computation, while straightforward, constitutes the most involved part of the calculation of the absorptive parts of the form factors as given by the right-hand side of eq. (5.1). The details of this calculation will not be given here, but are collected, for the interested reader, in appendix B, together with the notation. Here we simply display the final expressions of JHEP02(2019)049 the form factors at two-loop accuracy in the low-energy expansion, and comment on this result with a few remarks: • The origin of the coefficients ∆α +,S and ∆β +,S is discussed after eq. (B.7) and their expression given explicitly in eqs. (B.37) and (B.38).
• The absorptive parts of the functionsK (λ;0,1) i (s), i = 2, 3, which are given by the functions k i (s) defined in eq. (B.13), develop an imaginary part for s > 4M 2 π when M K takes it physical value. This feature is at the origin of some of the general properties of the form factors discussed at the beginning of this subsection, like the loss of real analyticity. More specifically, the latter follows from the analyticity properties of the second type of Feynman diagrams, with the vertex topology, depicted in figure 5, and which produce the contributions involving the functionsK (λ;0,1) i (s) for i = 2, 3.
• For the same reasons, the constants a +,S and b +,S are in general complex. Their imaginary parts are generated for the first time at the two-loop level, and are thus chirally suppressed. They arise through Feynman graphs with the vertex topology, but also through Feynman graphs without absorptive parts, of the type shown in figure 6. Moreover they will be proportional to the phase space for the K → πππ transition, which is also small. Due to this double suppression, these imaginary parts can be neglected in practice, and will, in particular, not impinge on the analysis in section 4, given the present (and, probably, future) statistical uncertainties of the data.
• Uncommon features, like circular cuts in the complex plane, intersecting the usual unitarity and left-hand cuts on the real axis, also show up in the analyticity properties of the partial-wave projections f Kπ→π + π − 1 (s) computed from the one-loop amplitudes. Their existence follows from the general analysis made in ref. [117], and they are discussed more specifically in refs. [82,115] and [112]. Figure 6. A Feynman diagram that does not contribute to the absorptive parts of the form factors W +,S (z) at two loops, but which gives imaginary parts to a +,S and to b +,S . The meaning of the lines and of the vertices is as in figure 2.

5.2
Comparing W +,S;2L (z) and W +,S;b1L (z) The one-loop expression for W + (z) does not give a very good description of the data. If the latter are fitted, as in section 4, with W +;1L (z), keeping a + as a free variable, the quality of the fit deteriorates tremendously. Typically, the value at the minimum of the χ 2 function is in the range between 300 and 500. The representation W +;b1L (z) thus constitutes a real improvement in the description of the data. Likewise, one may legitimately ask oneself to which extent the full two-loop expression W +;2L (z) constructed above would further modify this picture.
In comparing the full two-loop representations (5.8) of the form factors with the expressions (3.15) of ref. [63], we notice that the latter are essentially reproduced by the first line of eq. (5.8), which follows if in eq. (5.1) we restrict the pion form factor and the P -wave projections to their polynomial parts: taking, for a π V , the estimate a π V = 1/M 2 V given by vector-meson dominance. What is missing in order to completely reproduce the expressions of the form factors W +,S;b1L (z) is the term proportional to the product β +,S × a π V , which does not occur in eq. (5.8). This term actually represents a contribution of order three loops, whence its absence from eq. (5.8). Furthermore, if we go once more through the fits in section 4 with this term omitted from eq. (3.15), the outcomes for a +,S and b +,S are barely changed, so that, from a phenomenological point of view, this term is not important in the region of z corresponding to the phase space of the K → π + − transitions.
We may now proceed with the quantitative comparison between the full two-loop expressions given by eq. (5.8) and the form factors of ref. [63]. This comparison is shown in figure 7 for the modulus squared of W +;2L (z) and W +;b1L (z), for two sets of values for the parameters a + , b + and β + . This difference remains quite small, and in any case well below the present statistical uncertainties of the data. This feature is shared for other values of the parameters taken in the ranges discussed in section 4. Differences between W +;2L (z) and b +,S from data will thus not be modified in any substantial way if, instead of W +,S;b1L (z), one uses, in section 4, the full two-loop amplitudes W +,S;2L (z).

W + (z) beyond the low-energy expansion
While the two-loop representations of the form factors, or the truncated versions thereof proposed in ref. [63], provide an appropriate description of the experimental data in terms of two sets of parameters a +,S and b +,S , the latter cannot be predicted within the lowenergy framework considered in sections 2 and 5. In order to obtain predictions for them, it is necessary to set up a phenomenological description of the form factors that goes beyond the low-energy framework. If such a description of the form factors becomes available, the values of the constants a +,S and b +,S can then be obtained through the definitions given in eq. (3.19). The purpose of the present section is to proceed with such a phenomenological construction of the form factor W + (z). Building on the results obtained so far, we will propose a simple model for the form factor W + (z) that accounts for rescattering effects in the two-pion intermediate state beyond the framework set by the low-energy expansion, and, at the same time, provides a matching to the short-distance behaviour investigated in section 2.2. Consequently, our model will consist of three parts, W + (z) = W ππ + (z) + W res + (z; ν) + W SD + (z; ν). (6.1)

JHEP02(2019)049
Before going into the details of their respective evaluations, let us briefly describe the physical content of each of these parts. As already mentioned, the first part describes the contribution from the two-pion intermediate state to W + (z). It is constructed upon assuming, in analogy with the electromagnetic form factor of the pion F π V (s) [98,118], that it is given by an unsubtracted dispersion integral, The absorptive part consists of the two-pion spectral density ρ ππ + (s), and is obtained upon inserting a two-pion intermediate state in the representation of the form factor given in eq. (2.2), In order to evaluate this absorptive part, we require two ingredients: a representation of the pion form factor F π V (s), and a representation of the P -wave projection f Kπ→ππ 1 (s), which both extend to the whole energy range set by the cut singularity of W + (z). We will provide an explicit realization of these two quantities in subsection 6.1 below. For the time being, let us just notice that for the convergence of the unsubtracted dispersive integral it is sufficient that their product is bounded by, say, a constant for large values of s.
Above this lowest threshold, several intermediate states will contribute to the discontinuity of W + (z) in the 1-GeV region and beyond. Considering only two-meson intermediate states, the next thresholds will come from K + π − or K 0 π 0 intermediate states, followed by K 0K0 , K + K − , and so on. As the energy increases, the number of possible exclusive intermediate states grows, and they eventually merge into the inclusive contribution provided by the QCD continuum, as discussed in section 2.2. We will describe this process in terms of an infinite tower of equally-spaced zero-width resonances where M ∼ 1 GeV is the scale of the lowest resonance occurring in this tower. The main task facing us will be to find an appropriate Regge-type resonance model, which reproduces the correct QCD short-distance properties of the form factor. This issue will be adressed in section 6.2. In particular, the dependence of W res + (z; ν) on the shortdistance scale ν has to match the same dependence that appears in the third part on the right-hand side of eq. (6.1), simply given by the factorized contribution coming from the Q 7V operator, This last contribution is evaluated in section 6.3.

The contribution from the two-pion state
Assuming we have a representation of the form (6.2) at our disposal, the contributions a ππ + and b ππ + from the two-pion intermediate state to the coefficients a + and b + , respectively, are obtained through two sum rules that follow from the definitions given in eq. (3.19), . (6.7) As far as the convergence of the integral in (6.2) is concerned, it depends on the behaviour of both the pion form factor F π V (s) and the partial wave projection f Kπ→π + π − 1 (s) for large values of s. As already mentioned, in order to evaluate the sum rules in eq. (6.7), we need representations of the pion form factor F π V (s) and of the partial-wave projection f π + π − →K + π − 1 (s) that extend beyond their low-energy expansions. There exist several "unitarization" procedures that precisely allow to do this. One of them, the inverse-amplitude method (IAM) [119,120], gives quite reasonable results when applied to ππ scattering in the P -wave and to the pion form factor [121,122] in the region of the ρ meson. The result, obtained starting from the one-loop low-energy expression of F π V (s), provides a rather simple representation of F π V (s). Its phase is the phase of the P -wave projection of the ππ scattering amplitude, as required by Watson's final-state theorem, and it reproduces the one-loop expression of F π V (s) for small values of s. Furthermore, it exhibits a resonant behaviour at the expected value of the energy squared, i.e. s ∼ M 2 ρ . It is quite easy to understand in simple terms how this last property emerges from the representation (6.8). Neglecting at first the contribution from the pion loops, materialized by the term proportional to β in the denominator, this expression reduces to the well-known VMD form It has the expected pole at s = M 2 V ∼ M 2 ρ , and the real part of the pion-loop contribution will slightly change the real part of the pole position, while the imaginary part of the pion loop will move it from the real axis to the second Riemann sheet, leading to the resonant behaviour shown in figure 8.
We next turn toward the partial-wave projection f K + π − →π + π − 1 (s), to which we wish to apply the same procedure. Starting from its one-loop expression obtained in appendix B, and neglecting the contributions from the circular cuts, the IAM-unitarized version of eq. (B.15). The remaining terms give only a small correction in comparison. We thus end up with . (6.10) If we transpose the discussion that follows eq. (6.8) to the representation (6.10), we observe that, when the contributions from the pion loops are discarded, the "bare" pole is located at a value of s given by s pole = s 0 + (α + /β + )M 2 π . For the values given in eq. (3.17), i.e. α + /β + ∼ 7, this gives s pole ∼ 12M 2 π . In order to obtain a pole located at s pole ∼ M 2 ρ , one needs a higher value of the ratio α + /β + , say α + /β + ∼ 25. In view of the error bars of the numerical values of α + and β + , this can be achieved most economically upon increasing β + by about two standard deviations from its central value in eq. (3.17). Assigning even part of the effect to an increase of the absolute value of α + would represent a much more significant deviation from its value in (3.17). At this stage one might recall the discussion in section 4, where it was already noticed that upon keeping β + as a free variable to be fitted to the data, the outcome was favouring values deviating from the one in (3.17) by a similar amount. It would certainly appear as somewhat far-fetched to ground the reason for preferring larger values of β + on the simple characteristics of the IAM-unitarized partial wave f K + π − →π + π − 1 (s). Nevertheless, the concordance, on this issue, with the analysis presented in section 4 is definitely interesting and noteworthy. 2 In figure 9, we illustrate 2 For the sake of completeness, let us mention that the Dalitz plot of the K + → π + π + π − transition has been measured with high precision by the NA-48/2 Collaboration [123]. These results were published after ref. [104]. Using only ref. [123] does not allow for a separate determination of α+ and β+ but only of the ratio α+/β+. In terms of the Dalitz plot parameters g, h and k of ref. [123], one obtains α+/β+ ∼ −g/(h − 3k) = 6.53 (1) in agreement with the value obtained from eq. (3.17). the evolution of |f K + π − →π + π − 1 (s)| as given in eq. (6.10) for different values of β + and for s ≥ 4M 2 π . Notice also that the approximation considered in eq. (6.10) preserves Watson's final-state theorem, and the phases of F π V (s) and of f K + π − →π + π − 1 (s) should be identical. For the choice β = 1.11 (this value of β is itself chosen such as to match the phase of the pion form factor to the phase of the P -wave projection of the ππ amplitude, both obtained upon unitarization of their one-loop expressions by the IAM method), this is the case for β + = −0.85 · 10 −8 .
The numerical evaluation of the two sum rules in eq. (6.7) for these values of β and of β + then gives a ππ + = −1.58, b ππ + = −0.76. (6.11) In this approach, the overall negative sign of these numbers is driven by the negative sign of α + . The result for b + comes out relatively close to the values extracted from the data in section 4, whereas the absolute value of a + is about twice as large. Of course, in both cases there are other contributions which we need to discuss before being in a position of making a more definite statement. Let us finally notice that the same approach can also be implemented in order to evaluate the contribution from the two-pion state to a S and b S . It suffices to make the appropriate replacements in eq. (6.10). We first notice that a ratio α S /β S ∼ 25 lies (almost) within reach for the values indicated in eq. (3.18), e.g. α S ∼ −7.5 · 10 −8 , β S ∼ −0.3 · 10 −8 .
The resulting values of a ππ S and b ππ S would then also be negative, and about three times smaller in absolute value than those obtained for a ππ + and b ππ + in eq. (6.11).

Intermediate states with higher thresholds
In section 2.2, we have established the high-energy behaviour of the dimensionally renormalized form factor, see eq. (2.28). Reproducing this behaviour is necessary in order to obtain an amplitude that does no longer depend on the short-distance renormalization scale ν once the contribution from the local factorized operator Q 7V is added. This behaviour

JHEP02(2019)049
is, however, not reproduced by the contribution of the ππ intermediate state that we have just studied. It has therefore to come from the remaining infinite number of intermediate states, with higher and higher thresholds, and which eventually end up into the perturbative contribution of quarks and gluons at short distances. One might attempt to describe some of these additional intermediate states in a similar treatment as the one adopted here for the two-pion states. This holds, in particular, for the next thresholds, due to Kπ and toKK intermediate states. We leave such improvements for future work. Here, we will adopt a simpler point of view, where the additional intermediate states are described by zero-width resonance states. Such a picture would naturally emerge, for instance, from the perspective of the limit where the number of colours N c becomes large [124,125]. Contributions of this type are depicted on figure 10, and would produce a form factor with the following expression (6.12) The first sum runs over resonances with quantum numbers J P C = 1 −− , I = 1, S = 0. It starts here with the φ(1020), since the ρ meson is already contained in the contribution of the ππ intermediate states that we have discussed previously. The second sum runs over the resonances with quantum numbers J P = 1 − , I = 1/2, S = ±1, and starts with the K * (892). The strong couplings f V and g V are fixed, for instance, by the widths of the decays like φ → e + e − and K * → Kπ, respectively. Information on the weak couplings f V and g V is, however, not available. This represents the usual difficulty in obtaining reliable estimates of the low-energy constants in the weak sector through resonance saturation [93][94][95][96][97]. Another difficulty lies in the fact that the correct short-distance behaviour (2.28) can only be recovered upon considering an infinite number of resonances [125]. This second difficulty can be dealt with upon using available techniques [126][127][128] to construct resonance models with spectra of the Regge-type, and with residues that can be tuned such as to build up harmonic sums that can often be resummed exactly and, moreover, reproduce the prescribed asymptotic behaviour. We will present such a construction in the case of interest here. In order to set the stage, let us go back to the calculation in section 2.2. It teaches us that a dispersive representation of the type would fulfill the required conditions. Here A is a normalization constant, and M is a mass scale which corresponds to the onset of the perturbative continuum due to the quark loop, and which we would like to identify with a typical resonance scale M ∼ 1 GeV. Of course, the above dispersive integral does not converge, and we first need to consider the dimensionally regularized version of the spectral density ρ M (s), namely [129] ρ M (s; After renormalization in the MS scheme, one therefore finds The limit of large space-like values of z gives The short-distance behaviour given in eq. (2.28) is then recovered in this case with the choices

JHEP02(2019)049
Actually, the dispersive integral (6.13) with ρ M (x; D) can be done explicitly (w ≡ −zM 2 K /M 2 ): For the definition and properties of the hypergeometric function 2 F 1 , we refer the reader to ref. [130]. The properties (6.15) and (6.17) can then be directly recovered from those of this function, see e.g. [130, eq. 15.8.2], (6.20) It is possible to reproduce the salient properties of the simple model discussed above through a Regge-type resonance model of the form provided one can find a set of functions µ n (D) that satisfies the following requirements: Here ζ(x) is the Riemann ζ-function and Li y (x) denotes the polylogarithm function, defined as Li y (x) = n>0 x n n −y for |x| < 1 and arbitrary complex y, and by analytic continuation for other values of x. The integrand in the relation (6.25) has a pole at u = 0, coming from the pre-factor only, since the terms inside the square brackets are well behaved at u = 0, with Li −1 (x) = x/(1 − x) 2 . According to the Converse Mapping Theorem [131,132], from this pole at u = 0, located on the left of the fundamental strip 0 < c 2 < 4−D 2 , one deduces that ∞ n=1 a n D−4 2 remains finite as this limit is taken, the integrand has also a pole at u = 4−D 2 , due to the first term in the square brackets. In this case, the pole being located on the right of the fundamental strip 0 < c 2 < 4−D 2 , the Converse Mapping Theorem allows us to state that ∞ n=1 a n D−4 2 We thus conclude that we are able to build a Regge-type resonance model, defined by eq. (6.22), with 28) and which satisfies all the required properties. The part of the integral that remains finite in the limit D → 4 can be resummed explicitly, and we find with the coefficients c p (D) given by c 0 (D) = 1 and, for p > 0, by In the limit where D → 4, one obtains which indeed corresponds to eq. (6.29). Summarizing, we end up, after minimal subtraction, with the expression Accordingly, the corresponding contributions to a + and b + read (cf. eq. (3.19)), i.e.
a res The expression of b res + (ν) involves the slope at the origin λ + of the form factor f K ± π ∓ + (s). It is defined in eq. (A.1), and numerical values for both f K ± π ∓ + (0) and λ + are given in eq. (A.3).
We draw attention here to the fact that only the order O(α 0 s ) contribution to the short distance behaviour of W (s; ν) is reproduced by the resonance model. This means that eq. (6.34) satisfies both (2.14) and (2.28) at order O(α 0 s ) only, which implies that some scale dependence will remain. We could improve on this aspect upon including O(α s ) corrections, which are almost completely determined from the renormalization-group analysis in section 2.2. We would then need to build a resonance model that reproduces both the correct ln(−s/ν 2 ) and ln 2 (−s/ν 2 ) behaviours at high energy. We leave such an improvement for a future study.

JHEP02(2019)049
6.3 The contribution from the factorized Q 7V operator It remains to evaluate the last contribution, due to W SD + (z; ν). From its definition in eq. (6.6), combined with eq. (3.19), we obtain right away Values for C 7V can be found in ref. [89]. In particular, we use C (NDR,HV) 7V (1 GeV)/α = (−0.037, 0.000). Here we have set τ = 0, i.e. we have identified C 7V with z 7V , and we have chosen Λ (4) MS = 300 MeV. In order to investigate the dependence of the sums a res + (ν)+a SD + (ν) and b res + + b SD + (ν) on the short-distance scale ν, we use the lowest-order evolution equations, neglecting the mixing with the penguin operators, i.e.
Here α s stands for the running QCD coupling for N f = 3 active flavours, and for the numerical evaluation, we use Λ Collecting the various contributions to a + and to b + from the model proposed in this section, we end up, according to eq. (3.19), with the following expressions:  Figure 11. The evolution of a + and b + with respect to ν in both NDR and HV schemes, for M = 1 GeV.
Numerically, for 1 GeV ν 2 GeV, M = 1 GeV, and considering only contributions from C 1 and C 2 , one obtains The residual ν-dependence of a + and b + is depicted in figure 11. This result should mainly be viewed as a first serious attemp to evaluate a + and b + . Improvements, for instance on the description of W ππ + (z), or the inclusion of QCD corrections in W res + (z), are clearly required before realistic error bars can be assigned to the values displayed in eq. (6.42). We come back to these issues in the next section. Nevertheless, we find it encouraging that the outcome of the rather simple approach followed in the present section comes rather close to the values obtained from the experimental data in section 4. This is particularly true for b + .

Summary, conclusions, outlook
In this final section, we wish to summarize the content of this article, going successively through the main aspects of the issues that have been addressed, roughly following the list of items given at the end of the introduction. For each item, we provide conclusions and/or critical remarks, as well as an outline of perspectives for future improvements.
7.1 Extracting a +,S and b +,S from recent data Our first concern was to establish the values of the constants a +,S and b +,S that are provided by recent experimental data. In the case of K ± → π ± e + e − , we have shown that a combined fit to the data on the decay distribution from the two high-statistics experiments BNL-E865 and NA48/2 clearly favours the solution where a + and b + are both negative, with |a + | and |b + | comparable in size.
We have also performed fits to the data keeping, in addition to a + and b + , the curvature β + of the K ± → π ± π + π − Dalitz plot as a free parameter, and have found that somewhat JHEP02(2019)049 smaller (in absolute value) values than those obtained by direct determinations from K → πππ data are preferred.
In order to make comparisons between different approaches more convenient, we have introduced intrinsic definitions of the coefficients a +,S and b +,S in terms of the values of the form factors and of their slopes at the origin z = 0.

The low-energy expansion of the form factors to two loops
The extraction of a +,S and b +,S from data was done using the expressions W +,S;b1L (z) of the form factors given in ref. [63]. Our second concern was then to establish whether two-loop corrections not accounted for by these expressions might have an effect on these determinations.
We have provided the full two-loop expression of the form factors, taking only the singularities due to two-pion intermediate states into account, in analogy with ref. [63]. These two-loop expressions are based on the complete one-loop expression of the partialwave projections f Kπ→π + π − 1 (s), which include also ππ rescattering in the crossed channels. These features are not accounted for by the expressions of W +,S;b1L (z) given in ref. [63].
From a numerical point of view, these effects turn out to be quite small in the region of z corresponding to the phase space for the K → π + − decays. In practice, one may thus use the simpler form W b1L (z) in order to analyze the data, instead of the full two-loop expression W 2L (z), without significant impact on the determinations of a +,S and b +,S . This also strongly suggests that still higher-order corrections are quite small and are likely not to modify the picture in any substantial way.
We have also pointed out that the existing one-loop calculations suggest that kaon loops could possibly have a sizeable effect on a + and a S . However, making a quantitative statement on this issue would require a complete two-loop calculation in the usual framework of three-flavour chiral perturbation theory, including the computation of those Feynman graphs that do not exhibit non-trivial analyticity properties.

Contribution from the two-pion state
We have addressed the phenomenological evaluation of the contribution from the twopion state to a + and b + upon writing an unsubtracted dispersion relation for the form factor W + (z), from which sum rules for a + and b + can be obtained. The absorptive part of the dispersion relation is provided by the electromagnetic form factor of the pion F π V (s) and by the P -wave projection f Kπ→π + π − 1 (s) of the K ± π ∓ → π + π − amplitude. For these, we have used a very simple approach, where these two quantities are constructed through unitarization with the inverse amplitude method of their one-loop expressions in the chiral expansion. Nevertheless, this simple description leads to numerical results that lie in the ballpark of the values extracted from data. Constructing or using more realistic representations of F π V (s) and of f Kπ→π + π − 1 (s) is certainly an aspect where improvements are possible.
Indeed, there exist in the literature more elaborate representations of the electromagnetic form factor of the pion that describe data in a wide range of momentum transfer, see for instance refs. [109,118,[133][134][135][136][137][138][139] for a representative sample. In the case of the JHEP02(2019)049 partial waves f Kπ→π + π − 1 (s), recourse to a similar data driven description is unfortunately not possible. Moreover, the simple IAM unitarization procedure we have considered does not appropriately account for their full analyticity structure. More involved methods, like numerical implementation of the Khuri-Treiman representation [140], which has been used in other instances, see for instance refs. [112,141,142], are available and would represent a significant improvement in the description of these partial waves beyond the low-energy region. It would also be interesting to see whether the problem caused, within the IAM, by the too small value of α + /β + persists when a different unitarization method is considered. Another interesting possibility would be to also include the KK intermediate states into the dispersive representation, which however requires a two-channel analysis [143].

Matching with the short-distance regime
Whereas the form factors W + (z) and W S (z) are clearly dominated by low-or intermediateenergy physics, the time-ordered product of the electromagnetic current with the lagrangian density for |∆S| = 1 transitions is singular at short distances and needs to be renormalized. This renormalization is implemented through the operator Q 7V and its Wilson coefficient C 7V (ν). As a result, the form factors behave, in the asymptotic Euclidean region, as ∼ ln(−s/ν 2 ), where ν is the renormalization scale. At the phenomenological level, this short-distance behaviour results from the pile-up of more and more complicated intermediate states, with higher and higher thresholds, until the region where the QCD continuum sets in is reached. We have described this process through a, necessary infinite, set of zerowidth resonances. In the absence of QCD corrections, we have shown that it is possible to adjust the couplings of these resonances such as to reproduce the correct high-energy behaviour. Working only at lowest order leaves a rather strong sensitivity to the subtraction scale. Extending the resonance model in order to account also for the order O(α s ) QCD effects in the high-energy part, which are to a large extent known from the renormalization group argument given in section 2.2, would probably reduce this dependence on the short-distance scale.

Conclusion
This study was undertaken with the aim of exploring the possibility to achieve a determination of the constants a +,S and b +,S , describing the decay distribution of the K ± (K S ) → π ± (π 0 ) + − decay modes, such as to assess, through confrontation with present and forthcoming experimental data, the amount (if any!) of violation of lepton flavour universality in the kaon sector. We hope that our study demonstrates that such a determination based on a phenomenological approach is possible, with a reasonable amount of theoretical work, and that there is room for improvement in several of the aspects that contribute to it. Such an endeavour would then be complementary to existing and future efforts to address this issue through numerical simulations of QCD on the lattice [48][49][50]. We plan to come back to some of the aspects involved in such a phenomenological determination and discussed above in future work.

JHEP02(2019)049
Acknowledgments This work started while the three authors were attending the "NA62 Kaon Physics Handbook" meeting at the Mainz Institute for Theoretical Physics. The authors would therefore like to express their gratitude to MITP, its Director and its staff for their warm hospitality and for partial financial support. M.K. and D.G. thank the INFN-Sezione di Napoli and the Università di Napoli Federico II for their hospitality during several stays that made the completion of this work possible, as well as E. de Rafael for interesting and informative discussions. One of us (D.G.) also thanks the CPT for its hospitality. We also thank A. Nath for his active participation in some earlier stage of this project, E. Goudzovski for informative discussions, as well as M. Hoferichter, S. Kettell and L. Tunstall for useful correspondence.

A Numerical values
In this appendix, we provide the numerical input values for several quantities that have been used in the text. For the reader's convenience, some of them have been gathered in table 3. In addition, we also need the values of the form factors f K ± π ∓ + (s) and f K S π 0 + (s), and of their derivatives, at the origin s = 0: (A.1) We have taken their values from the analysis of ref. [144], which gives These values belong to the ranges determined from data in ref. [145].

B The computation of ψ +,S (s)
In this appendix, we compute the function ψ +,S (s), which describes the one-loop correction to Re f Kπ→π + π − 1 (s), the real parts of the P -wave projections of the amplitudes for the processes Kπ → π + π − , where Kπ stands either for K + π − or for K S π 0 , see eq. (5.2). In order to obtain ψ +,S (s), one thus first needs to construct these amplitudes at one loop. This can be done, for M K < 3M π , within an iterative construction, following the method  Table 3. Numerical values used for the various input parameters. The first seven entries are taken from ref. [58], and the values of the constants g 8 and g 27 come from ref. [1]. The fit of ref. [104] provides the values for the K → πππ Dalitz-plot parameters α 1 , . . . ξ 3 , which are given in units of 10 −8 .
that has been described several times, e.g. in refs. [146][147][148][149] and [110,115]. Indeed, up to and including two loops, the two amplitudes in question have the general structure The absorptive parts of the functions W 0,1 (s) are given in terms of the absorptive parts along the right-hand cut of the lowest S and P partial-wave projections of the corresponding Kπ → π + π − amplitudes: Similar expressions hold for the remaining functions W 0,1;t (s) and W 0,1;u (s), involving the absorptive parts of the amplitudes in the crossed s ↔ t and s ↔ u channels, respectively. If only two-pion intermediate states are considered, then the absorptive parts of the various functions appearing in the formula (B.1) will be given, at one loop, by the products of tree-level amplitudes for ππ scattering and for the processes Kπ → π + π − , which are simply first-order polynomials in the Mandelstam variables. In particular, the lowest-order Kπ → π + π − amplitudes are expressed in terms of the Dalitz-plot parameters α 1,3 , β 1,3 and JHEP02(2019)049 γ 3 , in the nomenclature of refs. [103] and [104]. The lowest-order ππ scattering amplitudes are likewise expressed in terms of the two subthreshold parameters α and β [146,147]. The projection on the P partial wave is given by  (s), given in eq. (5.4), is a polynomial of first order in s. At next-to-leading order, the contribution involving the polynomial P(s, t, u) can be given the following parameterization:  We have checked that they agree with the ones given in [115].
One may notice that the functions k 2 (s) and k 3 (s) involve λ

JHEP02(2019)049
We can now evaluate the coefficients ∆α +,S and ∆β +,S defined in eq. (B.10). With the numerical input provided in table 3 (we do not distinguish between the charged and the neutral pion masses), we then obtain