Improved Standard-Model prediction for K L → ℓ + ℓ −

: We present a comprehensive calculation of the K L → γ ∗ γ ∗ form factor in dispersion theory, using input from the leptonic decays K L → ℓ + ℓ − γ , K L → ℓ +1 ℓ − 1 ℓ +2 ℓ − 2 , the hadronic mode K L → π + π − γ , the normalization K L → γγ , and the matching to asymptotic constraints. As key result we obtain an improved determination of the long-distance contribution to K L → ℓ + ℓ − , leading to the Standard-Model predictions Br [ K L → µ + µ − ] = 7 . 44 +0 . 41 − 0 . 34 × 10 − 9 , Br [ K L → e + e − ] = 8 . 46(37) × 10 − 12 , and more stringent limits on physics beyond the Standard Model. We provide a detailed breakdown of the current uncertainty, and delineate how future experiments and the interplay with lattice QCD could help further improve the precision.

For K L → ℓ + ℓ − the roles become reversed, and the CP -conserving LD contribution now proceeds via the C amplitude in Eq. (1.1).Normalizing to the K L → γγ decay, the  decay rate can be written as where the dominant imaginary part comes from γγ intermediate states [22] Im . (1.6) The real part is typically decomposed into the loop function that corresponds to the oneloop result in ChPT and a local contribution χ(µ) [12,23,24] Re where χ(µ) receives both LD and SD contributions in the SM, 1 see Figs. 1 and 2. Fixing the sign following the arguments from Refs.[12,23,30], the SD part is known to about 4% precision [31,32] χ SM SD = −1.80(6), (1.8) see App.A, with uncertainties dominated by CKM matrix elements.The experimental results [6,[33][34][35][36] Br[K L → µ + µ − ] = 6.84 (11)  would thus allow one to derive BSM constraints if the LD contribution to χ(µ) could be calculated.
of the K L → γ * γ * form factor informed by the radiative decays K L → ℓ + ℓ − γ [38,39] and the asymptotic behavior expected from a partonic calculation [12].More recently, also first steps towards a lattice-QCD calculation of χ(µ) have been taken [40][41][42][43], and proposals were developed to extract SD information from the time evolution of the K → ℓ + ℓ − decay rate [44][45][46].In view of the experimental prospects for future precision measurements of K L decays in phase 2 of HIKE [47,48] and at KOTO II [49,50] (see also Ref. [51,52]), it is thus timely to revisit the calculation of the K L → γ * γ * form factor using methods that were recently developed in the context of pseudoscalar-pole contributions to hadronic light-by-light scattering in the anomalous magnetic moment of the muon [53][54][55][56][57][58].
To this end, we analyze the K L → γ * γ * form factor in dispersion theory, which allows us to not only include the leptonic modes 2 in its reconstruction from data, but, crucially, also the hadronic channel K L → π + π − γ.In addition, we perform the matching to asymptotic constraints via a dispersive representation of the respective loop integrals, which allows for a better separation of their low-and high-energy components.This strategy follows closely our calculation of the π 0 → γ * γ * transition form factor [57][58][59][60][61], which was mainly motivated by the pion-pole contribution to hadronic light-by-light scattering, but later applied also to specific channels in hadronic vacuum polarization [62][63][64][65][66] and, most importantly in the context of K L → ℓ + ℓ − , to an improved SM prediction for π 0 → e + e − [29].In addition, the calculation presented here profits from the experience gained with the analog form factors in η, η ′ → γ * γ * [67][68][69][70] and The outline of our analysis is as follows.We first derive a dispersive representation of K L → γ * γ * in Sec. 2, which is then matched to asymptotic constraints in Sec. 3 to arrive at the final representation given in Sec. 4. The resulting SM prediction for K L → ℓ + ℓ − is derived in Sec. 5, leading to the BSM constraints discussed in Sec.6, before concluding in Sec. 7. Details of the SD contribution in the SM, the dispersive representation of the loop integrals, and the integration kernels for K L → ℓ + ℓ − are discussed in the appendices.

Leptonic processes
Following the conventions of Ref. [1], the amplitude for K L → γ * γ * can be written as since the other two Lorentz structures, associated with scalar amplitudes a(q 2 1 , q 2 2 ), b(q 2 1 , q 2 2 ), only lead to CP -violating contributions.Accordingly, the scalar function c(q 2 1 , q 2 2 ) defines the K L → γ * γ * transition form factor of primary interest to us.Its normalization is determined via the on-shell process The spectrum for the singly-virtual decays with z = q 2 /M 2 K , q 2 the virtuality of the dilepton pair, f (z) = c(q 2 , 0)/c(0, 0), and spectral function Similarly, the doubly-virtual process 2 )/c(0, 0), according to where λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc).However, because of the kinematic constraints and statistics of these decays, both processes can, in practice, only provide useful information about the slope parameter b K L defined as ). (2.6) Among the two popular parameterizations of the form factor, Bergström, Massó, and Singer (BMS) have proposed [74,75] where ) with C K * = 2.5 and β K * = 1 for the resulting BMS model [75].The form in Eq. (2.8) is constructed in such a way that A K * (0) = 0, which reflects the fact that the decay K → πγ is not allowed by angular momentum conservation, and thus an on-shell K * ≃ Kπ → γ transition cannot occur either.Phenomenologically, this is an important feature of the K * contribution, as it can be sizable for the slope of the form factor, without altering the normalization.The singly-virtual decays constrain α K * = −0.217(34) (dominated by Ref. [76]) and α K * = −0.158(27) (dominated by Ref. [77]), for the electron and muon channel, respectively.There is also a measurement from the doubly-virtual process where ℓ 1 = ℓ 2 = e [78], which finds α K * = −0.14(22) with a large uncertainty after assuming factorization for the form factor.The world average, α K * = −0.205(22) [6], then leads to (2.9) Second, D'Ambrosio, Isidori, and Portolés (DIP) have proposed the parameterization [38] often subject to the high-energy constraint 1 + 2α DIP + β DIP = 0.In addition to the singly-virtual modes [76,77], the parameter α DIP was also extracted from the ℓ 1 = e, ℓ 2 = µ channel [79] assuming β DIP = 0.The global average, α DIP = −1.69(8)[6] is dominated by the K L → e + e − γ mode, leading to b K L = 2.81(13) GeV −2 , consistent with the extraction (2.9) from the BMS model within uncertainties.
Our aim is to improve model parameterizations of c(q 2 1 , q 2 2 ), by developing a dispersive representation for this form factor in analogy to Refs.[57,58].In particular, this allows us to profit not only from the leptonic decay data, but also from data on the related hadronic process K L → π + π − γ, to which we turn next.

