Exact NLO matching and analyticity in b → sℓℓ

Exclusive rare decays mediated by b → sℓℓ transitions receive contributions from four-quark operators that cannot be naively expressed in terms of local form factors. Instead, one needs to calculate a matrix element of a bilocal operator. In certain kinematic regions, this bilocal operator obeys some type of Operator Product Expansion, with coefficients that can be calculated in perturbation theory. We review the formalism and, focusing on the dominant SM operators O1,2, we perform an improved calculation of the NLO matching for the leading dimension-three operators. This calculation is performed completely analytically in the two relevant mass scales (charm-quark mass mc and dilepton squared mass q2), and we pay particular attention to the analytic continuation in the complex q2 plane. This allows for the first time to study the analytic structure of the non-local form factors at NLO, and to calculate the OPE coefficients far below q2 = 0, say q2<~−10GeV2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {q}^2\underset{\sim }{<}-10{\mathrm{GeV}}^2 $$\end{document}. We also provide explicitly the contributions proportional to different charge factors, which obey separate dispersion relations.

Exclusive b → s decays such as B → K ( ) and B s → φ have been on the focus point of theorists and experimentalists for some time, due to the potential they provide for tests of the Standard Model (SM). While the interest for such decays dates back to the era of the B-factories (which provided some of the first measurements), a renewed interest has been triggered by the measurements at the LHC, most prominently the ones by the LHCb collaboration. Starting with the "P 5 Anomaly" [1,2], and followed by a larger pattern of "tensions" of different degrees in the landscape of angular and dilepton-mass-squared distributions in B (s) → {K ( ) , φ}µ + µ − modes [3,4], these measurements (in inseparable association with theoretical work) have guided the community during the LHC era. More precise experimental studies are part of the programs for the LHC upgrade [5] and Belle-II [6], and there is little doubt they will lead to new discoveries. The question is whether these discoveries will involve Beyond-the-SM (BSM) or QCD/hadronic physics. While this is subject to the personal inclination of the reader, both outcomes are truly interesting. The exclusive b → s decays belong to the class of "rare" FCNC transitions which are loop-, CKM-and GIM-suppressed in the SM. This leads to branching fractions of the order of 10 −6 , and which could be easily altered by BSM physics lifting any of such suppression mechanisms. However, it is increasingly evident that large deviations with respect to the SM are not present, and as such, rare decays are no longer smoking guns of BSM physics. Thus we need to test SM predictions more precisely. This is now possible due to the large statistics collected at the LHC (with more than 2K selected B → K µµ events in Run 1 by LHCb), but it also implies that theory predictions with uncertainties below ∼ 10% are necessary, with model dependence reduced to the minimum.
Theory predictions for B → M + − observables depend on non-perturbative hadronic matrix elements of two types: "local" and "non-local" form factors (e.g. [7]). Contributions to the amplitude from semileptonic ([sΓb][¯ Γ ]) or dipole ([sσ µν P R b]F µν ) operators are exactly factorizable and proportional to local form factors -matrix elements of local fermionic currents -to all orders in QCD (but to the leading order in QED effects). These local form factors are known relatively well and can be calculated with Light-Cone Sum Rules (LCSRs) or Lattice QCD (LQCD) methods, both agreeing well with each other [8][9][10][11][12][13]. On the contrary, contributions from four-quark operators such as [sγ µ P L c][cγ µ P L b] are proportional to non-local form factors, more precisely, the matrix elements of time-ordered products of a four-quark operator and an electromagnetic current. The calculation of these non-local form factors is highly non-trivial and relies inevitably on some type of operatorproduct expansion (OPE) [14][15][16]. In this way, the complicated non-local form factors can be written in terms of simpler hadronic matrix elements, multiplied by coefficients that can be determined though a perturbative matching calculation. These simpler hadronic matrix elements are either local form factors, or matrix elements of bi-local operators defined on the light-cone, which can be expressed in terms of meson light-cone distribution amplitudes.
The matching to the leading (dimension-three) operators in the OPE can be extracted from the perturbative partonic calculation of the matrix element, which has been known up to order α s (two loops) for some time [17][18][19], albeit not in full analytic form in the JHEP04(2020)012 two relevant variables: q 2 (the dilepton squared invariant mass) and m c (the charm-quark mass). Only recently the necessary analytic calculation of the two-loop master integrals involved has been achieved [20], and applied to the problem at hand [21].
We have repeated the full analytic two-loop calculation independently, and checked the results of ref. [20], which we confirm. The explicit and independent check of this calculation is the first result of this paper. But we have done the calculation in a way that lays out the analytic structure of the results more explicitly, and imposing an analytic continuation which is more convenient for the dispersive analysis (see refs. [7,22]). The results in this form allow us to study the branch cut discontinuities and compare them with the expectations derived from unitarity, as well as to test all the analytic singularities of the two-loop amplitude by explicitly checking a dispersion relation. This is the second result of this paper. Finally, the dispersion relation formalism is an important tool to extend consistently the calculations in the LCOPE region (negative q 2 ) to the physical region at q 2 > 0. Under certain simplifying assumptions, this dispersion relation can be separated in pieces multiplying difference quark charge factors. For that purpose the NLO contributions to the OPE coefficients must also be separated in this way, but this separation has not yet been given explicitly. We do give separate contributions to the OPE coefficients to be used in the separated dispersion relations, which is the third result of this paper.
We start in section 2 by reviewing the theoretical framework and fixing the conventions and the notation. In section 3 we give the details of the analytic NLO matching calculation. In section 4 we address the issue of the numerical evaluation of the NLO functions, which requires some care due to the presence of Generalized Polylogarithms (GPLs) up to weight four. We also compare our results with the ones in the literature, and provide explicit numerical results at various kinematic points in the LCOPE region. In section 5 we discuss the analytic properties of the results and prove the structure of singularities by means of a dispersion relation. We then explain how to separate the NLO matching coefficients into the two contributions proportional to different charge factors. We conclude in section 6. The various appendices include additional information on: A. The attached Supplementary material files which contain all our results in electronic form as well as codes for numerical evaluations; B. The list of the relevant Master Integrals that appear in the calculation of the two-loop diagrams; C. The list of different weights appearing in the GPLs in the results; and D. A few examples on fixing the integration constants that arise in the calculation of the two-loop Master Integrals.
2 Theoretical framework 2.1 Set-up: Weak Effective Theory and conventions B decay amplitudes are calculated within the Weak Effective Theory (WET) where the SM particles with EW-scale masses have been integrated out. The WET lagrangian then contains QCD and QED interactions, and a tower of higher dimensional local operators which is typically truncated at dimension six [23,24]. The part of the WET Lagrangian JHEP04(2020)012 which is relevant for the contributions discussed in this paper is: We use the following conventions: In our calculation of NLO corrections from O 1,2 , the scheme dependence of m b is a higher order effect. We will neglect the strange quark mass throughout the paper.

