Bd,s → γℓℓ¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \gamma \mathrm{\ell}\overline{\mathrm{\ell}} $$\end{document} decay with an energetic photon

We calculate the differential branching fraction, lepton forward-backward asymmetry and direct CP asymmetry for Bd,s → γℓℓ¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \gamma \mathrm{\ell}\overline{\mathrm{\ell}} $$\end{document} decays with an energetic photon. We employ factorization methods, which result in rigorous next-to-leading order predic- tions in the strong coupling at leading power in the large-energy/heavy-quark expansion, together with estimates of power corrections and a resonance parameterization of sub- leading power form factors in the region of small lepton invariant mass q2. The Bd,s → γℓℓ¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \gamma \mathrm{\ell}\overline{\mathrm{\ell}} $$\end{document} decay shares features of the charged-current decay Bu → γℓν¯ℓ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \gamma \mathrm{\ell}{\overline{\nu}}_{\mathrm{\ell}} $$\end{document}., and the FCNC decays B → K∗ℓℓ¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {K}^{\left(\ast \right)}\mathrm{\ell}\overline{\mathrm{\ell}} $$\end{document}. As in the former, the leading-power decay rates can be expressed in terms of the B-meson light-cone distribution amplitude and short-distance factors. However, similar to B →K∗ℓℓ¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {K}^{\left(\ast \right)}\mathrm{\ell}\overline{\mathrm{\ell}} $$\end{document}, four-quark and dipole operators contribute to the Bd,s → γℓℓ¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \gamma \mathrm{\ell}\overline{\mathrm{\ell}} $$\end{document} decay in an essential way, limiting the calculation to q2 ≲ 6 GeV2 below the charmonium resonances in the lepton invariant mass spectrum. A detailed analysis of the main observables and theoretical uncertainties is presented.


Introduction
The semileptonic flavour-changing neutral-current (FCNC) decays driven by the quarklevel transitions b → q ¯ with q = d, s and = e, µ, τ are important tests of the dynamics of the quark and lepton sectors in the framework of the standard model (SM) and constitute very sensitive probes of nonstandard effects.
The exclusive B → K ( * ) ¯ decays are experimentally most easily accessible due to their comparatively sizeable branching fractions of O(10 −6 ). The angular distributions of their three-and four-body final states are among the prime targets of LHCb and Belle II for = e, µ, because they offer many CP-symmetric and CP-asymmetric observables. However, JHEP12(2020)148 the hadronic nonperturbative effects entering the theoretical treatment are complex ranging from form factors to resonant contributions, preventing currently theoretical calculations with percent accuracy. In contrast, the purely leptonic decays B q → ¯ allow theoretical control of hadronic effects at the percent level [1], because they depend -apart from higher-order QED corrections [2,3] -only on the B-meson decay constant, which can be computed in QCD lattice calculations with sub-percent accuracy [4]. On the other hand the helicity suppression of the two-body final state leads to tiny branching fractions of O(10 −9 ) for = µ, requiring huge data samples at LHC in order to perform measurements with precision comparable to the prediction. Nevertheless this decay has been observed by LHCb [5], CMS [6] and ATLAS [7] with a rate compatible with SM predictions.
Compared to the strong research activities on the aforementioned decay modes, the radiative leptonic decays B q → γ ¯ have received relatively little attention. The additional photon in the final state implies suppression by the electromagnetic coupling, but lifts the very different helicity suppression for = e and = µ present in B q → ¯ . The decay rates of B q → γ ¯ become comparable for the electron and muon final states and can be almost O(10 −8 ) for q = s [8][9][10][11]. Further, the angular distribution offers complementary observables to test the short-distance couplings. So far there are no experimental studies of the B q → γ ¯ decays, which for small photon energies also constitute a background process in the experimental analyses of B q → ¯ in the dilepton-invariant mass side-band [5,12,13].
Besides testing the SM, the theoretical interest in B q → γ ¯ also derives from the structure of hadronic effects. The limit of large photon energy E γ Λ QCD , where Λ QCD denotes the strong interaction scale, allows for a systematic treatment of nonperturbative effects when adopting the framework of QCD factorization and soft-collinear effective theory (SCET). As pointed out in [14], in this limit the leading nonperturbative corrections become universal for the decays B u → γ ν , B q → γγ and B q → γ ¯ , relating thus hadronic effects in b → u ν , b → qγ and b → q ¯ , respectively, and offering a programme for the combined analysis of long-and short-distance quantities in these decays. Other approaches are available in the literature relying on the B q → γ form factor calculation with dispersive methods and quark models [8,9,11,15]. The case of very soft photons due to initial-state radiation has been considered in the framework of heavy-hadron chiral perturbation theory in [12].
The factorization programme has been carried out for B u → γ ν starting within QCD factorization (QCDF) [16,17] and has then been extended to the two-step matching in SCET [18,19] at leading power (LP) in Λ QCD /m b , including next-to-leading order (NLO) radiative QCD corrections and resummation of the associated large logarithms. Next-toleading power (NLP) corrections at lowest order in QCD have been considered in [20], studying the consequences for the form factor symmetry relation. The NLP corrections have been scrutinized further with sum rule calculations [21][22][23][24].
The decay B q → γ ¯ appears as a hybrid of the charged-current decay B u → γ ν , which is driven by tree-level electroweak b → u ν transitions and the FCNC decays B → K ( * ) ¯ , which receive contributions from semileptonic b → s + (γ, ¯ ) as well as from hadronic b → s + (g, qq) transitions. When the final-state photon carries large energy relative to the strong interaction scale, the non-hadronic final state of the B q → γ ¯ decay JHEP12(2020)148 enables the calculation of the relevant form factors in terms of the B-meson light-cone distribution amplitude (LCDA) at LP in an expansion in the energy of the photon, similar to the decay B u → γ ν [18,19]. On the other hand, the contribution of b → sqq fourquark operators to B q → γ ¯ introduces issues familiar from B → K ( * ) ¯ decay, such as the presence of charmonium resonances in the ¯ invariant mass spectrum that limit the applicability of factorization methods to this decay [25,26].
Previous work on B q → γ ¯ [8][9][10][11]27] has usually focused on providing theoretical calculations of the (differential) decay rate in the entire kinematically allowed regions in terms of B → γ * form factors. This approach inevitably requires some amount of theoretical modelling of the form factors, and neglects "non-factorizable" O(α s ) corrections from the b → sqq four-quark and chromomagnetic dipole operators, which are known to be sizeable for B → K ( * ) ¯ decays [25].
In contrast, our focus is on the kinematically more restricted region of photon energy E γ Λ QCD , while making maximal use of factorization methods applicable in this limit. We perform a SCET analysis of B q → γ ¯ including O(α s ) QCD corrections at LP. We further include local NLP corrections in O(α 0 s ) explicitly, whereas non-local NLP contributions are parametrized following [20,22]. The analysis is more involved for the FCNC b → q ¯ transitions compared to the tree-level decay b → u ν due to the extended weak operator basis (2.1). Parts of our analysis are related to earlier works on B u → γ ν , B q → γγ and B → K ( * ) ¯ within QCDF [14,25,28] or SCET [18][19][20]22]. The factorization approach then allows us to express the decay amplitudes in terms of only the B-meson light-cone distribution amplitude (LCDA) at LP in the heavy-quark/large-energy expansion. The modelling of form factors is necessary only at sub-leading power in this expansion. The result is expected to be valid for sufficiently broad bins in invariant mass q 2 of the lepton pair below the charmonium resonances, q 2 6 GeV 2 , and above the narrow light-meson resonances. The latter restriction arises, since estimates following [29] imply that global duality is violated not only by charmonium, but also by the light-meson resonances. In order to extend our calculations to observables local in q 2 , we incorporate the light-meson resonances in our model for the sub-leading power form factors, which confirms the estimates of duality violation. This together with numerical cancellations of various otherwise dominant leading-power effects leads us to conclude that theoretical predictions of B q → γ ¯ observables in the low-q 2 region are plagued by large theoretical uncertainties.
The outline of the paper is as follows. In section 2 we detail the conventions for the effective weak interaction Lagrangian and parametrize the B q → γ ¯ amplitude in terms of hadronic tensors and the associated form factors. We employ a mix of QCD factorization and SCET techniques to calculate in section 3 the hadronic tensors at LP including O(α s ) radiative corrections and the summation of logarithms, and at NLP in O(α 0 s ). We next discuss issues related to the interpretation of the rescattering phase from factorization, duality violation and resonances. We then parametrize the single NLP form factor, which cannot be computed in factorization, and discuss the model that will be used for the numerical analysis. Section 4 is devoted to a detailed discussion of the various contributions to the B q → γ ¯ decay amplitude, which exhibits sizeable corrections and cancellations in the q 2 region of interest. We show the differential branching fraction, lepton JHEP12(2020)148 forward-backward asymmetry and direct CP asymmetry as functions of q 2 and integrated in various bins. We conclude in section 5. Further details on definitions, conventions, final-state radiation and the B-meson LCDA are given in appendices.

∆B = 1 effective theory
The low-energy effective theory of electroweak interactions in the SM for ∆B = 1 semileptonic b → q ¯ (q = d, s) decays is We adopt the basis proposed in [30,31] for the operators P i . The normalization factor t contains products of elements of the quark mixing matrix λ (q) for U = u, c, t. The term proportional to λ (q) u leads to tiny doubly Cabibbo-suppressed CPasymmetries in b → s transitions, but is not negligible in b → d transitions. The Wilson coefficients C i (ν) are evaluated at the hard scale ν of the order of the b-quark mass m b , after renormalization group (RG) evolution [32,33] from the electroweak scale. There are four-quark operators, the so-called charged-current operators P c,u 1,2 and the QCD-penguin operators P 3,4,5,6 . The remaining operators in the SM are the semileptonic operators with P L/R = (1 ∓ γ 5 )/2 and the dipole operators where m b is the running b-quark mass in the MS scheme at the scale ν. The QED × QCD covariant derivative is chosen as [25], 1 and Q = −1 for leptons, Q u = +2/3, Q d = −1/3 for quarks.