K
γ denoting the photon energy in the kaon rest frame.As a first step, we rewrite these conventions in terms of Mandelstam variables [68,69] and the center-of-mass angle leading to the identifications Moreover, the difference of the z i is related to z t according to . (2.15) The two amplitudes are typically decomposed into "inner-bremsstrahlung" (IB) and "direct emission" (DE), where the former only contributes to the electric amplitude E(z i ) which is thus CP violating, but still significant due to the infrared enhancement.The remaining contribution can be expanded in terms of multipoles where the odd electric and even magnetic multipoles are CP violating.This can be contrasted with the expected partial-wave expansion [82], e.g., with derivatives of the Legendre polynomials P ′ ℓ (z), so that requiring CP invariance indeed only leaves the odd partial waves in M (z i ) (and vice versa for E(z i )).To constrain the CP -conserving K L → γ * γ * process via the ππ singularities, we are thus mainly interested in the first multipole of M (z i ), which takes the role of the partial wave f 1 (t) in Refs.[68,69].
Experimentally, M 1 (z 3 ) was studied in Refs.[83][84][85] for q 2 = 0, using a ρ-pole ansatz of the form [86] , and the rest of the normalization was expressed in terms of the K S → π + π − decay [87,88] 3 Assuming no interference with the IB contribution, a branching fraction of Br[K L → π + π − γ (DE)] = 2.84(11) × 10 −5 is given in Ref. [6] (based on the K L → π + π − measurement from Ref. [89] and the ratio of IB vs. DE branching fractions from Refs.[83][84][85]).Finally, there is a measurement of the branching fraction Br[K L → π + π − e + e − ] = 3.08(20) × 10 −7 [90].Apart from branching fractions, experiments also performed analyses of the double differential decay rate ) , which facilitate spectral-shape studies of the magnetic amplitude M .In this regard, a representation such as Eq.(2.19) cannot be justified from a dispersive perspective, since the sum of a contact term and a ρ pole violates Watson's final-state theorem [91].To improve upon existing parameterizations the next step thus consists of restoring unitarity as a minimal requirement.