Local and non-local form factors in exclusive
To the leading non-trivial order in QED, the effective theory amplitude for the exclusive decayB → M + − , with M an undetermined meson (or hadronic state in general [13]), is given in terms of local and non-local form factors [7,13,25]: up to terms of O(α 2 ). Here q 2 is the invariant squared mass of the lepton pair and L µ i are leptonic currents, L µ V (A) ≡ū (q 1 )γ µ (γ 5 )v (q 2 ). In this amplitude we have neglected contributions from other local semileptonic and dipole operators that are not relevant in the SM, as well as higher order QED corrections, but it is exact in QCD. All nonperturbative effects are contained in the "local" and "non-local" form factors F (T )µ i and H µ , with This paper deals with the non-local form factors H µ (q, k), defined by the following matrix element: where j µ em = q Q qq γ µ q, with q = {u, d, s, c, b}. This corresponds to the matrix element of the non-local operator: which is the focus of the following discussion. JHEP04(2020)012

Operator Product Expansion for non-local form factors
A reliable calculation of H µ (q, k) is very important for phenomenology and a challenge for theory. At low hadronic recoil, q 2 ∼ m 2 b , the dx integral in eq. (2.6) is dominated by the region x ∼ 1/m b , and a local OPE exists for the operator K µ (q) [15,16]: where we have indicated the contribution of operators of dimension three (according to the counting in ref. [16]), and the ellipsis denotes contributions of operators of higher dimension d > 3, with OPE coefficients that are suppressed by This equation defines the OPE coefficients ∆C 7,9 . At large hadronic recoil, and below the on-shell branch cuts, q 2 0, the dx integral in eq. (2.6) is instead dominated by the region 1 x 2 ∼ 1/(4m 2 q − q 2 ), which allows for a light-cone OPE (LCOPE), where local operators with an arbitrary number of covariant derivatives along the relevant light-cone direction contribute at the same order [14]. The structure of the LCOPE coincides with the local OPE at dimension three, and therefore eq. (2.7) is also true at q 2 0. The power corrections are, however, different. Power corrections to both OPE expansions have been discussed in e.g. refs. [14,16,22].
Given eq. (2.7), the non-local form factors (2.5) are determined by the OPE coefficients and the local form factors: with the ellipsis denoting contributions from subleading terms in the (LC)OPE. Thus, the effect of the non-local contribution H µ in the amplitude (2.3) at this order in the OPE expansion can be absorbed into "effective" Wilson coefficients C eff 7,9 (q 2 ) = C 7,9 + ∆C 7,9 (q 2 ). These effective Wilson coefficients are scheme and scale independent. The same structure arises to all orders in QCD in the "factorization approximation", where all interactions between the charm loop and the constituents of the external mesons are neglected. However the OPE formalism beyond the leading order includes all non-factorizable contributions, which appear to be phenomenologically very relevant [26].

Structure of the OPE matching calculation
The OPE coefficients ∆C 7,9 (q 2 ) are calculable order by order in perturbation theory through a matching calculation. The easiest way to perform this matching is to equate the matrix elements of partonic states at each order in α s : We shall refer to the matrix element M µ (q) in the left-hand side as the "QCD amplitude" and the one in the right-hand side M µ OPE (q) as the "OPE amplitude". A perturbative calculation of the QCD amplitude leads to an expression of the form: (2.10)

JHEP04(2020)012
In ref. [17] the functions F (7,9) i (q 2 ) were calculated at low q 2 and the results were represented as expansions in the small parameters q 2 /m 2 b , z ≡ m 2 c /m 2 b and q 2 /(4m 2 c ). In ref. [18] the functions F (7,9) i (q 2 ) were calculated for the high q 2 range and the results were given as an expansion in z. In section 3 we describe the calculation of these NLO functions F (7,9) i (q 2 ) in a fully analytic form for z and q 2 . The full results are discussed in section 3.7.