B q → γ ¯ form factors
The transition amplitude for the decay of the B q meson is where is the polarization vector of the photon, and k its momentum, related to the dilepton momentum q ≡ p + p = p − k. When working to lowest non-vanishing order in JHEP12(2020)148 the electromagnetic coupling α em = e 2 /(4π), but to all orders in QCD, the amplitude can be written in the form where j µ = Q γ µ denotes the leptonic electromagnetic current. The quantities T µν i and S (i) ν define hadronic matrix elements, which encode all QCD effects at lowest non-vanishing order in the electromagnetic interaction. The calculation of these hadronic tensors is one of the main purposes of this work.
We briefly consider the terms involving S (i) ν . The structure of the purely leptonic tensor multiplying these terms shows that they correspond to final-state radiation (FSR) of the photon from the lepton-pair. The hadronic matrix elements S (i) ν are B-meson to vacuum matrix elements, which must be proportional to the B-meson momentum p ν . The contraction of the vectorial leptonic tensor with p ν vanishes, p ν L µν V = 0, while the corresponding contraction of the axi-vectorial leptonic tensor is proportional to the lepton mass. The only non-vanishing FSR contribution is therefore due to P 10 [34] and helicitysuppressed. The QCD effects are entirely described by the B-meson decay constant, since The FSR contribution is tiny except in very small phase-space regions, and can be safely neglected in the later numerical analysis. For further details on the FSR contribution we refer to appendix B. Our main concern are therefore the hadronic tensors T µν i , sometimes referred to as structure-dependent (SD) contributions. The individual SD contributions of each operator are contracted with the vector or axial-vector lepton currents . (2.7) The hadronic tensors T µν i (k, q) are the correlation functions 6, 8 (2.11) JHEP12(2020)148 Figure 1. The three possible contractions (A-type left, B-type middle) of the four-quark operators P 1,..., 6 at O(α 0 s ). The insertion shown on the right is power-suppressed. The B-type insertion is power-suppressed when For the first two diagrams, there exists a counterpart where the real/virtual photon is attached to the heavy-quark line (double line), which, however, is power-suppressed. The external photon line with momentum k symbolizes the electromagnetic current j µ f (x) in the definition (2.11) of the hadronic tensor T µν i , the virtual photon including its decay into the lepton-pair graphically represents the current j µ f (y). The B smeson state is represented by its leading partonic constituents, a b-quark and a light anti-quark.
with the electromagnetic quark-current Here f is summed over all five active quark flavours. The hadronic tensors T µν 9,10 are time-ordered products of the weak currents with the electromagnetic current in direct analogy with the charged-current decay B u → γ ν . The contribution from P 7 , T µν 7 = T µν 7A + T µν 7B , has been split into the part T µν 7A , which accounts for the emission of the final-state photon from the constituents of the B-meson, and the part T µν 7B , where the on-shell photon is emitted directly from P 7 , while the electromagnetic current produces the virtual photon with momentum q, which decays into the lepton pair. Quite generally, we refer to these two different contributions as A-and B-type, respectively. 2 Similarly, the hadronic tensors T µν 1,...6 for the four-quark operators contain A-type and B-type parts, depending on which of the two electromagnetic currents is inserted into a B-meson constituent quark line, and further an annihilation-type contribution, where neither of the on-shell and virtual photon is emitted from the initialstate quarks. See figure 1 for the lowest order diagrams in α s representing these three contributions.
The hadronic tensors T µν i (k, q) can be parametrized in terms of scalar form factors. Their number is determined by exploiting the Ward identities (e.g. [35]) that hold upon contraction with k and q, together with the algebraic relations k µ T µν 7B = q ν T µν 7A = 0. The Ward identities and the transversality of the on-shell photon, · k = 0, which implies that µ is a transverse index, JHEP12(2020)148 lead to the Lorentz decomposition which contains two helicity form factors F (i) L,R for each operator P i . We use the convention ε 0123 = −1, 3 for definitions of g µν ⊥ and ε µν ⊥ we refer to appendix A. Below we will employ QCD factorization techniques to factorize the hadronic matrix elements T µν i (k, q) and the corresponding form factors F In this limit the right-helicity form factors F L due to the left-handedness of the weak interaction and helicityconservation of QCD at high energy. One often defines vector F A between the transversity form factors holds at LP in the heavy-quark/large-energy expansion.

Factorization of form factors
QCD factorization can be applied to the hadronic process if the on-shell photon is very energetic E γ Λ QCD . It is most intuitive to work in the rest frame of the B q meson, where the three-vectors k and q are back-to-back and the constituent b and q quarks of the B q meson will have soft residual momenta of order Λ QCD .
We introduce the four velocity v (v 2 = 1) of the B q -meson and a pair of light-like vectors n − and n + , with n 2 − = n 2 + = 0, n − n + = 2, v = (n − + n + )/2, such that The photon energy E γ = (n + k)/2 = (m 2 Bq − q 2 )/(2m Bq ). For later convenience we define y ≡ 2E γ /m Bq with 0 ≤ y ≤ 1 − 4m 2 /m 2 Bq . We refer to momenta r with large component n + r = O(E γ , m b ) in the direction of n − and small r 2 as collinear. Anti-collinear momenta have this property with n + ↔ n − interchanged. Hence the momentum k µ of the on-shell final-state photon is collinear. The nature of the lepton-pair or virtual photon momentum q µ depends on whether the real photon energy E γ is close to its maximal kinematically allowed value m Bq /2, corresponding to y = 1. Since q 2 /m 2 Bq = 1 − 2E γ /m Bq = 1 − y, we consider two scalings for q 2 : 1) anti-hard-collinear (hc) q 2 ∼ E γ Λ QCD or m b Λ QCD , in which case the momentum q µ is dominated by its component proportional to the light-like vector n µ + , and 2) hard (h) q 2 ∼ m 2 b , in which case the photon energy is large, but y is not parametrically close to 1 in the heavy-quark limit.

JHEP12(2020)148
The LP A-type and B-type insertions. The grey blob denotes the effective hard vertex of the operators P 1,...,8 .
In practice, we will be interested in the region with the restriction to q 2 6 GeV 2 , as originally introduced for B → K ( * ) ¯ [25], to avoid the charmonium resonance region. This upper limit lies somewhat in between these two scalings. We will construct the form factors (2.13) in an expansion in powers of where the LP contribution includes the resummation of NLO radiative QCD corrections, for which we employ SCET. The NLP contribution is included to LO in QCD. Our result will be accurate to these orders for both possible scalings of q 2 .

Form factors at LP
The SCET framework allows for a systematic decoupling of fluctuations with hard virtualities of order m 2 b in the matching QCD → SCET I and fluctuations with hard-collinear virtualities of order E γ Λ QCD , m b Λ QCD in the matching SCET I → SCET II . The SCET framework at LP [18][19][20] can be directly applied to the semileptonic operators P 9,10 and there is no difference whether q 2 ∼ hc or q 2 ∼ h. The factorization of the form factors of the other operators i = 1, . . . , 8, which admit A-and B-type contributions is more complicated, and we provide the results below, together with some explanation of their derivation.

Effective hard vertex
We begin the discussion of the first matching step to SCET I by assuming that q 2 is hard and the insertion is of the A-type. A LP contribution to the hadronic tensors T µν i is obtained only when the separation x between the electromagnetic current j µ f (x) and the operator P i (0) is of order 1/ E γ Λ QCD . In other words, the propagator joining the two vertices in the left diagram of figure 2 must have hard-collinear virtuality, since a hard propagator would lead to power-suppression. Integrating out the hard scale therefore results in an effective flavour-changing vertex, represented by the grey circle in figure 2. The hadronic tensors appearing in (2.5) are therefore matched to SCET I as follows:

JHEP12(2020)148
where H i (q 2 ) denotes the hard matching coefficients, and the hadronic tensor for i = 7 is T µν 7A . All operators match to the same correlation function of the SCET I heavy-to-light current q hc γ ν ⊥ P L h v and the representation of the SCET I electromagnetic current. In the above equation the Wilson coefficients C i are evaluated at the scale ν ∼ m b after evolving them to ν from the electroweak scale. The matching coefficients H i are determined at the hard scale µ h ∼ m b in the matching QCD → SCET I . In addition to q 2 , they depend on the scale ν through the ultraviolet divergences of the QCD diagrams, such that the dependence on ν cancels in the product C i (ν)H i (ν, µ h ), as well as on the scale µ h through the infrared divergences of the QCD diagrams. The µ h -dependence cancels with the dependence on µ h of the SCET I correlation functions in The O(α s ) corrections to the hard functions H i (q 2 ) are identical to those for B → K ( * ) ¯ (see figure 2 in [25]). For the four-quark operators, they arise from two-loop diagrams. It is customary to absorb the four-quark operator contributions into effective Wilson coefficients of P 7−9 , and we follow this practice here. Hence we define where m b denotes the MS-scheme b-quark mass at the scale ν. 4 The effective hard functions including the O(α s ) correction are then given by 8 + C 3−6 -terms , where we have suppressed all arguments q 2 , µ h , ν. Here α s is evaluated at µ h . The effective Wilson coefficient combinations read as follows [30,31]:

JHEP12(2020)148
We use the definition of the function Y (q 2 ) from [25]. The function h(q 2 , m q ) [31] depends on the light quark masses m u,d , which are set to zero, or the charm-quark pole mass m c . For the semileptonic operators and the parts from the four-quark operators contained in the effective Wilson coefficients above, the NLO QCD corrections are contained in the matching coefficients C T 1 (q 2 ) of the heavy-to-light (axial-) vector and tensor QCD currents to the SCET I current q hc γ ν ⊥ P L h v [36]. We use the notation of [37] and provide the NLO expression here: , and N c = 3 in QCD. The NLO QCD corrections from the charged-current operators P u,c 1,2 and P 8 are given by F (7u,7c,9u,9c) 1,2 (q 2 ) from [38][39][40] and F 7,9 8 (q 2 ) from [25,39], respectively. The analogous corrections from P 3,4,5,6 labelled "C 3−6 -terms" in (3.6), (3.7) are suppressed by the small penguinoperator Wilson coefficients and have been neglected here. 5 The terms proportional to λ (q) u in V eff 7,9 give rise to direct CP violation, which is doubly-Cabibbo suppressed for q = s and leads to tiny CP asymmetries, whereas for q = d CP-violating effects can be larger.
To summarize the discussion up to this point, we note that the A-type contribution to B q → γ ¯ amplitude (2.5) after matching to SCET I at LP is obtained in the form (3.14) The derivation of this result assumed that the lepton-pair virtuality is hard. When q 2 is anti-hard-collinear, the structure of the result remains the same. However, the effective vertices, which are functions of q 2 /m 2 b , could now be evaluated at q 2 = 0, since the ratio q 2 /m 2 b is power-suppressed. More importantly, the effective vertex V eff 7 (q 2 ) becomes powerenhanced due to the photon pole 1/q 2 , relative to which V eff 9,10 (q 2 ) should be dropped as power-suppressed. This, however, results in a physically unacceptable approximation. For example, making the same approximation for the exclusive B → K ( * ) ¯ decay, the lepton forward-backward asymmetry would disappear and the predicted branching fractions and angular distributions in the small-q 2 region become unrealistic. The reason is that the magnitude of the Wilson coefficients C 9,10 is about an order of magnitude larger than C 7 , so that it is more appropriate to count C 9,10 ∼ m 2 b /q 2 × C 7 . In other words, the expansion in 1/m b should be performed for the C 9,10 and C 7 terms separately, and for each term the JHEP12(2020)148 LP should be kept. 6 In the same spirit, no harm is done by keeping the q 2 -dependence in the effective vertices also at small q 2 . In this way, one obtains a smooth interpolation between hard and anti-hard-collinear q 2 , as (3.14) applies without modification to the antihard-collinear, small-q 2 region. This is also the procedure that has been adopted for the exclusive B → K ( * ) ¯ [25], which forms the basis of the QCD phenomenology of this decay.
Another subtlety needs to be mentioned here. When q 2 is anti-hard-collinear, the effective hard vertices are no longer guaranteed to be dominated by hard virtualities. This is obvious for the one-loop four-quark operator contribution shown in the left diagram of figure 1, where the quark-loop is then purely anti-hard-collinear. At O(α s ) the two-loop diagrams contributing to F (7u,7c,9u,9c) 1,2 contain hard and hard-collinear regions. It therefore appears that the previous treatment should be substantially modified, since only the hard regions are integrated out in the first matching step to SCET I . However, the SCET I correlation function in (3.14) is governed by (hard-)collinear physics. The anti-(hard-)collinear and the (hard-)collinear sectors of SCET I are already decoupled after integrating out the hard modes and performing the soft-decoupling transformation in SCET I . The soft Wilson lines from the decoupling transformation cancel, as the anti-collinear final state is colourneutral, hence the existence of the anti-hard-collinear physics leaves no impact on the remaining collinear and soft interactions. In consequence, when q 2 is anti-hard-collinear, we may simply assume that the hard functions H i (q 2 ) introduced above are in fact the hard and anti-hard-collinear functions. In principle, these functions may therefore contain formally large logarithms ln q 2 /m 2 b , which we cannot resum by not factorizing the hard and anti-hard-collinear physics properly. However, the existence of the limit q 2 → 0 shows that such logarithms are absent, at least at LP and O(α s ). We note that the same discussion applies to the standard treatment of the exclusive B → K ( * ) ¯ mode, although we are not aware of its explicit mentioning.
In addition to the A-type insertion, there is also a B-type contribution, see the right diagram of figure 2, in which the role of the real and virtual photon is interchanged. When q 2 is anti-hard-collinear, the separation between the flavour-changing and electromagnetic current is again of order 1/(E γ Λ QCD ), and one obtains a LP contribution. Since the lepton-pair originates from the electromagnetic current, there is no B-type contraction from the operator P 9,10 , and in fact no contribution from the corresponding effective vertices V eff 9,10 . Hence We note that the effective vertex V eff 7 (0) is now evaluated for k 2 = 0, 7 and the Fourier transform of the SCET I correlation function is taken with respect to q rather than k. The remarks above concerning the interpretation of the "hard" functions H i (q 2 ) when the argument q 2 is anti-hard-collinear, apply to the present case k 2 = 0. 6 In case of C eff 8 , one then counts F (7,9) 8 (q 2 ) as C7 and C9 terms, respectively. 7 To be precise, the definition now contains the hadronic tensor T µν 7B instead of T µν 7A , but with the factor 1/q 2 taken out, the result for V eff

JHEP12(2020)148
When q 2 is hard, the SCET I correlation function in (3.15) is not the correct expression, since the distance x ∼ 1/ q 2 is hard and should already have been integrated out in the matching to SCET I . In other words, the propagator connecting the two currents in the right diagram of figure 2 is far off-shell and should be contracted to a point, resulting in a local operator rather than a SCET I correlation function. Also the vertex V eff 7 (0) is not necessarily the correct one in (3.15), since it was obtained by matching the flavour-changing current with an on-shell out-going quark. On the other hand, the off-shell propagator provides power suppression of the B-type contraction in the hard-q 2 region relative to the anti-hard-collinear one, such that B-type insertions are NLP effects for hard q 2 , for which we aim only at O(α 0 s ) accuracy. We now argue that within the approximation that we include only LO O(α 0 s ) contributions at NLP, the expression (3.15) can be used even when where l is the momentum of the spectator quark. The LP approximation in both q 2 regions is obtained from keeping only the potentially large momentum components in the numerator. For hard q 2 this includes (n + q) The important observation is that this term vanishes, since µ is transverse and k = E γ n − , hence / n − γ µ ⊥ , / n − = 0. Thus, the previous equation simplifies to in both q 2 -regions. Eq. (3.17) would be obtained from (3.15) in the anti-hard-collinear q 2 region, which proves that we can smoothly extrapolate (3.15) into the hard region at Since we do not aim at O(α s ) accuracy at NLP, we can use (3.15) in the hard-collinear and hard q 2 region, even if in the hard region the O(α s ) correction to V eff 7 (0) is not the complete one at NLO in α s .
To summarize the first matching step, we obtain the LP B q → γ ¯ amplitude (2.5) in the form where the two terms are given by (3.14), (3.15). The second step consists of matching the SCET I correlation function in these terms to SCET II by integrating out the hard-collinear modes. Before turning to this task in the following subsection, we mention that the same correlation function appears for all operators P i from the electroweak effective Lagrangian. The RG evolution from the hard to the hard-collinear scale is therefore universally related to JHEP12(2020)148 the anomalous dimension of the LP A0-type SCET I heavy-to-light current q hc γ ν ⊥ P L h v [36] and equals the one that appears in the B u → γ ν decay. We therefore have The RG equation for U H (q 2 , µ hc , µ h ) and its solution to NLL 8 can be found in eqs.
in terms of the evolution factor defined there.

Hard-collinear function
The matching of SCET I to SCET II amounts to integrating out the hard-collinear modes in the SCET I correlation function No new calculation is required in this step. For r 2 = 0, it has first been explained and computed to the one-loop order for B u → γ ν in [18,19]. The generalization to r 2 = 0 needed here was worked out in [22]. We recall from [18] that the SCET I electromagnetic current j µ q,SCET I (x) relevant here consists of two power-suppressed pieces . Here (and above) q hc = W † ξ hc is the hard-collinear quark field multiplied by the QCD hard-collinear Wilson line. The second term on the right-hand side describes the conversion of the soft spectator quark into a hard-collinear quark at the photon vertex and counts as O(λ 2 ). The first term is only O(λ) suppressed, but contributes to the amplitude with an external soft quark only through the O(λ) suppressed quark-gluon vertex of the SCET Lagrangian, that converts a soft quark into a hard-collinear quark through interaction with a hard-collinear gluon. This part of the current contributes to the matching only from NLO in the α s expansion through hard-collinear gluon corrections to the photon vertex.
With the hard-collinear physics at the scale E γ Λ QCD , m b Λ QCD integrated out, the remaining nonperturbative soft physics is parametrized at LP in the matching by the leading-twist light-cone distribution amplitude of the B meson, φ + (ω). The matching equation reads (3.23)