Unitarity and dispersion relations
As a first step, we decompose c(q 2 1 , q 2 2 ) into functions with definite isospin and a term that collects the K * contribution as detailed in Fig. 3: . Different contributions to the K L → γ * γ * form factor in a dispersive approach.The gray blobs denote the respective hadronic matrix elements, and the black square the strangenesschanging K * -γ * coupling that is again mediated by vector-meson states. where and the two photons have isovector (v) and isoscalar (s) quantum numbers as indicated by the subscripts.Bose symmetry implies c vv (q 2 1 , q 2 2 ) and c ss (q 2 1 , q 2 2 ) to be symmetric under the exchange of the arguments, and for the mixed terms c vs (q 2 1 , q 2 2 ) = c sv (q 2 2 , q 2 1 ), etc. Phenomenologically, the necessity of the K * contribution follows from the large slope parameter (2.9 58,61] and η, η ′ → γ * γ [67-70] not far from the typical scale M −2 ρ .In the same vain, M (t, q 2 ) = M 1 (z 3 )/M 3 K , the (normalized) multipole for general photon virtuality q 2 , is decomposed as ( The two-pion cuts then lead to the unsubtracted dispersion relations where q π (x) = x/4 − M 2 π , F V π is the electromagnetic form factor of the pion, and 2 GeV denotes a cutoff that separates the low-energy contribution.

Factorization
The simplest solution to the decomposition (2.24) involves separating the final-state interaction into an Omnès factor [92] determined by the P -wave ππ scattering phase shift δ 1 (t), and parameterizing the remainder via polynomials with P v/s (t) for the isovector/isoscalar contribution, a parameter ϵ K * for the K * , and another function a v/s (q 2 ) to describe the dependence on the photon virtuality.In addition, where the imaginary part follows a Breit-Wigner description of the K * resonance. 4Using the ansatz (2.27), the c vv contribution in Eq. (2.25) can be rewritten as 2 )a vv (q 2 1 ), The first line implies the factorization a vv (q 2 ) = α 1 a v (q 2 ).Moreover, Eq. (2.25) remains invariant upon rescaling It is further instructive to consider the limit of a narrow resonance [60] a vv (q 2 ) → eF where the KFSR relation [96,97] has been used to trade the ρππ coupling in favor of the pion decay constant F π .Up to a normalization, a vv (q 2 ) thus reduces to a ρ propagator in the narrow-resonance approximation.
The isoscalar function a s (q 2 ) is dominated by ω and ϕ resonances, leading to the ansatz with a spectral function determined by a suitable parameterization of the imaginary parts following [57,58,61]. 5In the limit of a narrow resonance, this form becomes identical to Eq. (2.31) upon identifying c ρ = eF 2 π P v (M 2 ρ ).For the mixed isovector-isoscalar contribution we have 2 )a vs (q 2 1 ), where the difference to a vv (q 2 ) solely originates from the polynomial P v/s (t).Finally, for the isoscalar-isoscalar contribution we set in analogy to Eq. (2.32).For the actual application, we further need to determine the overall sign and the relative size of the various isospin components, a point to which we will return in Sec.2.4.
The remaining parts enter in the context of the K * contribution, the third diagram in Fig. 3. To evaluate this diagram, we analyze the double discontinuity arising from the two-pion and πK ≃ K * intermediate states, in analogy to 2π and 3π singularities in the context of ρ-ω mixing in Ref. [95].We obtain the result where the overall sign is related to the one in c vv (q 2 1 , q 2 2 ) and the same expression holds for c sK * (q 2 1 , q 2 2 ) with a vv (q 2 1 ) → a s (q 2 1 ).The contribution from the ρ can be derived from the double discontinuity, leading to which coincides with a vv (q 2 ), a vs (q 2 ) in the limit P v/s (t) = 1.The coefficients of the ω and ϕ contributions are parameterized in analogy to the BMS ansatz (2.8), with the propagators replaced by a dispersively improved Breit-Wigner prescription as defined in Eq. (2.32).The remaining free parameters are ϵ K * and β K * , where the latter is necessary to accommodate the slope parameter (2.9).Factorization-breaking contributions to M v/s (t, q 2 ) are generated when including lefthand cuts (LHCs), which are expanded in a polynomial in Eq. (2.27).Neglecting isospinbreaking contributions, the leading LHCs to K L → π + π − γ * arise from a pion pole (for M v (t, q 2 )) and two-pion intermediate states (for M s (t, q 2 )).Since the former involves a K L → ππ vertex, this contribution will be CP violating.The two-pion LHC in M s (t, q 2 ) conserves CP , but involves an anomalous γ3π vertex (and only appears in the isoscalar amplitude).Accordingly, the role of LHCs beyond their polynomial approximation should be negligible, especially in the limited kinematic domain accessible in the kaon decay.1. Fit to the K L → π + π − γ spectrum from Ref. [84], see main texts for more details.

Fits to
As anticipated in Sec.2.2, experiments have observed evidence for an E * γ -dependent formfactor modification in contrast to a pure-M 1 DE amplitude [83][84][85].This enables us to determine the polynomial in Eq. (2.27).Accordingly, we perform a linear-polynomial fit to the spectrum data from Refs.[84] by (2.37) In the singly-virtual decay, it is not possible to separate all isospin components, instead, the spectrum is sensitive to the parameter combinations In particular, it is useful to rewrite the amplitude as to isolate normalization and t dependence in the fit.Ideally, one should then perform a combined analysis of all data for the spectra of the leptonic decays K L → ℓ + ℓ − γ and the hadronic channel K L → π + π − γ, in order to disentangle the linear term in Eq. (2.39) from the K * contribution.Unfortunately, the original data for these spectra are rarely available: the best measurements of K L → e + e − γ [76] and 77] are only provided in terms of fits using the two models introduced in Sec.2.1, but the original data for the spectrum are not included in the publications.Similarly, we could only obtain the K L → π + π − γ data from Ref. [84] by digitizing Fig. 4 therein (we checked that a fit to the digitized data reproduces the fit parameters given in the paper), while the spectra from Refs.[83,85] are not accessible at all.Given the high degree of correlation between the linear polynomial and the K * piece, the constraints from the dilepton data are critical, in such a way that we do need to rely on the value of η K * inferred from the global average of α K * [6].Accordingly, we will also use the slope determination via the BMS parameterization (2.9), to ensure a consistent set of input quantities.data from Refs.[76,77,83,85] were available.With η K * constrained to its BMS value, we then fit the remaining parameters N ππ and bππ to the spectrum from Ref. [84], first, with normalization in arbitrary units.In a final step, N ππ is rescaled to reproduce the total DE decay rate [6], for which our fit produces an uncertainty estimate consistent with Ref. [84] after taking into account the correlations between the fit parameters.The fit results are shown in Table 1 and Fig. 4, and will serve as basis for the error estimate of the dispersive form factor in Sec. 5. 6