Analytic structure and dispersion relations
In order to discuss the analytic structure of the non-local form factors, it is convenient to perform a Lorentz decomposition and focus on invariant functions: where η µ λ are a set of orthogonal Lorentz vectors depending on q and k and H λ (q 2 ) are a set of invariant non-local form factors (see e.g. ref. [7]).
Once the non-local matrix elements H λ (q 2 ) are known in the OPE regions of the q 2 plane, it remains to use this information to extrapolate the results to the physical regions of interest, within the range 0 ≤ q 2 ≤ (M B − M M ) 2 . For this we need some information about the properties of the functions H λ (q 2 ) in the complex q 2 plane. The most important of such properties is the analytic structure (the structure of their analytic singularities), that is, the presence of poles and branch cuts. Assuming the principle of maximum analyticity, these singularities are fully determined by the on-shell cuts of the matrix elements (see e.g. ref. [7]).
The first thing to note is that, independently of the value of q 2 , the functions H λ (q 2 ) are complex-valued due to on-shell intermediate states in the p 2 channel, e.g. B → DD s → M λ γ * . The singularity structure associated with the variable q 2 will then apply separately to the real and imaginary parts of H λ (q 2 ): H (re) λ (q 2 ) and H (im) λ (q 2 ). Each of these two functions are then real for q 2 < 0, but develop imaginary parts due to on-shell states in the q 2 channel, for q 2 > 0. All these on-shell states must have the (QCD-conserved) quantum numbers of the e.m. current, which means that (in full QCD) they are necessarily multiparticle states. Therefore the singularities are branch cuts, one for each multiparticle state: B → M λ X 1−− → M λ γ * , with X 1−− = {ππ, πππ, KK, · · · , DD, DD * , · · · }. Each of these branch cuts starts at its corresponding threshold s th = {4m 2 π , 9m 2 π , 4m 2 K , · · · , 4m 2 D , (m D + m D * ) 2 , · · · }.
Given the analytic structure of the functions H λ (q 2 ), one can write a dispersion relation to relate the values of these functions at specific points to an integral over the branch-cut discontinuity [22]: is the discontinuity along the cut (the spectral function). The spectral function ρ λ (t) may, in certain approximations, contain poles below the multiparticle threshold, and thus JHEP04(2020)012 in such cases the parameter s th is assumed to lie below such poles. The subtraction at q 2 0 is implemented to ensure the convergence of the dispersion integral [22]. While this dispersion relation is completely general, we assume that q 2 0 is within the OPE region (thus H λ (q 2 0 ) = H OPE λ (q 2 0 )), and q 2 can be on the physical range, and thus the i prescription in the denominator is chosen such that for (real) q 2 > s th , the pole in the integrand is above the real axis. This prescription can be ignored if q 2 is away from the branch cut.
One can now separate the different contributions to the e.m. current in eq. (2.5), and write three different dispersion relations for H λ,sb , H λ,c and H λ,ud [22]. These three dispersion relations are equivalent to eq. (2.23), but with two qualifications: (1) the spectral densities also depend on the channel, ρ λ,sb , ρ λ,c and ρ λ,ud , and (2) the OPE functions The reason that the terms with Q s and Q b are not separated is because they are not separately gauge invariant (see section 3.3), while the terms with Q u and Q d do not receive contributions from the two-loop matching corrections discussed in this paper, and will also depend on the charge of the decaying B meson. The explicit separation into terms with different charge factors Q s/b and Q c is one of the results in this paper that was not available before. The two-loop contributions to H OPE λ,sb (q 2 0 ) and H OPE λ,c (q 2 0 ) will come respectively from diagrams {a, b}, and {c, d, e} in figure 2. Other contributions from CKM-suppressed operators (with u, d, s loops) will contribute to H OPE λ,ud (q 2 0 ) and H OPE λ,sb (q 2 0 ). These corrections are simpler than the ones discussed in this paper (since they contain one fewer mass scale) and can be found in analytical form elsewhere [27].
Up to this point the discussion is rigorous and exact, relying only on maximum analyticity and unitarity. The separation into different charge factors has been performed to implement a simplifying assumption when modelling the spectral densities, based on OZI suppression [7,22]. Up to OZI-suppressed effects, the QCD spectral densities ρ λ,sb , ρ λ,c and ρ λ,ud receive separable contributions from intermediate states {φ, KK, . . . }, {J/ψ, ψ(2S), DD . . . } and {ρ, ω, ππ, . . . }, respectively [7,22]. Therefore the dispersion relation can be divided into three separate ones [22]: For consistency with the adopted approximation we have assumed that the resonances below the multi-particle thresholds in each channel are stable, and indicated only these poles in the spectral densities. The ellipses denote the subsequent continuum contributions with open flavors (e.g. DD, D * D , · · · in ρ λ,c (t)). The flavor separation of the dispersion relations has some phenomenological advantages [7,22]. JHEP04(2020)012 NLO does not contribute to the NLO matching. Right: Contribution to f (9) NLO which is equal to f NLO . The contribution of this diagram to f (7) NLO vanishes (since f

OPE functions at NLO and cancellation of IR divergencies
The NLO functions h (7,9) NLO arise from the diagram in figure 1 (left). According to the matching equations (2.18), (2.19), only h (9) NLO is needed for the NLO matching. On the other hand, the contribution to the function f (9) NLO given in figure 1 NLO , since the LO matching expression ∆C 9,LO = f (9) LO ensures that the charm loop can be replaced by K µ OPE at this order in the perturbative expansion. Thus, the two contributions will cancel in the combination f This cancellation is important because these are the only two contributions which are IR divergent. As a result, the NLO contributions in eq. (2.21) are obtained by evaluating the five classes of diagrams in figure 2.

Two loop contributions to the QCD amplitude
The contribution to the QCD amplitude from any given set of Feynman diagrams in figure 2 can be written as Conservation of the e.m. current implies that V µ (i) has the structure of eq. (2.10): which is a consequence of the Ward Identity to be checked from the calculation. In the calculation of V µ (i) , we use the EOM for the quark spinors (keeping m b = 0 but setting m s = 0 here) to remove all factors of / p and / q, and we set p 2 = m 2 b and (p − q) 2 = m 2 s → 0. At the end one finds that V µ (i) has the form: JHEP04(2020)012 where A (i) , B (i) and C (i) are scalar functions of m b , m c and q 2 . On dimensional grounds, From these coefficients one can read off the functions f (7,9) (i) (q 2 ) and check the Ward Identity. From A (i) and B (i) one has: and the Ward Identity is respected if and only if the coefficients C (i) satisfy: This condition applies to gauge-invariant combinations and not to single diagrams. We will detail which are the gauge-invariant combinations below. We evaluate scalar quantities A where the numbers n i are integers (positive or negative), with N i = 7 j=1 n i j , the objects P i are propagators (see below), and the indices {i 1 , . . . , i 7 } depend on the class. In addition, d = 4 − 2 , andμ 2 ≡ µ 2 e γ E /4π, with µ the MS scale. Our choice of momentum routings fixes the first five propagators in each class, and the other two are chosen to be linear in loop momenta and such that the seven propagators form a linearly-independent set. The JHEP04(2020)012 complete set of propagators needed is: and the scalar integrals for each class are: j[a; n 2 , n 3 , n 4 , n 5 , n 8 , n 7 , n 9 ] , j[d; n 1 , n 2 , n 12 , n 4 , n 11 , n 6 , n 7 ] , j[b; n 2 , n 3 , n 4 , n 10 , n 11 , n 7 , n 9 ] , j[e; n 1 , n 2 , n 3 , n 4 , n 12 , n 7 , n 13 ] , is solved. In the following we describe the analytic calculation of the two-loop scalar integrals.

IBP reduction and Master Integrals
At this point, the result of each diagram is a function of many scalar integrals with many different tuples {n i 1 , . . . , n i 7 } in its class. We can now use integration-by-parts identities (IBPs) to reduce the set of scalar integrals appearing in each class to a small set of Master Integrals (MIs). For this purpose we use the Mathematica code LiteRed [28]. After reduction, the total number of two-loop MIs in each class is m i = {7, 9, 9, 15, 5} for i = {a, b, c, d, e}, respectively. These MIs are listed in appendix B, and collectively denoted by J i,k , with i = {a, b, c, d, e}, and k = 1 . . . m i for each i.
With the functions A (i) , B (i) , C (i) written in terms of MIs one can check the Ward Identity by verifying eq. (3.5), which holds analytically and explicitly in terms of the unevaluated MIs. This does not happen individually for each diagram, but for the following combinations: , and e 1 + e 2 + e 3 , according to the numberings in figure 2.
We now perform some simplifying operations on the master integrals. First, we express the integrands themselves in terms of the invariant variables on which the scalar integrals depend, which we choose to be For this purpose we note that there always exist two light-like vectors k 1,2 (k 2 1 = k 2 2 = 0) such that p = k 1 + k 2 and q = k 1 + s k 2 . Then p − q = (1 − s) k 2 , and the condition (p − q) 2 = 0 is automatically satisfied. In addition, k 1 · k 2 = m 2 b /2. Thus, expressing the integrands in terms of k 1,2 instead of p, q leads to the (dimensionless) scalar integrals j[i; {n i }] as explicit functions of (s, z).
Second, in order to be able to do a rational transformation to a canonical basis of master integrals (as explained below), for each set of diagrams (i) we make a change of JHEP04(2020)012 variables (s, z) → (x i , y i ), with x i = x i (s, z) and y i = y i (s, z) a set of functions that will be specified later. In terms of these new variables the dimensionless MIs are written as J i,k ( , x i , y i ).

Differential equations in canonical form and iterative solution
For each set of diagrams, we construct the system of differential equations: where a i,x , a i,y are m i ×m i matrices depending on , x and y. The derivatives of the MIs J i,k are performed by differentiating the integrands, which produce new scalar integrals, and then applying the IBP reduction again on these scalar integrals to express the derivatives ∂ x,y J i,k themselves in terms of the MIs J i,k . One can then read off the matrices a i,x and a i,y . A basis of Master Integrals is said to be "canonical" [29] if a x,y ( , x, y) = A x,y (x, y), with A x (x, y) and A y (x, y) two N × N matrices independent of . Given a canonical basis M , the differential equations have the form: Although not explicitly used in the following, we note that there is a matrixÃ( Once a canonical basis is found, the system of differential equations can be solved automatically order by order in . To keep the notation as simple as possible in this section, we will assume that all the master integrals in the canonical basis are regular in (if not, we redefine them by multiplying all of them with the same appropriate power of ). We then write the -expansion for the master integrals and the differential equations read: We first construct the general solution of the differential equation containing the derivative with respect to y. Using partial fraction decomposition, A y can be written in the form where A j y a set of constant matrices, and the quantities w j (x) are called the "x-dependent weights" (see appendix C). These differential equations can be solved iteratively due to the JHEP04(2020)012 structure of (3.13): etc., in terms of Generalized Polylogarithms (GPLs) [30], defined iteratively as [31] G(w 1 , . . . , w n ; 16) where 0 n denotes n consecutive zeroes. In each step of the iteration, integration constants (with respect to the y integration) are added, which however depend on the variable x; they are denoted as C n (x).
Using the fact that the GPLs in the above equations either tend to zero in the limit y → 0 or to log n x n! (when all n weights are zero), it is straightforward to derive ordinary differential equations for the C n (x) quantities, obtaining The matrix A x evaluated at y = 0 has, after partial fraction decomposition, the form where A k x is again a set of constant matrices, and the quantities w k are now constant weights (see appendix C). The solutions of the differential equations for C n (x) again are determined iteratively: where C n on the right-hand side are constants with respect to both variables. Thus, the problem of calculating the MIs is reduced to find a canonical basis and to fix the integration constants, which is a much more tractable challenge. In order to find a canonical basis for each set J i,k of MIs, we use the mathematica program CANONICA [32].

JHEP04(2020)012
This code is able to look for transformations that involve rational functions of the arguments. For this reason, the right set of variables (x i , y i ) must be found for each case before using this program. Starting from our original variables s = q 2 /m 2 b and z = m 2 c /m 2 b , we define, for each diagram set i, the variables x i and y i : In terms of these variables and with the help of CANONICA, we are able to find linear transformations such that the MIs M i,k constitute a canonical basis for each set i = {a, c, d, e}. For set b, the situation is somewhat more complicated: There is a linear transformation involving rational functions of the arguments x b and y b for the MIs J b,1−6 and this six-dimensional block can be treated in a straightforward way, but the complete nine-dimensional problem contains complicated square roots of these variables in the transformation matrix to the canonical basis and in the matrices A x and A y which define the differential equations in this basis. Similar as after eq. (4.46) of ref. [20], we introduced the variables t b and v b to rationalize these roots: For this reason the results for the MIs J b,7 , J b,8 , and J b,9 involve GPLs with arguments t b and/or v b . We stress that the chosen variables x i have the properties that they tend to zero when z goes to infinity. Similarly, the variables y i (as well as t b and v b ) go to zero for s → 0 (when i = b, c, d, e) and for s → 1 (when i = a), independently of the value of z. In these limits, the functions G(. . . ; x i ), G(. . . ; y i ), G(. . . ; t b ) and G(. . . ; v b ) can be expanded in a straightforward way for the small values of x i , y i , t b and v b , respectively. This turns out to be very useful when fixing the integration constants in the following section, because we will heavily make use of the asymptotic properties of the originals integrals J i,k in the limit where x i and/or y i , t b , v b go to zero.

Fixing integration constants and analytic continuation
Once the canonical basis is found and the general solution of the differential equations in this basis is constructed, we have to fix the integration constants. To this end we transform JHEP04(2020)012 in a first step the MIs back to the original basis by making use of the transformation matrices T i (i.e. eq. (3.21)). The constants are then determined by either computing the MIs J i,k in the various classes i at a particular kinematical point for which the calculation is simple, or by using asymptotic properties in the limit z → ∞. These properties follow in a straightforward way from the heavy mass expansion (HME) of a given integral [33].
We explain this in some detail for the nine MIs in class c: it turns out that only the integral J c,1 , which is simply a product of two one-loop tadpole integrals, has to be calculated explicitly. In the limit for large m c (m c m b ) the other eight integrals can be naively Taylor expanded in the external momenta and in m b . Note that in the present situation the only subdiagrams in the sense of the HME are just the full diagrams (i.e. the full MIs) and therefore the naive Taylor expansion is justified. The leading power n in the m c -expansion of a given integral J is then identical to the mass dimension of the integral, where the mass dimension is an even integer; the structure of J is where K is a constant prefactor and P is a polynomial of the indicated arguments. The GPLs in the general solution for the MIs (from the differential equations) can be easily expanded for large z and small s in class c. Very often, the expanded solution for a given integral contains higher powers in m c than that determined from the HME argumentation. The requirement that these terms are absent allows to determine some of the integration constants. From the HME structure it is also clear that only even powers of m c can be present; this fact fixes the remaining integration constants. It is worth emphasizing that all constants can be fixed by the explicit knowledge J c,1 in class c and the structure of the powers in m c . The explicit HME evaluation of the MIs is not even necessary.
For classes {b, d, e} the fixing of the integration constants is done in the same way as in class c: only a small number of simple one-loop integrals have to be calculated explicitly; again the GPLs in the results for the MIs (from the differential equations) can be easily expanded for large z and small s and all constants can be fixed. A few examples on the fixing of integration constants in classes c and e are given in appendix D.
We now turn to class a. Due to the variable y a = 1/ 1 − 4z/(1 − s), we need to use the behavior of the solutions of the MIs near s = 1 (not at s = 0 as in the other classes) and again for z → ∞. Apart from heavy mass expansion arguments (which are the same as in the other classes), we need to calculate directly the three integrals J a,1 , J a,4 and J a,5 (which all factorize into two one-loop integrals), in order to fix the integration constants. Among them, only J a,4 depends on s. The explicit result reads complex plane cut along the positive real s-axis, having a discontinuity on this axis. However, when expanding the GPLs in the solution of the differential equations for J a,4 around s = 1, we find a regular behavior, which is due to the fact that the solution in terms of GPLs with argument y a represents a different analytic continuation. In order to obtain an analytic continuation with the branch cut along the positive real s-axis (see sections 2.5 and 5) we need to consider the differential equations for the upper and the lower s-half planes separately. In particular, we have to fix the integration constants for the two pieces separately. In this way, the branch cuts in all classes appear along the positive real axis, starting at s = {0, 4z, 4}, depending on the class. These branch cuts will be analyzed in detail in section 5.
Our final results for all MIs in the Feynman region have been checked numerically using Sector Decomposition as implemented in SecDec [36,37].

Counterterm contributions
For the renormalization we will follow closely ref. [17], and therefore we prefer to stick to the notation of that paper within this section: (3.26) The set of operators O 1 -O 10 is given in eq. (2) of ref. [17], while O 11 and O 12 are evanescent, that is, they vanish in d = 4 dimensions. Although there is certain freedom in the choice of the evanescent operators (e.g. one may add terms of order ), it is convenient to use the same definitions as in ref. [34] in order to combine our matrix elements with the Wilson coefficients calculated there: The renormalization constants δZ ij are written as  (3.30) The counterterm contributions to the functions F (7,9) i due to the mixing of O 1,2 into fourquark operators are denoted by F ct (7,9) i→4 quark , and are related to the one-loop matrix elements of four-quark operators by  (7,9) i→9 , and given by [17] F ct(7) for which the strong coupling renormalization constant Z gs is needed: In addition to the contributions from operator mixing, there is a contribution from the renormalization of the charm quark mass. This is taken into account by replacing m c with Z mc · m c in the one loop contributions given in eq. (2.11). Note that in this paper we are using the pole mass definition of m c , characterized by the renormalization constant We have checked that the sum of the divergent parts of all these counterterm contributions is identically opposite to that of the unrenormalized matrix elements, thus proving the cancellation of ultraviolet divergences.

JHEP04(2020)012
On the other hand, the finite part of the counterterm contributions, which we denote by F ct(j) i (i = 1, 2; j = 7, 9), contribute to the renormalized NLO functions F (j) i . Besides working out the exact results for the counterterm contributions F ct(j) i in terms of GPLs, we have also separated the different contributions proportional to the different charge factors Q s,c,b , since they renormalize the different contributions to F (j) i with different analytic structure. It turns out that the only contributions proportional to Q s,b to F ct(j) i come from the mixing O i → O 4 , specifically from the one-loop matrix element of O 4 with an s-or b-quark in the loop, and thus these contributions are easy to isolate. In the end, our results for the counterterm contributions are given by the sum of three pieces:

Results for renormalized matching coefficients at NLO
Collecting all the pieces, the final results for the matching coefficients ∆C 7,9 (q 2 ) in eq. (2.8) at NLO are given by where f LO (q 2 ) is given in eq. (2.11) and the renormalized NLO functions F (7,9) 1,2 (q 2 ) are the sum of the contributions from the two-loop diagrams a through e and the counterterm contributions: F (3.40) The complete analytic results for the functions F  The coefficients ∆C 7,9 (q 2 ) can also be split in the two different contributions ∆C (c) 7,9 (q 2 ) and ∆C (sb) 7,9 (q 2 ) proportional to the charge factors Q c and Q s,b respectively, and contributing to the functions H OPE λ,c (q 2 ) and H OPE λ,sb (q 2 ) discussed in section 2.5. For this separation we refer to section 5 below.

Numerical evaluation of GPLs
For the fast numerical evaluation of the GPLs we use the C++ ginac package [35] interfaced with Mathematica. In particular, we use the ginac multiple polylogarithm G, to evaluate the GPLs with unit argument and the last weight non-zero: G(w 1 , . . . , w n ; 1) , with w n = 0 . (4.1) When w n = 0, the GPL with arbitrary (non-zero) argument is obtained from the identity while the GPL with zero argument is zero. This part is implemented by the Mathematica interface. In order to evaluate the cases with w n = 0 we need to eliminate all the "trailing zeroes" in the GPLs, which refer to any string of consecutive zeroes at the end of the weight list, e.g., G(1, −2i, 0, 0 ; 3 + i) has two trailing zeroes. Reexpressing the GPLs in terms of new GPLs without trailing zeroes is also done by the Mathematica interface, recursively in the number of trailing zeroes, by means of the following formula: This provides a complete algorithm for the evaluation of any GPL. For convenience, we provide our C++/Mathematica bundle (with front-end package GPL.m) as an ancillary file supplementing this paper (see appendix A.1 for details). All our numerical results have also been reproduced using Maple, which includes a built-in function for GPLs. However, the evaluation within Maple is significantly slower that the one provided by GPL.m. In order to properly evaluate our expressions, we consider separately the GPLs with arguments x i or y i . For GPLs with argument x i , we numerically evaluate x i by adding a small negative imaginary part to z, typically of order 10 −12 . For GPLs with argument y i , on the contrary, we evaluate the x i dependent weights in the limit in which the small imaginary part on z tends to zero; the arguments y i are calculated by taking z real from the beginning and by adding a small positive/negative imaginary part (typically of order 10 −8 ) to s, when s lies on the real axis.

Numerical evaluation of NLO corrections and tests
Once the numerical evaluation of the GLPs has been addressed, the numerical evaluation of the NLO functions F j i (s, z) is relatively simple. We use the Mathematica package FFNLO.m, which is appended as Supplementary material to this paper (see appendix A.2 for details). This program makes a prior list of all the GPLs appearing in the functions to be evaluated, evaluates them only once using GPL.m, and then substitutes the values in the functions. [17] at low-q 2 (solid orange line) and the ones of ref. [18] at high-q 2 (dashed purple line). Note that we have plotted the results of refs. [17,18] beyond their region of applicability. In these plots we have set z = (0.29) 2 and = 10 −8 .

JHEP04(2020)012
In addition, it takes into account the sign of Im(s) correctly, as the functions F j i(a) have a different form in the upper or lower complex-s plane due to the double fixing of boundary conditions (i.e. section 3.5). The prescription for z is fixed as described above.
We have tested the results against those in refs. [17,18], finding very good numerical agreement with tables 1 and 2 in both papers. As already mentioned, the results of refs. [17,18] apply specifically to the low-q 2 and high-q 2 regions respectively. In figure 3 we have plotted these results within and beyond their respective regions of applicability and compared them with the analytic results obtained in this paper. We find an excellent agreement within the appropriate regions. Deviations with respect to the low-q 2 results occur starting around s −0.4. Thus, for the calculation of the OPE matching coefficients in this region it may be advisable to use the results given in the present paper.

Selected results at different values of s and z
The results for the NLO functions F (7,9) 1,2 (q 2 ) are intended to be used to calculate the function H µ in the OPE region, by means of eqs. (2.8), (3.37) and (3.38). For the determination of exclusive b → s amplitudes at large hadronic recoil, this OPE region corresponds to the region of negative q 2 [7,22]. As the m c dependence for these values of s is mild, a quadratic interpolation of the values at these three points will represent this dependence accurately enough.
5 Study of the analytic structure at NLO

Singularities of the NLO functions
The matching coefficients ∆C 7,9 (q 2 ) will mimic the analytic structure of the non-local form factors H λ (q 2 ) discussed in section 2.5. In this case the analytic singularities are due to on-shell intermediate partonic states in the b → s amplitude, producing branch cut discontinuities in both variables q 2 and (q + k) 2 . This structure can be observed explicitly in the analytic results for ∆C 7,9 (q 2 ) calculated here, where the contribution from each diagram to each singularity can be checked.
The expected singularity structure is the following. First, the analytic structure of each of the diagrams as a function of complex s ≡ q 2 /m 2 b can be chosen to have a branch cut on the positive real line above some specified (perturbative) threshold: s > s th , where the threshold depends on the diagram. In addition, some diagrams are real on the real line below the threshold, while some are complex-valued. This is due to the fact that some of the diagrams (the ones that are complex) contain on-shell cuts in the variable p 2 b ≡ (q +k) 2 , which we fix to p 2 b = m 2 b from the start. According to their (expected) analytic structure, the set of diagrams can be classified in four groups: 1. Diagram b 2 : Branch cut for s > 4, real for s < 4. The rest of the diagrams, a 1,3 and b 1,3 do not have branch cuts in the variable s because the photon couples to the external legs of the diagram. Note also that the specific threshold (4m 2 s /m 2 b , 4m 2 c /m 2 b or 4m 2 b /m 2 b ) can be determined from the charge coupling (whether the diagram is proportional to Q s , Q c or Q b ). This relates to the discussion in section 2.5, and applies also to the counterterm contributions.
From the explicit results obtained here for the contribution to ∆C 7,9 from each group of diagrams and counterterms, we can check this analytic structure. This is done in two steps: 1. Checking explicitly that the discontinuity lies where it is expected, and that the values of each contribution below threshold is real or complex as predicted.
2. Checking appropriate dispersion relations, thus supporting the absence of further singularities besides the expected branch cuts. This is done by checking, for each diagram class, the following equation:  Table 1. Values for the functions F (7,9) 1,2 (q 2 ) at negative q 2 , for three choices of z = m 2 c /m 2 b . The renormalization scale has been fixed to µ = m b . These numbers do not depend on whether one includes an infinitesimal positive or negative imaginary part for s.

JHEP04(2020)012
for any two points {s 0 , s 1 } in the complex plane. Any additional singularities will (generically) produce extra contributions beyond the integral in the r.h.s., and thus the fact that this dispersion relation holds is consistent with the absence of additional singularities anywhere on the complex plane, away from the real interval [s th , ∞).
Concerning the discontinuities along the real axis, figure 4 and figure 5 show the contribution to the form factors for each diagram class, evaluated above and below the real axis, for a reference value of z = 0.1. We see that the results obey the branch cut structure laid out above. Since the contributions from diagrams b, d and e are real below threshold, the branch-cut discontinuity is purely imaginary, as can be seen from the plots. On the contrary, the contributions from diagrams a and c are complex-valued below the thresholds since they have on-shell cuts in the variable p 2 b . This leads to a complexvalued branch-cut discontinuity (with a non-zero real part) in the ranges 0 < s < 4z and 4z < s < 1 respectively.
Besides explicitly confirming the expected branch-cut structure of the two-loop contributions, we find two features that we consider noteworthy: • The discontinuities in diagrams a and c become purely imaginary for s > 4z and s > 1, respectively.
• The contribution from diagrams c features a pole on the real axis when approaching the point s = 1 from the negative imaginary plane. This pole is related to an anomalous threshold.
The same structure of branch cuts is found for the various counterterms: discontinuities starting at s > 0, s > 4z and s > 1 for F ct (7,9) i,Qs (s), F ct (7,9) i,Qc (s) and F ct (7,9) i,Q b (s) respectively. We refrain from showing the corresponding plots for brevity.
Concerning the dispersion relation, we have checked that eq. (5.1) is satisfied with good numerical accuracy separately for all diagram classes, each with its corresponding threshold. To give an example, we consider F 2,(b) (s) with z = 0.1. As discussed above, this function contains a branch cut starting at s th = 4. We find that its discontinuity can be fitted approximately by Disc F Using this fit (for the sake of rapid integration) we find, for example taking s 1 = −3 + i and s 0 = −1 − 2i in eq. (5.1): 2,(b) (t) (t + 3 − i)(t + 1 + 2i) = 0.0894966 − 0.160839 i .

JHEP04(2020)012
As another example including a point at s 0 > 0: For s 1 = −1 and s 0 = 0.7, we find: again showing that the dispersion relation is very well verified. For applications with s 0 on the cut, the dispersion integral must include the prescription (t−s 0 −i ) in the denominator of the integrand, in order to regulate the pole (cf. eq. (2.23)). Thus, numerically the value taken for will determine the precision with which the discontinuity and the dispersion integral are evaluated.

OPE coefficients with flavor separation
At this point we can collect the separate contributions to the OPE coefficients ∆C 7,9 (q 2 ) proportional to the charge factors Q c and Q s/b . Denoting these two contributions by ∆C  (7) i,Qc , i(e) + F ct (9) i,Qc , i(a) + F i(b) + F ct (7) i,Qs + F ct (7) i,Q b , (5.9) i(a) + F (9) i(b) + F ct (9) i,Qs + F ct (9) i,Q b , (5.10) where in (5.7) we have omitted the term F i,(e) = 0. These OPE coefficients will contribute separately to the functions H OPE λ,c (q 2 ) and H OPE λ,sb (q 2 ) appearing in the two different dispersion relations in eq. (2.25). As discussed above, they have the proper analytic structure with branch cut discontinuities starting at s > 0 and s > 4z, for ∆C The corresponding results for F (j) 1,x are qualitatively similar. The conclusion is that, within the LCOPE region q 2 < 0, the contribution proportional to the charge factor Q c is in most cases a few times larger than the one proportional to Q s/b . , in the q 2 < 0 region. In these plots we have set z = (0.29) 2 .

Conclusions and outlook
The determination of non-local effects in exclusive b → s processes is of great phenomenological interest, but very challenging theoretically. These effects are associated with the matrix element of a bi-local operator (cf. eq. (2.5)), which is significantly more complex than the usual "local" form factors that govern the naively-factorizable part of the amplitudes (such as the ones arising from semileptonic and electromagnetic dipole operators). The current approach to non-local effects is to write an OPE for the bi-local operator in a kinematic region where the OPE converges (even if unphysical) and then to extrapolate the results to the physical region using analyticity or dispersion relations. At the level of the OPE, the non-local matrix element can then be expressed in terms of simpler form factors, and OPE coefficients that are determined from a perturbative matching calculation. The leading OPE coefficients have been known up to NLO for some time, but only in certain expansions on q 2 and/or z = m 2 c /m 2 b [17,18]. Here we have presented a recalculation of these two-loop contributions, fully analytic in both variables. This calculation has made use of the formalism of differential equations in canonical form, and the results are expressed in terms of Generalized Polylogarithms up to weight four. A particular attention has been put in obtaining an analytic continuation of the Feynman integrals with the desired singularity structure; for this purpose, special care is needed in fixing the integration constants in the solution of the differential equations. Numerically, our results agree with previously known expanded results within their range of applicability, but deviate notably for q 2 −10 GeV 2 .

JHEP04(2020)012
With the fully analytic results at hand, we have been able study the analytic properties of the non-local form factors, and we have confirmed the expectations from unitarity. In particular, we have verified the dispersion relations and checked the absence of singularities beyond the branch cuts from intermediate states in the q 2 channel.
In addition, we have presented the complete set of results separated into contributions proportional to different charge factors. This allows to study the extrapolation to the physical region separately for cc states, ss and bb states, and light states [7,22].
While the contributions from the operators O 1,2 considered here are the dominant ones in the SM for b → s transitions, it would be interesting to complete this calculation including the full set of four-quark operators in the general Weak Effective Theory [24]. This is important for an improved analysis beyond the SM [38], and also for the case of b → d transitions, where the up-quark contributions are not CKM suppressed [39]. The results for the renormalized two-loop functions F (7,9) 1,2 , as well as the separate contributions from each diagram class F We note that F170 = F270 = 0.
These routines operate by first collecting a list of the different GPLs that appear, in order to evaluate each GPL only once. This leads to a huge increase in the speed of the evaluation.

B List of Master Integrals
In this appendix we collect the list of all Master Integrals (MIs) J i,k that appear in the calculation of the two-loop diagrams a-e in figure 2. The notation is described in section 3.2. For diagrams a there are 7 MIs:

D Explicit examples for fixing integration constants
We first consider the master integral from diagram e with four propagators, i.e. J e,5 . Solving the corresponding differential equations in the canonical basis and then transforming the solution to the ordinary basis we get, for the

JHEP04(2020)012
where c 2 is an integration constant. Imposing the condition that J (−2) e,5 is nonsingular for s → 0 (which is equivalent to y e → 0), we get c 2 = − 1 2304π 4 , leading to (D.5) J c,2 has three propagators. J c,7 also has three propagators but one of them is squared. This means that J c,2 ∼ z and J c,7 ∼ z 0 for large z. Or in terms of x c , J c,2 ∼ x −2 c and J c,7 ∼ x 0 c when x c → 0. Imposing these conditions, we find c 1 2 = − 1 4096π 4 , leading to Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.