JHEP12(2020)148
Since both cases will be needed, we employ the convention that n · r denotes the large component n + r, when r is a collinear momentum and n − r when it is anti-collinear. 9 Concretely, for the A-type contribution, we need J(n + k, 0, ω) with n + k = 2E γ , while for the B-type one the relevant function is J(n − q, q 2 , ω) with n − q = m Bq and The scale-dependent quantities F Bq , J, φ + are assumed to be evaluated at the hardcollinear scale µ hc , and the scale argument has been omitted in (3.23). Their definitions are as follows. The B-meson LCDA φ + (ω) [42,43] is the Fourier transform of the HQET matrix element where [tn − , 0] denotes a straight soft Wilson line connecting the light-like separated points 0 and tn − . It is customary to relate the HQET B-meson decay constant F Bq to the scaleindependent decay constant of full QCD f Bq , introduced in (2.6), and use the latter as an input. The relation is The RG evolution factor U −1 F from µ hc to µ h can be obtained from the appendix of [20] Finally, the hard-collinear matching function ("jet function") reads [22] J(n · r, r 2 , ω; with n · r = r 2 /n · r. We note that the second line vanishes for r 2 = 0. Hence for the A-type insertions, the convolution J(ω) ⊗ φ + (ω) in (3.23) can be expressed as in terms of the inverse (λ Bq ) and the first two inverse-logarithmic moments (σ of the B-meson LCDA [20]. However, for the B-type insertions the expression is more complicated and the entire function φ + (ω) must be known to evaluate the convolution integral. JHEP12(2020)148

Final factorized form
Putting together (3.14), (3.15) and (3.23), we obtain the following compact result for the LP B q → γ ¯ amplitude (2.5): Note that the amplitude contains g µν ⊥ + iε µν ⊥ and not g µν and the so-called form-factor-symmetry relation F A as a consequence of helicity conservation in the heavy-quark and large-energy limit.
For completeness, we also give the left-handed form factors. For this purpose we define and find We recover the result for B u → γ ν [20] from (3.34), by setting C 9 = 1 and C i = 0 otherwise.

Form factors at NLP
As for the case of B u → γ ν [20] we include NLP Λ QCD /E γ , Λ QCD /m b corrections to the above B q → γ ¯ form factors, but we aim only at O(α 0 s ) accuracy at NLP. Such power corrections arise from three sources: • The coupling of the real or virtual photon to the heavy quark.
• Power corrections to the (anti-) hard-collinear light-quark propagator in the LP O(α 0 s ) contributions, that is, power corrections to the SCET I correlation function (3.21) at tree level.

JHEP12(2020)148
• Annihilation-type insertions of the four-quark operators, such that the real and the virtual photon are attached to the quark loop, see the right diagram in figure 1.
The first two effects are also present in B u → γ ν . The chromomagnetic dipole operator P 8 can be ignored, since its insertions involve at least one power of α s . The NLP contributions from the semileptonic operators P i (i = 9, 10) can be taken directly from the B u → γ ν calculation [20]: At NLP the right-helicity form factor is non-vanishing, but this "symmetry-breaking" contribution (as it implies By this we mean that at this order they can be expressed in terms of the B-meson decay constant. On the other hand, the power correction to the left-helicity form factor cannot be factorized and is parametrized by an unknown function, the "symmetry-preserving" soft form factor. We also refer to these contributions as "non-local" power corrections. For q = u relevant to B u → γ ν it has been calculated with QCD sum rules [21,22,24]. The definition in (3.35) is such that in the SU(3)-flavour symmetry limit of QCD [20] for B u → γ ν simply by the ratio of the electric charges of the spectator quarks. We remark that we consistently set the strange-quark mass to zero in our analysis, which would otherwise give m s /(2E γ ) corrections to the above expressions.
The NLP contributions of P 7 from the A-type insertion are We note that the emission of the photon from the heavy quark is symmetry-preserving for P 7 . For the B-type insertion of P 7 we find The parameterization of the NLP B-type insertion requires another soft form factor ξ Bq (E γ ). The SCET I correlation function (3.21) can be considered as a function of n · r and r 2 , where n · r is the large component of the (anti-)collinear momentum r. To parametrize the NLP correction, we may introduce the soft form factor ζ Bq (n · r, r 2 ) of two variables. Similar to the general hard-collinear function J(n · r, r 2 , ω) in (3.23), the soft-form factors in (3.37), (3.39) are then the special cases The effect of the four-quark operators at O(α 0 s ) is two-fold. First, the NLP terms from the left and middle diagram of figure 1 replace C 7 and C 9 multiplying F above by C eff 7 and C eff 9 . Second, they allow annihilation-type contractions (right diagram in figure 1), which appear only at NLP, and must be computed separately. The complete set of such diagrams is depicted in figure 3. Contributions of this type were computed for the case of two real photons, B q → γγ in [28]. The authors of this paper also showed that no infrared singularities appear in the two-loop O(α s ) corrections to these diagrams, hence allowing for a quantitative interpretation of the lowest-order one-loop diagrams.
We computed the one-loop annihilation contribution for the B q → γ ¯ situation of one virtual and one real photon attached to the quark loop. As expected, the result is local such that the nonperturbative hadronic physics can be expressed in terms of f Bq . We express the result in the form where the minus (plus) sign refers to L (R). The functions f V,A (y) read Here y = 2E γ /m Bq as before. We further introduced the mass ratio z Q = m 2 Q /m 2 Bq , where Q = u, d, s, c, b denotes the various quark flavours. In the numerical evaluation we set

Discontinuity, duality and validity of the factorization approach
To proceed to the numerical predictions for the B q → γ ¯ decay rates, we need a model for the generalized soft form factor ζ Bq (n · r, r 2 ), or the two single-variable form factors in (3.41). The form factor ξ Bq (E γ ) = ζ Bq (2E γ , 0) that appears in A-type contributions could be computed with QCD sum rules as for B u → γ ν , for which sophisticated results including radiative corrections and higher-twist effects already exist [21][22][23][24]. This method applies to the transition to a real photon or a photon with Euclidean virtuality k 2 < 0, but not to the case in the B-type contribution. This form factor develops an imaginary part and resonances, which cannot be fully described with large-energy factorization methods. In the following we provide some general considerations on the factorization calculation of B → γ * form factors in the physical region q 2 > 0. We take note of a recent computation of these form factors with QCD sum rules [44], which however applies to larger q 2 than of interest here.
We first note that the B q → γ ¯ decay amplitude contains a discontinuity already at LP from two sources. First, in the A-type contribution from V eff 7,9 (q 2 ), in lowest order given by the cut through the quark loop in the left-most diagram of figure 1. This discontinuity is similar to the one for B q → V ¯ and its implications for short-distance and factorization calculations are well understood. For the b → s transition it leads to prominent charmonium resonances in the q 2 -spectrum and large global parton-hadron duality violation for integrated or binned spectra, which limit the short-distance calculation to q 2 6 GeV 2 . On the other hand, the light-meson resonances related to light-quark loops cause negligible amounts of parton-hadron duality violation when sufficiently wide bins in q 2 are considered, as explained in [29]. The second source of a discontinuity appears in the B-type contribution from real intermediate states in the correlation function T µν (q) (3.21) in the physical region q 2 > 0 (see right diagram in figure 2). In the following we are concerned with this discontinuity.
From (3.30) or (3.33) we obtain where we used that the discontinuity is restricted to ω < q 2 /m Bq to set the upper integration limit. The second line holds in the tree-level approximation J(m Bq , q 2 , ω) = 1. Since φ + (ω) ∝ ω for ω → 0, the discontinuity survives as q 2 → 0, but is negligible relative to JHEP12(2020)148 the real part of F (7-eff,LP) L , which develops the photon pole. This partonic discontinuity must be interpreted as dual to a large number of continuum hadronic intermediate states in the sense of parton-hadron duality. When q 2 ∼ E γ Λ QCD is hard-collinear or larger, this is certainly the case, since the invariant mass of the hadronic intermediate state becomes parametrically large in the heavy-quark/large-energy limit. On the other hand, for small q 2 ∼ Λ 2 QCD , the partonic discontinuity becomes unreliable. Note however, that it also becomes power-suppressed, since ω * ∼ Λ 2 QCD /m b when q 2 ∼ Λ 2 QCD , hence φ + (ω * ) ∼ 1/m b as opposed to φ + (ω * ) ∼ 1/Λ QCD for hard-collinear q 2 . In reality, these parametric estimates do not work well. For example, when q 2 ∼ m 2 φ is near the φ meson resonance, m φ formally counts as Λ QCD , but ω * ∼ 200 MeV is closer to the QCD scale Λ QCD than to Λ 2 QCD /m b . We next compare the leading-power B-type contribution A type−B defined in (3.15) to the amplitude generated by saturating the hadronic tensor with a single vector resonance V of mass m V and width Γ V . Employing a standard spectral representation of T µν 7B (k, q) (2.10), we obtain (3.48) The vector meson decay constant and B → V transition form factors are defined by and respectively, with q = p − p . The factor a (q) V arises from the quark flavour wave function of the resonance. For the cases of interest below -the φ, ρ and ω mesons -the nonvanishing constants are a (s) The constant c V in (3.48) collects these flavour factors as well as the electric charges from the matrix element 0 j µ f V (q, ε) of the electromagnetic current, with values −1/3, −1/2, 1/6 for V = φ, ρ, ω. We further used T 2 (0) = T 1 (0) to obtain (3.48). Eq. (3.48) should be compared to the first and last line of (3.30). Hence, we obtain the resonant amplitude from the LP B-type amplitude by the substitution (3.52)

JHEP12(2020)148
With , we find that the resonant amplitude is suppressed by (Λ QCD /m b ) 2 in the heavy-quark limit when q 2 is hard-collinear and by Λ 2 QCD /(m b Γ V ) in the resonance region. This suggests that we can add the resonant amplitude to the parameterization of the soft form factor ζ Bq (n · q, q 2 ) without double counting of part of the short-distance contributions, since they are formally of lower order in the heavy-quark expansion. 10 In order to address the question whether global duality is violated by the presence of resonances, we consider the ratio of differential decay rates in a q 2 -bin. Since we are only interested in an O(1) estimate, we evaluate the resonance contribution in the narrow-width approximation, and the shortdistance contribution in the tree-level approximation J = 1 for the hard-collinear function.
We define the complex, q 2 -dependent inverse moment of the B-meson LCDA, and obtain The second line is obtained upon neglecting the q 2 -dependence of λ Bq (q 2 ), and neglecting q 2 /m 2 Bq in the integrand. The above estimate allows us to draw a number of important conclusions: • Since the resonance is localized while the short-distance amplitude is smooth, we would have expected the ratio to decrease as 1/q 2 max with the width of the integration interval for q 2 max q 2 min . Instead the impact of the resonance on the integrated decay spectrum decreases only logarithmically. This behaviour appears because the amplitude shows the photon-pole 1/q 2 enhancement.
• In the large-energy/heavy-quark limit R ∼ (Λ QCD /m b ) 2 , counting Γ V ∼ Λ QCD ; as expected the resonance is a sub-leading power correction. However, R can be strongly enhanced for narrow resonances. In [29] this mechanism has been identified as the source of large global duality violation for the inclusive B → X s ¯ decay. Here we find (with parameters as specified in section 4 and omitting the bin-size dependent JHEP12(2020)148 logarithm in this estimate) R ≈ 57 for the φ resonance contribution to B s → γ ¯ , and R ≈ 2.8 and 3.9 for the ρ and ω meson contribution, respectively, to B d → γ ¯ . Hence we conclude that global duality is violated by a large factor in B s decay due to the narrow width of the φ meson, and there is still an O(1) contribution from the ρ and ω resonance for the B d case. The short-distance contribution is only dominant, if the q 2 bin of dΓ/dq 2 does not include the resonance and is sufficiently in the hard-collinear region.
• Eq. (3.55) is similar to eq. (44) of [29], except that the factor f 2 here. The first factor arises from the production of the resonance from a local current, while here the resonance arises from the transition from an extended B q meson. Although both ratios scale as (Λ QCD /m b ) 2 in the heavyquark limit, the second is larger by nearly two orders of magnitude, since the B qmeson form factors and decay constant do not satisfy the heavy-quark mass scaling at m b ≈ 5 GeV. This explains why the ρ resonance makes a negligible contribution for the inclusive B → X s ¯ transition discussed in [29], but is O(1) for B d → γ ¯ .
Since the prominent lowest-mass resonance(s) may dominate any bin of the q 2 -distribution, which contains them, in particular for B s → γ ¯ , the short-distance calculation is valid only in a small q 2 region from approximately 2 GeV 2 to about 6 GeV 2 below the charmonium resonances. To extend the theoretical prediction to smaller q 2 , we propose an ansatz for the soft NLP form factor that includes the lowest resonance(s). From the above discussion we deduce that this can be done without double counting, but the ansatz departs from the rigour of the short-distance calculation.

Ansatz for the soft NLP form factor
The form factors ξ Bq (E γ ), ξ Bq (E γ ) are suppressed by a single power of Λ QCD /E γ in the large-energy/heavy-quark limit. As mentioned above we add to these form factors a resonance contribution, which is formally suppressed by another power, to extend the local description of the q 2 spectrum into the low-q 2 region. We further draw on the observation from [24] that the power-suppressed soft form-factor contribution tends to have opposite sign to the leading-power contribution. We therefore subtract from ξ Bq (E γ ), ξ Bq (E γ ) the leading-power contribution in the tree approximation multiplied by a parameter r LP , which formally scales as Λ QCD /E γ and for which we adopt the value r LP = 0.2 ± 0.2. This leads to the ansatz where λ Bq = λ Bq (q 2 = 0) is the standard q 2 -independent inverse moment of the B-meson LCDA. The first term in each expression involves the tensor form factor T Bq→V 1 (q 2 ) of B q → V transitions, for which we will employ the combination of LCSR and lattice results

JHEP12(2020)148
from [45] in the numerical evaluation below. For the form factor ξ Bq (E γ ), which appears in the A-type contribution, QCD sum rule calculations [21][22][23][24] exist including radiative corrections and higher-twist effects. However, for a coherent approximation of both types of form factors, we choose the same ansatz here. The two expressions above can be considered as special cases of the same resonance+factorization ansatz for the generalized soft form factor ζ Bq (n · r, r 2 ). We also remark that the tensor form factor literally appears only for the 7A contribution to F L . For i = 9, 10, the vector and axial-vector form factors would appear. However, all these form factors reduce to the universal, transverse vector meson form factor ξ Bq→V ⊥ (q 2 ) when radiative and power corrections are neglected [43,46]. Within this approximation, it is consistent to use T By default we only include the φ(1020) resonance in case of the B s → γ ¯ decay. The higher ss resonances have large masses and we assume that they can be subsumed in binned averages of the short-distance contribution. 11 For the B d → γ ¯ decay, we include V = ρ, ω.

Phenomenological analysis
The phenomenological analysis of B q → γ ¯ for q = s, d and = e, µ will use the results of the form factors presented in the previous section. The systematic factorization leads to a reduction of renormalization scheme dependencies when including NLO QCD corrections to the LP contribution. Further, the NLP contributions allow to study the breaking of the form factor symmetry and provide estimates of nonperturbative resonant and nonresonant contributions in the large-energy limit, i.e. the low-q 2 region. Many features of the factorization approach can be investigated at the level of amplitudes, for which we refer to section 4.1. The observables of phenomenological interest and our final results will be presented in section 4.2.
The numerical values of the SM and hadronic parameters, which will be used in the analysis, are collected in table 1. The strong coupling α s (µ) in the MS scheme is calculated from α s (m Z ) with n f = 5 using three-loop evolution, including quark flavour threshold crossings at the scale µ 4 = µ h (transition to n f = 4) and µ 3 = 1.2 GeV (n f = 3). The Wilson coefficients C i of the weak EFT are calculated at the electroweak matching scale µ W = 160 GeV and then evolved with the required accuracy (following [32,33] The bottom-and charm-quark masses, which enter the one-and two-loop matrix elements of the four-quark operators P 1,...,6 through the hard functions in (3.6), (3.7) and (3.11), are usually chosen to be the pole masses. The pole masses suffer from large JHEP12(2020)148    table 1 with two-loop expressions and an uncertainty spanned by using oneand three-loop ones. On the other hand, the MS masses are also not appropriate here, at least for the charm quark, since the hard functions then exhibit imaginary parts from unphysically small values of 4m 2 c . A good compromise is the potential-subtracted (PS) renormalization scheme [63], since the PS mass has a well-behaved relation to the MS mass while being numerically closer to the physical thresholds. We therefore use the PS masses as parameters and perform the scheme conversion of the hard functions from the pole to the PS scheme for the masses. In practice, this has to be done only for the charm mass, such that the NLO corrections F We shall see below that the parameters of the B-meson LCDA cause by far the largest theoretical uncertainty from input parameters. This is not surprising, since the related charged current process B u → γ ν is usually advocated as a measurement of λ Bu ≈ λ B d . While λ B d may therefore be known more precisely in the future, there is no obvious measurement of λ Bs . There also do not exist studies of SU(3) breaking for this quantity, which requires us to make an educated guess. The value assumed in table 1 is obtained from the fact that the B-meson LCDA represents the distribution of light-cone momentum of the light-quark in the meson. A quark with larger mass is expected to have a larger JHEP12(2020)148 light-cone momentum on average. More details on the treatment of the B-meson LCDA in this analysis are provided in appendix C.
If not stated otherwise, below we provide results for q = s and = µ, because B s → γµμ is experimentally most easily accessible at LHCb.

B q → γ ¯ amplitude
The B q → γ ¯ amplitude can be the parametrized in terms of four helicity amplitudes, where V and A refer to the vector and axial-vector chirality structure of the lepton currents, respectively. With the aid of (2.5) and (2.13) the helicity amplitudes are given by with LP and NLP contributions according to (3.2). The axial-vector amplitudes A hA depend only on a single Wilson coefficient, C 10 . The transversity amplitudes have definite CP transformation properties facilitating the analysis of the time-dependence in neutral B-meson decays. The amplitude of the CP-conjugated decay A(B q → γ ¯ ) is given by replacing in (4.1) with all complex-valued fundamental couplings g (Wilson coefficients and CKM elements) complex conjugated. The dependence on the phase of the CP transformation CP B q = e iξ Bq B q of the B q -meson state will cancel in observables.

Amplitudes at LP
We start with the LP contributions to the amplitudes A LV and A LA . (Recall that A RV = A RA = 0 at LP.) In the evaluation of the LP amplitude (3.30) we drop systematically higher-order terms in α s than NLO. This concerns the weak-EFT Wilson coefficients, the hard functions and the SCET I evolution matrix U H (µ h , µ hc ) in (3.20), the static HQET decay constant F Bq (µ hc ) in (3.25), and the convolution of the jet function with the B-meson LCDA φ + (ω; µ hc ). Schematically, the NLO correction reads where each quantity X = X (0) + X (1) + O(α 2 s ) is expanded in α s at its corresponding scale. In case of the evolution factors U , " (0) " means LL instead of NLL.

JHEP12(2020)148
Previous studies of B → (X s , K, K * ) + ¯ [25,26,38,39] have shown that NLO QCD corrections to these decays are sizeable, especially from the four-quark operators P c 1,2 . Including them leads to a significant reduction of the hard renormalization scale uncertainty. We therefore first look at the scale uncertainties at LO and NLO of our calculation of the B q → γ ¯ mode. Here the LO approximation also implies only LL evolution in U H and U F , and the omission of O(α s ) corrections to the hard and jet functions.
The reduction of the dependence on µ h = {2.5, 5.0, 10.0} GeV upon going from LO to NLO is shown in figure 4 for A LV and A LA . The NLO QCD corrections are sizeable for V eff 9 and V eff 10 , where LO and NLO scale uncertainty bands do not overlap. These corrections have been neglected in previous predictions of B q → γ ¯ . Whereas A LA is real-valued and slowly varying over q 2 ∈ [4m 2 , 8.0] GeV 2 , A LV has a sizeable imaginary part already at LO QCD because of the B-type contribution i = 7Beff, in contradistinction to the decays B → (X s , K, K * ) + ¯ , where such contributions are suppressed by the small QCDpenguin coefficients. 12 Further, Re(A LV ) is dominated for q 2 2 GeV 2 by the photon pole in i = 7Aeff, 7Beff, which interferes destructively with the i = 9eff contribution for q 2 2 GeV 2 and leads to a zero crossing around q 2 0 ≈ (3.0 − 3.5) GeV 2 . The zero crossing is responsible for the sign flip of the lepton forward-backward asymmetry A FB (q 2 ), see section 4.2. NLO QCD corrections shift q 2 0 to slightly larger values and reduce the µ hdependence. Moreover, in the q 2 -region of the zero-crossing of the LP contribution, the total amplitude is also very sensitive to NLP corrections, see section 4.1.2.
The dependence of A Lχ on the hard-collinear scale µ hc is also shown in figure 4, varying µ hc = {1.0, 1.5, 2.0} GeV. Also here a reduction of the scale uncertainty can be observed after including the NLO QCD corrections. Compared to the µ h dependence, the residual µ hc dependence at NLO is smaller in A LV and slightly larger in A LA , but both scale dependencies are small compared to other theoretical uncertainties.
The charm-quark pair threshold causes a divergence of A LV in figure 4 at q 2 = 4(m PS c ) 2 ≈ 7.73 GeV 2 from the two-loop matrix elements of P c 1,2 , which signals a breakdown of factorization and restricts us to q 2 6 GeV 2 . The uncertainties due to the charm-quark mass on Re A LV (Im A LV ) for q 2 6 GeV 2 are at most 1.0% (2.0%) in the PS scheme. Note that the physical threshold is at q 2 = m 2 J/ψ > 4(m PS c ) 2 . One might therefore be tempted to adopt the pole scheme for the charm-quark mass, in which case the partonic threshold divergence occurs closer to the physical threshold. However, we find that the uncertainty on Re A LV (Im A LV ) from m c for q 2 6 GeV 2 is 5.0% (5.0%) in the pole scheme, much larger than in the PS scheme, due to the larger uncertainty in the mass value.
The variation of the bottom-quark mass shows that for q 2 6 GeV 2 the uncertainty on Re A LV (Im A LV ) is less than 0.1% (0.3%) in the PS scheme and less than 0.3% (1.0%) in the pole scheme, and similarly for Re A LA .
As already mentioned, a strong dependence of A Lχ (χ = V, A) on the first inverse moment λ Bs of the B-meson LCDA φ + (ω) should be expected, which is shown in figure 5. It represents the largest uncertainty compared to the tiny µ h and the µ hc uncertainties, and 12 In the case of charged b → q ¯ decays also charged-current operators P u 1,2 contribute to these so-called weak annihilation contributions [25,26]. They are CKM suppressed for q = s, but not negligible for q = d, as for example in B + → M + ¯ , M = ρ, π [26,64].

JHEP12(2020)148
also affects the location of the zero-crossing q 2 0 of Re A LV . To estimate the further modeldependence on the B-meson LCDA, we vary also the first (hatted) logarithmic moment σ (1) Bs = (0.0 ± 0.7), as described in more detail in appendix C. The uncertainty due to σ (1) Bs is rather large for the contribution of i = 7Beff, but overall small compared to the λ Bs uncertainty.

Amplitudes at NLP
The NLP corrections to the amplitudes are due to local and nonlocal A-type and B-type, as well as the local four-quark contributions in (3.42), see section 3.2 and figure 3.
The local and nonlocal A-type NLP contributions to A hV (h = L, R) and A LA in (3.35)-(3.38) due to the semileptonic (i = 9eff, 10) and dipole (i = 7Aeff) operators as well as the local contributions from the four-quark (4qu) operators in (3.42) are compared in figure 6 separately to the LP results. For the non-local form factor parameterization we use the central value r LP = 0.2. It can be seen that the contributions i = 7Aeff, 4qu have a photon pole for q 2 2 GeV 2 and that i = 4qu contributes an imaginary part to the leftand right-helicity amplitudes. The NLP corrections to A LA are within the expected size of a Λ QCD /m b ∼ (10−15)% correction relative to LP. In the case of A LV , this applies individually to the real part of i = 7Aeff and i = 9eff, amounting to 15% and (10 − 15)% relative to their respective LP contributions, whereas their imaginary parts of the LP and NLP contributions are tiny (9eff) or vanish altogether (7Aeff). The real and imaginary parts of the local NLP contribution i = 4qu turn out to be rather large, constituting +(20 − 25)% and −30% of the real part of the LP i = 7Aeff contribution, which has been chosen for comparison because of a similar photon pole. The NLP corrections are most important in the q 2 region where the LP contributions cancel and give rise to the zero crossing of Re(A LV ). In fact, the location of the zero, q 2 0 ≈ 4 GeV 2 is significantly shifted and the sum of these NLP corrections is sizeable compared to the LP part for q 2 ∈ [3, 5] GeV 2 around the zero crossing q 2 ∈ [3, 5] GeV 2 , especially for Re(A LV ).
The right-handed helicity amplitudes A Rχ (χ = V, A) are entirely NLP. A RV exhibits a zero at q 2 0 ≈ 3 GeV 2 due to the interference of the i = 4qu, 9eff contributions, whereas the i = 7Aeff part is not enhanced at q 2 → 0, see (3.38). The amplitude A RA is small.   over V in (3.57). We assume the unknown decay constants and form factors of the higher φ, ρ and ω resonances to be a fraction

JHEP12(2020)148
of those of the lowest ones and use the same value r q f T = 0.3 for all higher resonances V > V . Since the higher resonances lie in the range of q 2 , where we might expect the continuum contribution calculated in factorization to already be dual to the contribution from the relatively broad resonances, the value of r LP should be increased (implying a smaller continuum part in the sum of LP and NLP according to (3.57)). Hence, we introduce separate parameters r A LP and r B LP , where the first is always kept at r A LP = 0.2 ± 0.2, but the second should depart from this value in the "3(4) Res" model, at least for q = s, where the effect of resonances is more pronounced.
In imaginary parts. The zero crossing of the real part for q = s is shifted to q 2 0 ≈ 4 GeV 2 when compared to the LP result of A LV in figure 4. For q = d, the zero-crossing occurs at the lower value q 2 0 ≈ 2.5 GeV 2 .

B q → γ ¯ observables
From the amplitude (4.1) we obtain the two-fold B q → γ ¯ decay-rate distribution in terms of the kinematic variables q 2 ∈ [4m 2 , m 2 Bq ] and cos θ ∈ [−1, 1], where θ is the angle between the B q -meson direction and the lepton momentum in the dilepton center-ofmass frame. The q 2 -dependent angular coefficients are expressed in terms of the transver-JHEP12(2020)148 sity amplitudes as follows: and Eq. (4.7) refers to the square of the so-called structure-dependent (SD) amplitude. We provide the lepton-mass suppressed parts of the double-differential width due to the FSR amplitude of P 10 and its interferences with the SD amplitude in appendix B. The differential decay width and the normalized lepton forward-backward asymmetry are where the dots denote lepton-mass suppressed terms from the FSR contribution. The angular observables a(q 2 ) and c(q 2 ) differ only by lepton-mass effects, such that their difference is lepton-mass suppressed. The lepton-mass suppression renders the experimental measurement of these parts of the two-fold differential decay width challenging. Moreover, the lepton-mass suppressed FSR contributions in appendix B cannot be disentangled from them and introduce a dependence on C 10 , which is absent in a(q 2 )−c(q 2 ). Further, the FSR JHEP12(2020)148 contribution introduces a nonanalytic cos θ dependence as can be seen in (B.6) and (B.7). In consequence there are the two main observables, dΓ/dq 2 and A FB (q 2 ), in the part of the phase space where 4m 2 q 2 . The two-fold decay-rate distribution for the CP-conjugated decay B q → γ ¯ is given by the quantities a, b and c without bars, which are obtained from (4.8)-(4.10) by using the CP-conjugated transversity amplitudes A ⊥χ and A χ (χ = V, A). The CP-transformation properties (4.4) of A ⊥χ imply a minus sign in front of b(q 2 ) in correspondence with other B → V ¯ decays (V = K * , φ, . . .) [65]. The combination of the two-fold decay-rate distributions of the decay (4.7) and the CPconjugated decay (4.14) allows for two CP-averaged and two CP-asymmetric observables. We consider the CP-averaged rate and the normalized lepton-forward-backward asymmetry as well as the CP rate asymmetry, defined by We do not study the CP-asymmetry of the normalized lepton-forward-backward asymmetry defined as . We note that an untagged sample actually depends on the CP asymmetry of the lepton forward-backward asymmetry A CP . On the other hand the measurement of the lepton forward-backward asymmetry The q 2 dependencies of the branching fraction dB/dq 2 , A FB (q 2 ) and A CP (q 2 ) for B s → γµμ and B d → γµμ are shown in figure 9 when including successively the various higherorder QCD and NLP corrections discussed in section 3. We recall that the calculation is valid in the "low-q 2 " region q 2 6 GeV 2 or equivalently E γ 2.1 GeV (y 0.8). In the figure, the NLP corrections have been separated into 1) purely local ("loc"), 2) nonlocal A-type ("A") as given in (3.56) and 3) nonlocal B-type ("B") as given in (3.57). At LP, the NLO QCD corrections -see "LP, LO" vs. "LP, NLO" -decrease (increase) dB/dq 2 for q 2 3.0 GeV 2 (q 2 3.0 GeV 2 ) and shift the position of the zero-crossing of A FB (q 2 ) by about ∆q 2 0 ≈ +0.5 GeV 2 towards higher values. Further, the NLO QCD corrections substantially change the CP asymmetry A CP (q 2 ). The CP asymmetry is always tiny for q = s, because it is suppressed by the Cabibbo angle, but for q = d it can be reach −10% locally below the zero-crossing at q 2 3.5 GeV 2 .
The local NLP corrections lead to a decrease of about (15 − 20)% of dB/dq 2 for q 2 3.0 GeV 2 , i.e. in the region where the individual LP contributions i = 7Aeff, 7Beff, 9eff cancel, and give a large shift ∆q 2 0 ≈ +0.8 GeV 2 to the zero of the A FB (q 2 ). The inclusion of A-type nonlocal corrections ("NLP loc + A") almost cancels the effect of "NLP loc" in dB/dq 2 , but not so in A FB (q 2 ) and A CP (q 2 ) for the chosen input r q f T = 0.3 and r LP = 0.2, leaving A CP (q 2 ) for q = d at about −20% at very low q 2 and +5% at higher q 2 . The NLP nonlocal B-type corrections for the default model with only lowest resonances affects only the very-low q 2 region, with strong impact from the φ(1020) on dB/dq 2 for q = s.   Table 3. The q 2 -integrated branching fractions 10 9 × B(B s → γµμ) and 10 11 × B(B d → γµμ) in four / five bins in the various approximations shown in figure 9. "NLP all" corresponds to "NLP loc + A + B" in figure 9. The last four columns list the uncertainties due to 1) renormalization scale uncertainties, 2) B-meson LCDA parameters and 3) the continuum fraction r LP as well as "total" when adding them in quadrature.
The lowest resonances also affect the A FB (q 2 ) and A CP (q 2 ) locally in q 2 . In total, the zero-crossing q 2 0 is significantly increased by NLP contributions compared to "LP, NLO". The difference in q 2 0 between B s and B d decays of about 0.5 GeV 2 is due to a different value of λ Bq , but also non-negligible contributions ∝ λ (q) u for q = d. At very low q 2 the distributions shown in figure 9 are strongly affected locally by the resonances from the B-type contributions, and therefore prone to the uncertainty in modelling the corresponding form factor. In the spirit of global parton-hadron duality, we investigate whether such effects become averaged once integrated over sufficiently large bins in q 2 . The q 2 bins are chosen in the low-q 2 region, with upper bin boundary q 2 max = 6 GeV 2 to avoid large impacts of charmonium resonances. 13 The lower bin boundaries are chosen as q 2 min = {4m 2 µ , 1.0, 2.0, 3.0, 4.0} GeV 2 , excluding q 2 min = 1.0 GeV 2 for q = s because it is at the center of the φ(1020) peak region. 14 The results for B(B s → γµμ) and B(B d → γµμ) are listed in table 3 for the same approximations of the various curves as in figure 9.
The block of the first three corrections in table 3, "LP, LO" through "NLP, loc", are found within systematic approximations and assumptions, and can be considered as model-independent. Omitting for a moment the nonlocal B-type contributions and com- 13 The bin q 2 ∈ [4m 2 µ , 8.64] GeV 2 is shown for comparison with [10]. 14 See [11] for predictions in this q 2 bin.