Vector meson dominance
The K L → π + π − γ channel only determines the sum of the isovector and isoscalar parameters in the form of Eq. (2.38) as they enter in the singly-virtual form factor, but not the relative weights between the different isospin contributions.In such a situation, the vector-meson-dominance (VMD) picture often provides useful information to constrain the relative strength of different isospin components in the doubly-virtual form factor as well, see, e.g., Refs.[72,101,102].
The SU (3) pseudoscalar matrix |ϕ⟩ transforms under CP as CP |ϕ⟩ = −|ϕ⟩ T .The phase in front is fixed by the CP quantum numbers of the neutral states, e.g., from CP |π 0 ⟩ = −|π 0 ⟩.Thereafter, the K L state is defined as , in line with the prescription applied in Ref. [5].After identifying the phase convention in this way, we can write down the simplest trace that correctly projects out the K L state, where |V ⟩ is the vector-meson nonet matrix and we assume ideal mixing for ω and ϕ.Adding the relative factors from the V γ couplings [72,101,103], the relative weights between the different contributions to the form-factor normalization read c vv (0, 0) : c vs (0, 0) : c ss (0, 0) = 1 : − 1 3 : 5 9 . (2.42) The same calculation also predicts the relative size of ω and ϕ contributions to a s (q 2 ) and a ss (q 2 1 , q 2 2 ), in particular Similarly, the trace Tr {λ 6 , |V ⟩}|V ⟩ , (2.44) in combination with the V γ couplings, predicts the relative size of the K * -γ transitions These ratios differ from the ones assumed in the BMS model (2.8), 1 : 1/9 : 2/9, which are obtained when removing the singlet from |V ⟩ in Eq. (2.44), as done in Ref. [74].The two variants give an identical decomposition into isoscalar and isovector contributions, but the relative size of ω and ϕ differs.Although the resulting change in the K * parameterization (2.35) scales with M ϕ − M ω , we observe that the effect is numerically significant.Accordingly, since we need to rely on ϵ K * from the BMS fits, we employ the same decomposition of ω and ϕ, leading to the form anticipated in Eq. (2.35).Finally, while the normalizations c vK * (0, 0) = c sK * (0, 0) vanish by construction, we can consider the ratio which is constrained in the dispersive approach in terms of a vv (q 2 ) and a s (q 2 ).Here, VMD predicts the same relative weights of ρ, ω, and ϕ as in Eq. (2.45), which does not agree with c ϕ = 0 in a s (q 2 ), i.e., in a similar way as in the comparison to the BMS parameterization we observe an ambiguity in the separation of ω and ϕ contributions.The relative normalization  Ultimately, the impact of these VMD assumptions on the final uncertainty proves relatively minor, for the following reasons: the singly-virtual form factor is probed directly in terms of K L → ℓ + ℓ − γ and K L → π + π − γ, so that VMD arguments only enter in the doubly-virtual direction.Second, the ambiguities related to ω and ϕ only affect the internal decomposition of the isoscalar contribution, and, e.g., the three variants c ϕ = ±2c ω , c ϕ = 0 for the K * affect Re A µ (M 2 K ) at a level below 0.1.In general, we find that the result for the reduced amplitude Re A ℓ (M 2 K ) is rather robust against variation of the relative VMD weights, see Sec. 5, once the constraints from the form factor normalization (2.2) and slope (2.9) are imposed.As discussed in Sec. 4, only one of the relative normalizations (2.42) needs to be provided to determine all parameters, so that choosing either one of them serves as a way to assess the impact of the VMD assumptions, and the corresponding variation will be included in the uncertainty estimate for the dispersive part of the form factor quoted in Sec. 5.

Matching to asymptotic constraints
The asymptotic form of c(q 2 1 , q 2 2 ) was constrained in Ref. [12] using arguments from a partonic calculation, which we will incorporate here in our dispersive representation.The key result reads with tree-level Wilson coefficients normalized as The loop functions originate from diagrams with cand u-quark loops as depicted in Fig. 5, where either both photons couple in the loop, I(q 2 1 , q 2 2 ), or just one, T (q 2 i ), and result from the difference of the two flavors. 7However, we note that Eq. (3.1) as in Ref. [12] is valid for deeply-virtual kinematics and up to power-suppressed terms of O(ξ), ξ = m 2 u /m 2 c .In order to establish a matching between the low-energy dispersive contributions and the asymptotic constraints, we start from the full expressions including the dependence on m u , and separate the low-and high-energy contributions by introducing a cutoff in a dispersive representation of the respective loop integrals, see App.B for details.

Two-point function
The starting point for T (q 2 ) is given by the usual vacuum polarization from a quark loop, i.e., It is clear that T (r) does not contribute to the form factor normalization, a property that should be preserved when removing the low-energy contributions.Moreover, for q 2 → ∞ the difference of cand u-quark integrals decreases as 1/q 2 , which should also be reflected by the final representation.
The definition of T (q 2 ) in Ref. [12] follows when dropping the log ξ term and the dependence on ξ in the integral (3.2), i.e., T (r) [12] = 2 This loop function has a good asymptotic behavior for r → ∞ but displays a sensitivity to infrared scales.In particular, T (r) is now divergent for r → 0, as can be made explicit by rewriting Eq. (3.3) as For the asymptotic matching, we thus need to impose integration cutoffs in the dispersive representation for T (q 2 ) where σ q (s) = 1 − 4m 2 q /s, see Eq. (B.4).Here, the integration cutoffs s u and s c are related to ensure the correct asymptotic behavior, e.g., for the canonical choice s c = 4m 2 c one has s u = e 5/3 m 2 c ≃ 5.29m 2 c .

Three-point function
As detailed in App.B, the full loop function can be reconstructed from Refs.[105,106].The resulting expression is well-defined for all virtualities, thus correctly reproduces the partonic contributions to the real-photon kinematics obtained in Refs.[107][108][109], i.e., it simplifies to where r 3 = M 2 K /m 2 c .We also reproduce the result of Ref. [12], by taking r 3 = 0 and combining the uand c-quark contributions, (3.9) Although the limit ξ → 0 is finite, the final expression in Eq. (3.9) is not well-defined for all virtualities, and we again need to isolate the low-and high-energy contributions by imposing a cutoff in the dispersive representation.
As a first step, we analyzed the discontinuities of I c (q 2 1 , q 2 2 ), which shows that even for general r 3 the expression reduces to a simple scalar loop integral with C 0 as defined in Eq. (B.7).Moreover, the analysis in App.B shows that upon introducing a cutoff s u in the dispersion relation, no infrared singularities occur and the limit m u → 0 can be taken, in such a way that i.e., no non-trivial contributions remain from the u-quark loop.For the c-quark loop, setting M K = 0 permits a very simple double-spectral representation of the form . (3.12) A generalization of the de-facto single-variable dispersion relation in the second line to the case M K > 0 is given in Eq. (B.12), while a full double-spectral representation requires a contour deformation according to Eq. (B.20).We observed that the numerical difference between the M K = 0 simplified form and the full representation (B.20) is negligible for the current application.
4 Final representation To derive the final representation we first need to establish the sign convention for the K L → γ * γ * form factor c(q 2 1 , q 2 2 ).A detailed discussion of the sign of the SD contribution relative to the LD part is provided in Ref. [12], under the assumption that the dominant contribution to A µν [K L → γγ] arises from the π 0 pole c(0, 0) indeed close to the experimental value |c(0, 0)| = 3.389(14) × 10 −9 GeV −1 .Accordingly, it seems safe to assume that the signs of A[K L → γγ] and A[K L → π 0 → γγ] coincide. 8The sign of G 8 cannot be extracted from experiment, but it can be determined by considering hadronic matrix elements of four-quark operators in the ∆S = 1 Lagrangian, ⟨π 0 |H W |K L ⟩, in a factorization assumption, which leads to G 8 < 0 [30,111].This conclusion agrees with large-N c considerations [23,24], and we will therefore adopt c(0, 0) = −3.389(14)× 10 −9 GeV −1 (4.2) as boundary condition.The final representation is then decomposed as where the asymptotic contribution follows as given in Sec. 3 with s c = 4m 2 c , s u = e 5/3 m 2 c , and Wilson coefficients C i (µ) including the QCD renormalization group (RG) running from Ref. [11] to resum the large logarithms (α s log M W /µ) n and α s (α s log M W /µ) n to all orders in α s .The scale is set to µ 2 = |q 2 1 + q 2 2 |/2, and the C i (µ) are kept constant below µ cut = 2 GeV.Since the RG corrections to the linear combination C 2 + 3C 1 happen to be large, canceling a significant part of the tree-level result, it is worthwhile to also consider the next-to-leading-order RG corrections from Ref. [11] to obtain a more realistic estimate of the perturbative uncertainty, despite the dependence on the scheme for γ 5 .We choose the result in the 't Hooft-Veltman scheme [112] as our central value, and the variation to the leading-order RG and naive dimensional regularization as an estimate of the uncertainty.c asym (0, 0) then contributes only 2.1% to the form factor normalization, while the correction to the slope is negligible.
In our final representation, the dispersive part is decomposed according to 8 At leading order in SU (3) ChPT the π 0 and η contributions cancel (assuming the Gell-Mann-Okubo mass formula), but any realistic η-η ′ mixing scheme [23,110] implies a destructive interference between η, η ′ and thus π 0 dominance.The sign of the LD contribution is also discussed in Ref. [21].
To determine the number of independent parameters, we first observe that Accordingly, the quantities in Eq. (2.38) (with results collected in Table 1), the total normalization (4.2), and one of the relative VMD weights given in Eq. ( 2.42) allow one to separate the individual terms in and choosing either one of the three relative VMD ratios serves as an estimate of the theoretical uncertainty.The product a s (0)a vs (0) is resolved by employing Eq. (2.47).
Next, we observe that a vv (q 2 )a vv (0) + a vs (q 2 )a s (0) = a ρ (q 2 )N ππ + a ρ (q 2 )q 2 + āρ b ππ , (4.11) which, together with the normalizations a vv (0), a vs (0), a s (0), and a ss (0, 0), is sufficient to determine the singly-virtual dependence in the first line of Eq. (4.5), since a s (q 2 ) and a ss (q 2 1 , q 2 2 ) do not involve further free parameters.In contrast, a vv (q 2 ) and a vs (q 2 ) are more difficult to disentangle.The definitions of N ππ , b ππ in Eq. (2.38) and of a vv (0), a vs (0) in Eq. (4.7) imply the singular system of equations so that not all parameters N s/v , α s/v can be determined.For instance, a given solution remains invariant under Under this transformation, the resulting form for a vv (q 2 ) changes by