JHEP12(2020)148
paring "NLP loc + A" in the first two bins of B(B s → γµμ) shows that the photon pole at q 2 → 4m 2 µ in the bin q 2 ∈ [4m 2 µ , 2.0 GeV 2 ] contributes almost 90% of the rate. In this approximation the rate is of the order 4 · 10 −9 , comparable to the branching fraction B(B s → µμ) of the non-radiative mode. When adding the nonlocal B-type contributions, however, the φ(1020) contribution dominates over these "short-distance" corrections and enhances the rate by a factor of three to 12 · 10 −9 , see "NLP all", and is responsible for about 70% of the signal in the bin [4m 2 µ , 6.0 GeV 2 ]. The photon pole constitutes also about 80% of the rate of B(B d → γµμ) in q 2 ∈ [4m 2 µ , 6.0 GeV 2 ] when omitting the nonlocal Btype contributions due to the ρ and ω resonances, which however have only a small impact compared to the φ in B(B s → γµμ), as expected from the discussion of duality violation in section 3.3. They enhance the rate by a factor 1.2 to about 1.8 · 10 −10 , which is comparable to the SM prediction of B(B d → µµ) ≈ 1 · 10 −10 . The comparison of the first bin [4m 2 µ , 6.0 GeV 2 ] to the second bin shows that the branching fractions are reduced by a factor of about 40 and 8 for q = s and q = d, respectively, once the photon pole contribution is excluded, making them an order of magnitude smaller than for the corresponding non-radiative rare decays, B q → µμ. Note, however, that for = e the B q → eē rates are strongly helicity suppressed by (m e /m µ ) 2 , whereas the rates of B q → γeē are identical to B q → γµμ in the SM up to phase-space effects for very low q 2 .
Comparison of the value B(B s → γµμ) = (8.4±1.3)·10 −9 in the q 2 bin [4m 2 µ , 8.64 GeV 2 ] given in [10] with table 3 shows that our prediction is about a factor 1.5 larger and within the given errors there is some tension. Here we used the QCD factorization approach to compute the A-and B-type contributions, allowing us also to include at LP the NLO QCD corrections, contrary to [10]. In addition, we include in our model of the nonlocal NLP A-and B-type contributions also a continuum contribution. On the other hand we do not include a model for charmonium resonances (but rather stay away from them by restricting q 2 6 GeV 2 ), which could be responsible for some differences when q 2 approaches 8 GeV 2 . Similar comments apply to [11] who predict B(B s → γµμ) = 7.79 · 10 −9 in the bin q 2 ∈ [4m 2 µ , 6.0 GeV 2 ] and B(B d → γµμ) = (1.02 ± 0.16) · 10 −11 in the bin q 2 ∈ [1.0, 6.0] GeV 2 , roughly a factor of 1.6 and 2.3, respectively, smaller than our results in table 3. We attribute these differences to the difference between the QCD factorization computation of the B → γ * form factors and the more model-dependent approaches used in [10,11], and differences in the numerical values of the T Bq→V 1 form factors in the parameterization of the nonlocal NLP B-type contribution relative to [11].
We calculate the theoretical uncertainties from 1) scale variation, 2) B-meson LCDA parameters and 3) the modelling of the nonlocal A-and B-type contributions, as discussed for the amplitudes in section 4.1. They are listed for the best approximation "NLP, all" in table 3 and shown in figure 10 for the q 2 -distributions dB/dq 2 , A FB (q 2 ) and A CP (q 2 ). The default 1(2)-Res model is used and the uncertainties are added successively in quadrature.
The largest uncertainty is due to the LCDA parameters, amounting to about +70 −30 % for B s → γµμ in the bin q 2 ∈ [2.0, 6.0] GeV 2 (see table 3 moment of the B-meson LCDA, λ Bq , since roughly the LP is proportional to 1/λ Bq . 15 It is even larger for B d → γµμ in that q 2 bin, because here the relative variation of λ B d = (350± 150) MeV is larger and the minimal value in the variation is 200 MeV to be compared to λ Bs = (400±150) MeV. The variation r LP = 0.2±0.2 of the continuum fraction in the form factors model for the nonlocal A-and B-type contributions contributes another ±(35−45) % to the uncertainty when going to q 2 min 2.0 GeV 2 , beyond the region of the lowest resonances. The large uncertainty of the power-suppressed contributions does not come as a surprise given that the continuum is modelled as a fraction of up to 40% of the LP amplitude. 15 The dependence on the renormalization scales µ h,hc reduces when including NLO QCD corrections to the LP amplitude. At NLP there would be a strong uncancelled µ h dependence of m b C eff 7 in the model of the nonlocal A-and B-type contributions in (3.56) and (3.57). We neglect this particular scale dependence, which is a consequence of our chosen model, and find that the remaining µ h dependence of the local NLP corrections is small. The scale uncertainties is sub-leading, less than 15 %, compared to the parametric uncertainties. The ratios of observables with different lepton flavours provide interesting tests of lepton-flavour universality. The deviations from unity are tiny in the SM, and hadronic uncertainties cancel to large extent. Here we briefly discuss the ratio As shown in figure 12, these ratios take values between (0.98−1.00) for both q = s, d, except for q 2 0.4 GeV 2 . Within the framework of our calculation, lepton-mass dependence arises only from β through the phase-space measure, not from the amplitudes themselves. In table 4 we give the ratio for selected bins in q 2 . R qγ is always very close to unity except when q 2 min is below 1 GeV 2 . The considered sources of uncertainties cancel to very high degree, and in view of the JHEP12(2020)148  remaining tiny errors one should keep in mind that further not included higher order QCD and QED as well as NLP corrections can contribute at this level. So far we discussed the so-called instantaneous distributions at the time of the production of the neutral B q meson without accounting for B q -B q mixing. The time evolution due to mixing between production and decay further complicates the predictions of experimentally measurable observables. Thus, besides tagged and untagged there can be time-integrated or time-resolved measurements, all of which require a specific theoretical treatment of the time-dependence. In particular time-resolved measurements allow to measure effective lifetimes as well as direct and mixing-induced CP-asymmetries.
The modifications due to time-dependence might not be negligible for B s mesons, but are certainly small compared to the current parametric uncertainties. The B s → γ ¯ observable that is most likely to be measured first is the time-integrated branching ratio hence we focus on this to explore the phenomenological impacts of the B q -meson width difference. At the time t = 0 of the production of the B q meson the time-dependent decay rate dΓ[B q (t = 0) → γ ¯ ]/dq 2 = dΓ/dq 2 is given by the instantaneous one in (4.12) and correspondingly for the CP-conjugated decay dΓ[B q (t = 0) → γ ¯ ]/dq 2 = dΓ/dq 2 . In the case of B s → γ ¯ , the time-dependence can be included for the time-integrated branching fraction by the simple replacement in a(q 2 ) and c(q 2 ), where y s ≡ ∆Γ s /(2Γ s ). This holds in the SM where the tiny λ (s) u term can be neglected and no other weak phases than those in V tb V * ts contribute. The width difference of B s is already well measured to y s = 0.065 ± 0.003 [66]. For B d → γ ¯ , however, the λ These changes correspond to factors 1.006, 1.013, 1.024 and 1.027 with respect to the instantaneous branching ratios. The changes are even smaller than the O(y s ) ≈ ∓ 6.5% corrections of the individual transversity amplitudes squared (4.19) due to cancellation between the two in the branching fraction. Such a cancellation should be expected, since at LP the two transversity amplitudes are equal, and in this case the width difference affects the sum of the two amplitudes by the factor 1/(1 − y 2 s ), which deviates from unity only by O(y 2 s ). Overall, the width difference results in tiny modifications of the instantaneous rate compared to parametric and theoretical uncertainties. The time-dependence of other observables can be worked out as well if required, but we refrain from doing so here.

Summary and conclusions
In this paper we provided the first analysis of the rare radiative B q → γ ¯ decays in the kinematic regime of large photon energy with QCD factorization techniques, which construct a systematic expansion in inverse powers of large photon energy and heavy quark mass. As in the case of B u → γ ν , a two step matching on SCET I and SCET II leads to a systematic decoupling of hard and hard-collinear degrees of freedom. At leading power (LP) in 1/E γ , 1/m b this allows us to include the sizeable O(α s ) QCD corrections, which have previously been omitted. The neutral current radiative decays are more involved than B u → γ ν , since the real photon is not only emitted from the constituent quarks of the initial B-meson (A-type emission), but also directly from an operator in the weak EFT (B-type emission). Power-counting shows that the B-type emission also contributes at LP.
Further, we analyze the next-to-leading power (NLP) corrections and provide all local corrections, which can be calculated unambiguously. We also investigate the nonlocal corrections to A-and B-type emission from subleading-power form factors. A new feature compared to factorization calculations of B u → γ ν is the NLP form factor in the timelike region, which is affected by resonances in the region of very small dilepton mass, for which we adopted a parameterization in terms of resonances and a continuum term. We note that the well-known partial cancellation of the contributions from the operators P 9 and P 7 to the transversity/helicity amplitudes, which is also responsible for the zero at q 2 0 in the lepton-forward-backward asymmetry, makes the NLP corrections very important. This together with the resonance contamination in the very low q 2 region, the restriction JHEP12(2020)148 to q 2 6 GeV 2 to avoid the charmonium resonance, and the large uncertainty from the Bmeson LCDA parameters leads to the conclusion that precise calculations of rare radiative B q → γ ¯ decay observables are difficult in practice.
The analytical results allow us to provide predictions in the Standard Model for the q 2 -differential and q 2 -integrated CP-averaged branching fraction (dB/dq 2 ), lepton-forwardbackward asymmetry (A FB ), rate CP-asymmetry (A CP ), and lepton-flavour ratios. We quantified the dominant uncertainties from renormalization scales, the parametric dependencies from the B-meson light-cone distribution amplitude (LCDA) as well as the modelling of nonlocal NLP corrections. These studies show that the rate for the B s decay mode is strongly increased by the φ resonance due to the B-type contribution of the operator P 7 at very low q 2 2.0 GeV 2 , but becomes small once restricting to q 2 > 2.0 GeV 2 . Our main results are summarized in table 3  in two q 2 -bins, of which for q = s the first is dominated by resonances. In the case of B s → γµμ the B denotes the time-integrated branching fraction, which accounts for the non-vanishing width difference of the B s system. Compared to previous estimates [10,11], the branching fractions from the QCD factorization calculation performed here are roughly a factor of two larger. We attribute these differences to the difference between the QCD factorization computation of the form factors and the more model-dependent and less complete parameterizations used in previous work as well as different form factor input in the nonlocal NLP B-type contributions used in [11]. The presence of the additional real photon lifts the helicity suppression of the purely leptonic B q → ¯ , which yields similar rates for final states with electrons and muons. The rate CP asymmetry A CP for q = s is tiny, but reaches between −20% and +5% locally in q 2 for q = d.
Note added. When this paper was completed, ref. [67] appeared, which presents the first theoretical estimate of λ Bs . The QCD sum rule calculation for this quantity and λ Bs /λ B d is in very good agreement with the values used here.

A Definitions and conventions
Throughout we use the common definitions with the property p ν L µν V = 0 and p ν L µν A ∝ f Bq m /m Bq . In consequence the FSR contributions vanish except the one from P 10 , which is helicity suppressed and proportional to f Bq . The FSR amplitude then takes the form in terms of the kinematic variables (q 2 , cos θ ), where A FSR ≡ 4 C 10 Q f Bq . The additional contributions to the two-fold differential decay width from the so-called structure-dependent (SD) contributions in (4.7) consist of the interference of SD × FSR

JHEP12(2020)148
where for the CP-conjugated decay the transversity amplitudes in (B.6) should be replaced according to (4.4). Both are lepton-mass suppressed as can be seen by the overall factor (1 − β 2 ) and become important only for very low q 2 ∼ 4m 2 1 GeV 2 . Further they have a nonanalytic dependence on cos θ given by the factor (1 − β 2 cos 2 θ ) −n with n = 1, 2, respectively, which becomes most pronounced in regions of the phase space when q 2 4m 2 and hence β → 1. Then the collinear divergence in the limit cos θ → ±1 is regulated by the finite lepton mass m = 0. The nonanalytic cos θ dependence prevents a simple polynomial dependence on cos θ as given in (4.7), but due to the lepton-mass suppression, this is of concern only for small regions in the phase space, such that the phenomenologically most important observables without lepton-mass suppression remain the differential decay width and the lepton forward-backward asymmetry.

C B q -meson LCDA: RG evolution and model
The factorization approach leads to a convolution of the jet function with the nonperturbative B-meson LCDA φ + (ω; µ) in the LP amplitudes (3.30). In A-type insertions, at NLO in QCD, the convolution (3.28) is related to the inverse and first two logarithmic moments, whereas in B-type insertions the knowledge of the functional form of φ + (ω; µ) is required even at LO in QCD. Apart from theoretical constraints, little is known about φ + , because of the absence of stringent phenomenological constraints. In fact, the radiative chargedcurrent decay B u → γ ν in the same kinematical region as B q → γ ¯ is considered as the prime decay to determine the inverse (λ Bu ) and the first logarithmic (σ (1) Bu ) moments [20][21][22][23][24], and has been investigated by the Belle Collaboration [68]. There is a preference for small values of λ B u,d ∼ 0.2 GeV from hadronic two-body decays B → ππ, πρ, ρρ [69,70] in the framework of QCD factorization. QCD sum rules yield λ Bq ∼ 0.46(11) GeV [71].
Such values usually refer to a particular renormalization scale. Throughout we set µ 0 = 1 GeV as the initial scale, and RG evolution is used to evolve φ + to the hard-collinear scale µ hc , accounting for higher-order QCD corrections. The RG equations for φ + (ω; µ) and its moments involve a convolution in the momentum variable ω, which can be avoided at the leading-logarithmic order when going to dual (position) space [72], where the corresponding nonperturbative function η + (s; µ) = U + (s; µ, µ 0 )η + (s; µ 0 ) has autonomous scaling for each value of s.
We resort to the three-parameter model [24] φ + (ω; µ 0 ) = Γ(β) Γ(α) which permits to solve the transformation from dual to momentum space analytically, and in consequence also an analytic solution of the RG equation. The three parameters ω 0 , α and β are assumed to determine φ + (ω) at the scale µ 0 and the RG evolution is performed as given in [24]. Further

JHEP12(2020)148
which allows to associate ω 0 to λ Bq , whereas α and β determine the hatted logarithmic moments defined as The three-parameter model reduces to the exponential model [42] for α = β, which depends only a single parameter λ Bq = ω 0 , the inverse moment.
In the numerical analysis, we choose as default the exponential model, i.e. α = β. This fixes the logarithmic moments, as for example σ (1) Bq = (0 ± 6) to estimate the further model-dependence, using the three-parameter model (C.1). We use the following (α, β) tuples together with ω 0 :