.14)
The shift therefore vanishes at q 2 = 0 (by construction), for large q 2 (since then the integrand becomes odd under x ↔ y), and in the limit of a narrow resonance (since then the spectral functions force x − y to zero).In practice, we indeed see that the remaining ambiguity in a vv (q 2 ) and a vs (q 2 ), which, due to Eq. (4.11), can only affect the doubly-virtual form factor, is completely negligible for reasonably chosen parameter space, so that for all practical purposes the entire functional form of c disp (q 2 1 , q 2 2 ) is predicted.The scale Λ ππ in Eq. (4.7), beyond which the linear polynomials are continued as constants, is chosen at Λ ππ = M K .With this assignment, our representation saturates the normalization sum rule predicted by the VMD ratios (2.42) at a level around 115%.In practice, to estimate uncertainties we impose the sum-rule normalization exactly and choose one of the three relative weights in Eq. (2.42) as predicted from VMD.The average of the three variants is taken as central value, and the spread as a measure of the systematic uncertainty.Similarly, we take c ϕ = 2c ω as central value for the K * piece, see Sec. 2.4, but include the variation to c ϕ = −2c ω and c ϕ = 0 in the error estimate.

Reduced amplitude
The reduced amplitude A ℓ (M 2 K ) in the prediction for the normalized branching fraction as defined in Eq. (1.5) is typically expressed as where q is the momentum of the kaon and c(q 2 1 , q 2 2 ) = c(q 2 1 , q 2 2 )/c(0, 0) is the normalized transition form factor.Therefore, we need to perform the integral for the final representation (4.3) to obtain A ℓ (M 2 K ).For the isospin part, we can evaluate its contributions following Refs.[29,113], e.g., for the isovector-isovector piece, we can write it as where K(x, y) is the integration kernel.For the K * contribution, we need to account for a new topology K(x, y, z) in addition to performing the integral over the pertinent spectral densities.The integration kernels are spelled out in App.C, and their reductions and numerical evaluations are performed using FeynCalc [114][115][116], LoopTools [117], and Package-X [118].The theoretical uncertainty of the dispersive representation is taken into account as follows: first, the matching threshold √ s m is varied in the range 2.0(2) GeV; fit uncertainties are taken into account by varying the parameters of Table 1 in their respective ranges (but the effect turns out to be small); the dependence on VMD assumptions is tested using the different variants for the VMD weights in the normalization as given in Eq. (2.42), see Sec. 4, and for the relative size of ω and ϕ contributions in the K * piece, see Sec. 2.4; lastly, the transition point Λ ππ above which the polynomials (2.37) in the evaluation of a vv (q 2 ) and a vs (q 2 ) are kept constant is varied between M K and 1 GeV.The dispersive uncertainty is then defined as the maximum difference to the central result.
The asymptotic contributions from the three-point and two-point functions read x 2 K(x, x), where the definition of K(x) is also collected in App. C. In practice, we estimate the asymptotic contributions by means of the approximate expressions given in App.C to be able to include the running of the Wilson coefficients in a straightforward manner.We have crosschecked that the approximate formula (C.2) indeed reproduces the results of the exact calculations (5.3) at tree level, up to tiny corrections that are entirely negligible compared to the dominant perturbative uncertainties. 9Therefore, we vary the scheme and the order of perturbation theory to estimate the asymptotic uncertainty: it is defined as the maximum deviation of calculations either in the naive dimensional regularization scheme or at the leading-log approximation from the central result evaluated in the 't Hooft-Veltman scheme.This turns out to be a conservative choice as it amounts to a 3% violation of the form factor normalization, well beyond the experimental precision (4.2).With these procedures for uncertainty estimates, we find the following imaginary parts of the amplitude Beyond the γγ imaginary parts, our analysis also takes the additional ππγ and 3πγ intermediate states into account via the dispersively reconstructed spectral functions, but the comparison of Eqs.(5.4) and (5.5) reaffirms that the imaginary parts are entirely dominated by the γγ cut [22].For the total LD contribution to the real parts, we obtain Re A LD µ (M 2 K ) = −0.50disp + 0.34 asym = −0.16(21) disp (27) asym (17) exp [38] total , Re A LD e (M 2 K ) = 31.99disp − 0.31 asym = 31.68(59) disp (73) asym (27) exp [98] total , where the errors refer to the uncertainties in the dispersive and asymptotic part of the representation, respectively, as well as the experimental uncertainty in the slope (2.9).Matching to Eq. (1.7), this translates to the low-energy constants χ µ (µ) = 4.96 (38), χ e (µ) = 8.0(1.0), evaluated at µ = 0.77 GeV.Accordingly, extracting χ(µ) at one-loop order is not sufficient to obtain the expected lepton-flavor-universal result [12,38,119], reflecting the impact of higher chiral orders [27,28].A similar effect has been observed for η, η ′ → ℓ + ℓ − [113], and the comparison to χ(µ) = 2.69(10) from π 0 → e + e − [29] highlights further that oneloop ChPT is not sufficient for an accurate phenomenology of pseudoscalar dilepton decays.
Adding the SD contribution (1.8), the complete real parts are 10 leading to the branching fractions (5.9) The final result for Re A µ (M 2 K ) can then be compared to experiment, Eq. (1.10), where we used Eq.(5.4) for Im A µ .This leaves room for a BSM component depending on the sign in Eq. (5.10).Choosing the negative value, our prediction thus displays a mild tension at the level of 1.7σ with experiment.For the electron mode, the experimental value reads Re A exp e (M 2 K ) = ±30.5(13.0), in agreement with Eq. (5.8), but with errors too large to derive meaningful BSM constraints.

Comparison to previous work
It is instructive to compare our results for the LD contribution to previous evaluations, see Fig. 6 and Table 2. 11 The simplest model, a VMD form factor with mass scale fit to reproduce the slope (2.9), gives 10 If the sign of c(0, 0) in Eq. (4.2) were positive, the LD contributions would change to 27)asym( 17)exp [38] total , Re A LD e (M 2 K ) = 30.98disp + 0.31asym = 31.29(59) disp (73)asym( 27)exp [98] total , and the total real parts would become where the uncertainty is propagated from b K L .The corresponding diagonal form factor is compared to our result in Fig. 6, which shows that the main differences to the VMD model concern a faster decrease for small Q 2 and a negative asymptotic tail.While the dispersive part in Eq. (5.6) comes out similar to the VMD result, the asymptotic contribution entails a positive shift that brings the net result closer to the experimental value.Ultimately, this enhancement arises from logarithmic corrections that lead to an asymptotic behavior proportional to log Q 2 /Q 2 , a feature that is absent in both the VMD and DIP models, with asymptotic limits proportional to 1/Q 4 and 1/Q 2 , respectively.
propagating the uncertainties from α DIP .The comparison to our result is again illustrated in Fig. 6, which shows that the higher central value of Re more than compensating for the missing log Q 2 in the asymptotic behavior and thereby even reversing the sign of the resulting value for Re A µ (M 2 K ).The corresponding low-energy constant χ DIP µ (µ) = 5.70(28) comes out close to χ [12] µ (µ) = 5.83(1.0)th (0.15) exp .This central value was obtained in Ref. [12] from a DIP parameterization with α DIP = −1.611(44) and an expansion in M K and m µ . 12Using the same value for α DIP , but keeping the full loop integral, we find χ µ (µ) = 5.42 (15) for this DIP contribution.The uncertainty is dominated by an estimate of high-energy contributions by varying a phenomenological ansatz within the range allowed by the partonic expression (3.1). 13In this work, we improved the matching to asymptotic constraints by separating the low-and high-energy parts of the loop functions via dispersion integrals, to the effect that the main uncertainty now arises from the RG corrections to the Wilson coefficients.However, the range for the high-energy contributions considered in Ref. [12] does suggest a correction to the DIP solution in the direction of the VMD form factor, which is consistent with our findings and would move the value of the low-energy constant closer to Eq. (5.7).
6 Consequences for physics beyond the Standard Model Re where the normalization takes into account the relative sign compared to the LD contribution in the same conventions as in Sec. 4 and App.A (similar expressions are available in the literature, see, e.g., Refs.[121,122]).The axial-vector coefficient maps onto the Gilman-Wise basis [123][124][125] , and the tree-level matching to SMEFT coefficients [126,127] reads with flavor indices prst = ℓℓ21.

Modified Z couplings
A special case is given by the test of modified Z couplings, defined by the effective Lagrangian [128] which maps to The SM contribution to Z ds , corresponding to the last two diagrams in Fig. Limits on Re Z ds can also be extracted from which is sensitive to both the real and imaginary part of Z ds [128,130,131] Evaluating the SM contribution with the input quantities summarized in Ref. [131] (relying on Refs.[132][133][134][135][136][137][138]), but with the CKM parameters as collected in App.A, we find Br[K + → π + ν ν]| SM = 8.2(5) × 10 −11 , in good agreement with Ref. [139].Assuming Im Z BSM ds = 0 and choosing again the sign that ensures agreement with the SM, we obtain (6.11) leading to a limit almost identical to Eq. (6.8).This observation nicely illustrates the complementarity of K L → µ + µ − and K + → π + ν ν, since a combined analysis of the two processes can be used to disentangle BSM contributions to the real and imaginary part of Z ds at a similar level of precision.In view of the projected experimental advances for K + → π + ν ν at NA62 and HIKE, this motivates commensurate efforts to improve the SM prediction for K L → µ + µ − as well.
-25 - In this work, we analyzed the K L → γ * γ * form factor in dispersion theory.To predict its low-energy part, we made use of data for K L → γγ and to determine normalization and slope, respectively, but, crucially, also for K L → π + π − γ, to constrain the 2π cuts in the spectral function for an isovector photon.Together with minimal narrow-width parameterizations for the isoscalar spectral functions as well as the K * ≃ Kπ contributions, we found that the resulting form factor, albeit mainly constrained by singly-virtual data, leaves remarkably little flexibility for doubly-virtual kinematics either.The representation for the full form factor was then supplemented by an asymptotic contribution, which derives from a dispersive formulation of the partonic calculation of K L → γ * γ * , to be able to separate low-and high-energy contribution by a suitable cutoff in the integrals.We observed that besides uncertainties in the dispersive representation and from the experimental value of the slope parameter, another important source of uncertainty arises from the perturbative running of the Wilson coefficients, which we studied including the resummation of subleading logarithms.The central results are given in Eq. (5.6) for the real part of the long-distance amplitudes Re A ℓ (M 2 K ), and in Eq. (5.9) for the corresponding branching fractions.The sum of the long-distance amplitude for K L → µ + µ − and the short-distance contribution in the SM agrees with experiment at the level of 1.7σ, and the improved uncertainty in Re A µ (M 2 K ) allows us to formulate more stringent limits on physics beyond the SM.The uncertainty in the prediction of the long-distance contribution could be further improved if better data for the K L → π + π − γ and K L → ℓ + ℓ − γ spectra were available, this would allow one to reduce the dispersive uncertainty and the precision of the form-factor slope.In this context, we emphasize that it is unfortunate that the original data for the respective spectra from Refs.[76,77,83,85] are not accessible, as, especially for the dilepton data, we are thus forced to rely on the quoted parameters for fits using a specific model for the form factor.While we tried to minimize the impact of the corresponding assumptions by constructing our dispersive parameterization in a way consistent with these fits, they could have been avoided altogether if the data for the spectra themselves had been preserved.
The asymptotic uncertainty could potentially be improved using the interplay with lattice QCD [40][41][42][43]140], which could also help unambiguously determine the relative sign of long-and short-distance contributions.Further, it would be very helpful if lattice QCD could provide information on the relative weights between normalizations of different isospin currents, as this would allow one to test the vector-meson-dominance assumptions currently needed for the prediction of the doubly-virtual form factor. Already with the improvements implemented in this work, the final uncertainty in Eq. (5.8) is getting close to experiment (5.10), so that with future advances in theory it should become possible to match the current experimental precision, especially, if improved input for the K L → π + π − γ and K L → ℓ + ℓ − γ spectra became available.As we argued in Sec.6, the corresponding constraints on physics beyond the Standard Model are complementary to the ones obtained from K → πν ν, motivating further efforts to improve the SM prediction for K L → µ + µ − .and recast into the dispersive representation In the matching to short-distance constraints, we use a variant of Eq. (B.2) in which m u → 0 and lower cutoffs s q applied in the dispersion integrals, to separate the parts already accounted for via the low-energy contributions.This leads to written in a form that applies for space-like q 2 and expanded for asymptotic values in the last step.To ensure that the entire contribution drops as 1/q 2 , s u and s c therefore cannot be varied independently.In particular, one has s u > s c , e.g., for s c = 4m 2 c one needs to set s u = e 5/3 m 2 c .

B.2 Three-point function
In the literature, the three-point function is represented in the form [105,106] where r 1 = q 2 1 /m 2 q , r 2 = q 2 2 /m 2 q , r 3 = M 2 K /m 2 q .Starting from this parameterization, we calculated the discontinuity in r 1 , which takes a very simple form disc r 1 I q (r 1 , r 2 , r Such a form can again be used to isolate low-energy contributions, i.e., by imposing a lower cutoff s u , the limit m u → 0 can be taken, and only I u (r 1 , r 2 , r 3 ) = −1 remains.For the c-loop contribution, the single-variable dispersion relation suggests a form .
(B.12)Such a form, however, does not strictly suffice to isolate the low-energy contributions, since the second variable q 2 2 can still take arbitrary values.For this reason, it is advantageous to formulate a double-spectral representation and impose cutoffs in both virtualities.
As a first step, we consider the limit r 3 = 0 of the discontinuity (B.Unfortunately, a simple representation such as Eq.(B.16)only applies for r 3 < 0, otherwise, additional contributions from anomalous cuts need to be included [56,165,166].Such contributions arise because the branch points of λ 123 , r ± (r 1 , r 3 ) = ( √ r 1 ± √ r 3 ) 2 , (B. 19) can move through the cut between r + 2 and r − 2 onto the first sheet and require a deformation of the contour.In the case 0 < r 3 < 4, we find the representation Finally, in the case r 3 > 4, as relevant for q = u, the r± 2 become real again, but λ 1/2 (x, y, r 3 ) produces imaginary parts that correspond to the decay K → ūu.Such imaginary parts clearly need to be covered by the low-energy part of the representation and cannot be present in the asymptotic contribution.However, the scaling of the u-quark integral can already be read off from Eq. (B.15): introducing a cutoff s u in both integrations and evaluating the δ-function, the remaining integral becomes since s u prevents a potential infrared singularity for q 2 1 = q 2 2 = 0. Given that an analysis based on Eq. (B.20) yields the same scaling, we set I u (q 2 1 , q 2 2 , M 2 K ) = −1.

C Integration kernels
The integration kernels for the evaluation of the loop integral (5.1) are given as

Figure 6 .
Figure 6.The normalized diagonal form factor c(−Q 2 , −Q 2 ) (black dashed line with uncertainty band), in comparison to the VMD model and the DIP parameterization [38].

Table 2 .
Comparison of the LD amplitudes and LECs to VMD and DIP parameterizations, see main text for details.
[129]ves Z ds = λ t C 0 (x t )[129], and, together with the W box, determines the SD part of χ(µ) in Eq. (1.8).Choosing the sign in Eq.(5.11)in line with the SM, we obtain the constraint