Rare FCNC radiative leptonic $B_{s,d}\to \gamma l^+l^-$ decays in the Standard Model

We revisit rare radiative leptonic decays $B_{s,d}\to \gamma e^+e^-$ and $B_{s,d}\to \gamma \mu^+\mu^-$ in the Standard Model and provide the updated estimates for various differential distributions (the branching ratios, the forward-backward asymmetry, and $R_{\mu/e}$, the ratio of the differential distribution for muons over electrons in the final state). The new ingredients of this work compared to the existing theoretical analyses are the following: (i) we calculate all $B_d\to\gamma$ and $B_s\to\gamma$ form factors induced by the vector, axial-vector, tensor and pseudotensor quark currents within the relativistic dispersion approach based on the constituent quark picture; (ii) we perform a detailed analysis of the charm-loop contributions to radiative leptonic decays: we obtain constraints imposed by the electromagnetic gauge invariance and discuss the existing ambiguities in the charmonia contributions.


INTRODUCTION
Rare radiative leptonic B (s,d) → γl + l − decays are one of the flavor-changing neutral current (FCNC) decays: at the quark level, they are induced by b → {s, d} quark transitions, which in the Standard Model (SM) are forbidden at tree level. Such transitions occur via penguin and box diagrams containing loops and thus lead to small branching ratios, of order 10 −8 -10 −10 [1]. Possible contributions of new particles to the loops make these decays particularly sensitive to potential New Physics.
Obviously, other FCNC B-decay modes are good places to search for deviations from the SM. The focus of this paper is on rare radiative leptonic B d,s -decays.
The previous anylases employed various approximations for the decay amplitudes which excibit a rich structure of nonperturbative QCD effects. This work improves the existing analysis in several aspects: (i) We calculate all necessary B (s) → γ form factors at timelike momentum transfers using the dispersion formulation of the relativistic constituent quark model [33][34][35]. This approach proved to be very successful for the calculation of numerous meson-to-meson weak transition form factors [36]; in this work we apply this approach to the calculation of the B → γ transition form factors, taking into account rigorous constraints on the transition amplitude imposed by electromagnetic gauge invariance.
(ii) We derive the general gauge-invariance constraints on the charm-loop contributions to the B → γl + l − amplitude. We then perform a numerical analysis of charm-loop contributions in B → γl + l − decays, including nonfactorizable corrections, making use of the existing results for the B → V l + l − amplitude.
(iii) We present a detailed study of a number of observables in B s,d → γl + l − decays (the differential distributions, the forward-backward asymmetry, and R µ/e , the ratio of the differential distributions for muons over electrons in the final state, which has been recently emphasized in [30] as an interesting observable for radiative leptonic decays).
The paper is organized as follows: Section 2 briefly recalls the effective Hamiltonian for FCNC b → s, d transitions, and Section 3 describes various contributions to the B → γl + l − amplitude induced by H eff (b → (s, d)l + l − ). In Section 4, we discuss constraints on the B → γ transition amplitude imposed by electromagnetic gauge invariance. In Section 5, we study in detail contributions to the amplitude of radiative leptonic decay B → γl + l − induced by c-quark loops, including nonfactorizable effects, and derive rigorous constraints on these contributions imposed by electromagnetic gauge invariance. Section 6 recalls the differential distributions in radiative leptonic decays. Section 7 presents the analytic results for the transition form-factors within the dispersion approach based on constituent quark picture for all necessary B → γ form factors. Section 8 contains the numerical predictions for the necessary form factors and the observables. Section 9 summarizes our results.

THE b → d, s EFFECTIVE HAMILTONIAN
A standard theoretical framework for the description of the FCNC b → q (q = s, d) transitions is provided by the Wilson OPE: the b → q effective Hamiltonian describing dynamics at the scale µ, appropriate for B-decays, reads [37][38][39] [we use the sign convention for the effective Hamiltonian and the Wilson coefficients adopted in [40,41]].

1)
G F is the Fermi constant. The basis operators O b→q i (µ) contain only light degrees of freedom (u, d, s, c, and b-quarks, leptons, photons and gluons); the heavy degrees of freedom of the SM (W , Z, and t-quark) are integrated out and their contributions are encoded in the Wilson coefficients C i (µ). The light degrees of freedom remain dynamical and the corresponding diagrams containing these particles in the loops -in our case virtual c and u quarks -should be calculated and added to the diagrams generated by the effective Hamiltonian.
Necessary for theB s → γl + l − decays of interest are the following terms in (2.1) 2 (theB d → γl + l − case is obtained with the obvious replacement s → d) [26]: 2) The C 7γ part in Eq. (2.2) emerges from the diagrams in Fig. 1a,c with the virtual photon emitted from the penguin Notice that the sign of the b → dγ effective Hamiltonian (2.3) correlates with the sign of the electromagnetic vertex. For a fermion with the electric charge Q q e, we use in the Feynman diagrams the vertex iQ q eqγ µ qǫ µ .
As already noticed, the light degrees of freedom remain dynamical and their contributions should be taken into account separately. The relevant terms in H b→s eff are those containing four-quark operators:  [38][39][40][41][42]. In this section, we present the contributions to the B → γl + l − amplitude induced by operators (2.2) and (2.3) [26]. The B s → γ * transition form factors of the corresponding basis operators are defined as [25] . We treat the form factors as functions of two variables, F i (k ′2 , k 2 ): here k ′ is the momentum emitted from the FCNC b → q vertex, and k is the momentum of the photon emitted from the valence quark of the B-meson. The constraints on the form factors imposed by gauge invariance are discussed in Sect. 4.
A. Direct emission of the real photon from valence quarks of the B meson We denote as A (1) the contribution to theB s → γl + l − amplitude, induced by H b→sl + l − eff : the real photon is directly emitted from the valence b or s quarks, and the l + l − pair is coupled to the FCNC vertex (see diagrams of Fig. 1). It corresponds to the momenta k ′ = q, k = p − q, k ′2 = q 2 and k 2 = 0, and thus involves the form factors F i (q 2 , 0): with B. Direct emission of the virtual photon from valence quarks of the B meson Another contribution to the amplitude, A (2) , describes the process when the real photon is emitted from the penguin FCNC vertex, whereas the virtual photon is emitted from the valence quarks of the B-meson (Fig. 2). The amplitude A (2) has the same Lorentz structure as the C 7γ part of A (1) where now k = q, k ′ = p − q, k ′2 = 0 and k 2 = q 2 . The amplitude thus involves the form factors Fig. 3 gives diagrams for A Brems , the Bremsstrahlung contribution to theB s → γl + l − amplitude: Let us emphasize (see [26]) that the contribution of the operator O 9V to the Bremsstrahlung amplitude vanishes.

CONSTRAINTS ON THE TRANSITION FORM FACTORS
We now discuss the requirements imposed by the electromagnetic gauge invariance on the γ * |qO i b|B q (p) transition amplitudes induced by the vector, axial-vector, tensor, and pseudotensor weak currents. This discussion extends the discussion of [25] and includes also the case when the real photon is emitted from the FCNC b → q vertex. The corresponding form factors are functions of two variables, k ′2 and k 2 , where k ′ is the momentum of the weak b → q current, and k is the momentum of the electromagnetic current, p = k + k ′ . Gauge invariance provides constraints on some of the form factors describing the transition of B q to the real photon emitted directly from the quark line, i.e. for the form factors at k 2 = 0.
These form factors fully determine the amplitudes of the FCNC B-decays into leptons in the final state. For instance, the four-lepton decay of the B meson requires the form factors f i (k ′2 , k ′2 ) for 0 < k 2 , k ′2 < M 2 B . For the case of the B → γl + l − transition one needs the form factors f i (k ′2 = q 2 , k ′2 = 0) and f i (k ′2 = 0, k ′2 = q 2 ), where q is the momentum of the l + l − pair.

A. Form factors of the vector weak current
In case of the vector FCNC current, the gauge-invariant amplitude contains one form factor g(k ′2 , k 2 ): The amplitude is automatically transverse and is free of the kinematic singularities so no constraints on g(k ′2 , k 2 ) emerge.

B. Form factors of the axial-vector weak current
For the axial-vector current, the corresponding amplitude has three independent gauge-invariant structures and three form factors, and in addition has the contact term which is fully determined by the conservation of the electromagnetic current, ∂ µ j e.m. µ = 0: Here QB q = Q b − Q q is the electric charge of theB q meson and fB q > 0 is defined according to The kinematical singularity in the projectors at k 2 = 0 should not be the singularity of the amplitude, and therefore gauge invariance yields the following relation between the form factors at k 2 = 0: For the neutralB d,s mesons, the contact term is absent and therefore the form factor a 1 should vanish at k 2 = 0, a 1 (k ′2 , k 2 = 0) = 0. This relation is fulfilled automatically, as the two contributions, corresponding to the the photon emission from the valence b-quark and from the valence s, d-quark cancel each other at k 2 = 0. The amplitude of the transition to the real photon is described by a single form factor The transition amplitudes induced by the tensor weak current can be decomposed in the Lorentz structures transverse with respect to k α : The contact terms are absent in this amplitude as well as in the amplitude of the pseudotensor current. The kinematic singularity of the projectors at k 2 = 0 should not be the singularity of the amplitude, therefore Multiplying (4.6) by k ′ ν , we obtain the penguin transition amplitude Notice that the penguin amplitude contains only one combination of the form factors. Nevertheless, the requirement of the regularity of the amplitude (4.6) yields the constraint (4.7).

D. Form factors of the pseudotensor weak current
The transition amplitude of the pseudotensor weak current is given in terms of the same form factors as the amplitude (4.6), and, similar to (4.6), contains no contact terms: The kinematical singularity in the projectors at k 2 = 0 should cancel in the amplitude, again leading to the constraint Eq. (4.7). For the penguin pseudotensor amplitude we then obtain Notice that the contribution of the second Lorentz structure in (4.10) vanishes both for k 2 = 0 (because of the constraint Eq. (4.7): at k ′2 = 0, kp = kk ′ ) and for k ′2 = 0. However, it does not vanish for both k 2 , k ′2 = 0; therefore, the second Lorentz structure contributes to the amplitude of the four-lepton decays.
We can now build the bridge to the form factors which describe the real photon emission by the valence quarks defined in Eq. (3.1): denoting the momentum of the l + l − pair as q, i.e., setting k 2 = 0 and replacing k ′2 → q 2 , we obtain the form factors in Eq. (3.1) through the form factors g, a 2 , g 2 , g 1 (k ′2 = q 2 , k 2 = 0): The form factors describing the real photon emission from the penguin, are obtained by setting k ′2 = 0 and replacing k 2 → q 2 in the form factors g 1,2 (k ′2 , k 2 ): Let us notice that the form factor g 1 (q 2 , 0) should vanish at q 2 = M 2 B in order to kill the unphysical pole at q 2 = M 2

B
in the form factor F T A (q 2 , 0). We shall therefore perform an appropriate subtraction in the spectral representation for g 1 (q 2 , 0) to provide this property.

CHARM-LOOP CONTRIBUTIONS TO THE AMPLITUDE
Whereas heavy degrees of freedom (t, W , Z) have been integrated out when constructing the effective Hamiltonian for b-decays, light degrees of freedom, in particular c and u quarks, remain dynamical and their contributions in the loops should be taken into account separately.
We consider in this section the charm-loop contributions to the B s → γl + l − amplitude, which are related to the following matrix element: Here the quark fields are the Heisenberg operators in the SM, i.e. the corresponding S-matrix includes weak interactions of quarks. The matrix element (5.1) has the form dictated by the conservation of the vector charm-quark and the electromagnetic currents that requires k α H µα (k ′ , k) = 0 and k ′µ H µα (k ′ , k) = 0 (notice the absence of any contact terms): with the invariant form factors H i depending on two variables, H i (k ′2 , k 2 ). The singularities in the projectors at k 2 = 0 and k ′2 = 0 should not be the singularities of the amplitude H µα (k ′ , k), leading to the constraints Let us show that H 3 does not contribute to the B s → γl + l − amplitude: to obtain the latter, H µα should be multiplied by either ǫ α (k)lγ µ l or ǫ µ (k ′ )lγ α l. In each case, those terms in the H 3 -part of H µα containing k ′ µ or k ′ α vanish in the B s → γl + l − amplitude; the contribution of the "regular" structure k α k µ also vanishes because the form factor H 3 = 0 if k 2 = 0 or k ′2 = 0. (The situation is different for the transition into four leptons via two virtual photons, in which case the H 3 structure also contributes to the B 0 → l + l − l + l − transition amplitude). Now, let us consider the matrix element (5.1) at the lowest order in the weak interaction. Fig. 4 shows the diagrams representing the charm contribution to the B → γ * γ * amplitude. The diagram of Fig. 4a is generated by the squark part of the electromagnetic current j e.m. In addition, the B → γl + l − amplitude receives contributions from similar diagrams with the c-quark replaced by the u-quark. The latter, however, contain the CKM factor V ub V * us ≪ V cb V * cs and are therefore strongly suppressed compared to the charm contribution.

A. Charming penguins
The analytic expression for the charming-penguin diagram of Fig. 4a has the form cc .
Similar to the diagrams discussed in the previous Section, the diagram of Fig. 4a generates the following two types of contributions to the B s → γl + l − amplitude shown in Fig. 5: the c-quark emits the virtual photon whereas the s-quark emits the real one (a) and the c-quark emits the real photon whereas the s-quark emits the virtual one (b). Diagrams of Fig. 5 lead to the following contributions to the B s → γl + l − amplitude: The Lorentz structure of the amplitudes A (1,2) cc coincides with the Lorentz structure of the amplitudes A (1,2) , therefore the full charm contribution can be described as additions to the invariant amplitudes in Eqs. (3.2) and (3.4): The challenging task in the analysis of the charm-loop contributions given by H i (q 2 , 0) is the necessity to describe a wide range 0 < q 2 < M 2 B , including the region of charmonium resonances. Perturbative QCD cannot be applied here and nonperturbative approaches based on hadron degrees of freedom are necessary, see discussion in [42][43][44][45][46][47][48][49]. For H i (q 2 , 0) one may write dispersion representation in q 2 with two subtractions, similar to the B → K * l + l − amplitudes H 1,2 of [42]: where a i and b i are the (unknown) subtraction constants and the functions h i (q 2 ) describe the hadron continuum including the broad charmonium states lying above the DD threshold. The contribution to A (1) i given by the form factors H i (q 2 , 0) may be described as the correction to C 9V , i.e. by the . Obviously, this correction will be process-and Lorentz-structuredependent: it will in general be different for B → P l + l − decay and for B → V l + l − decay and different in the A V and A A amplitudes. Describing the nonfactorizable effects as a shift in C 9V is not particularly convenient in the region of small q 2 : whereas the full nonfactorizable correction amounts to a few percent at small q 2 [42], it explodes if expressed as a correction to the coefficient C 9V . Nevertheless, describing the charm-loop effects as corrections to the Wilson coefficients has one very important advantage: these corrections are obtained as the ratio of the functions H i (q 2 , 0) and the appropriate form factors (B → K * or B → γ). It is reasonable to expect that the effects related to the difference between the vector meson and the photon in the final state cancel to large extent in the ratios and that the corrections to the Wilson coefficients are approximately equal to each other for B → γl + l − and B → V l + l − . The accuracy of this approximation is expected to be at the level of 10%-20%, the typical accuracy of the vector meson dominance.
Similarly, the contributions to A (2) i given by H i (0, q 2 ) may be described as corrections to C 7γ : C 7γ → C eff 7γ (q 2 ) = C 7γ + ∆C 7γ (q 2 ). In principle, this correction is also non-universal and q 2 -dependent. However, the form factors H i (0, q 2 ) and F T i (0, q 2 ) have similar q 2 -dependences, as they contain contributions of the samess hadron resonances in the q 2 -channel. One therefore expects that the correction to the Wilson coefficient C γ (which is purely nonfactorizable, see the discussion below) may be taken q 2 -independent.
In view of these arguments, we will use the results from [42] for ∆C 9 and ∆C 7 obtained at low q 2 for our analysis.

Factorizable part of the amplitude Hµα
Before discussing the full charm-loop corrections to the Wilson coefficients, we present the results for factorizable contributions. Taking into account only factorizable gluon exchanges leads to where the expression in brackets is just the amplitude of (4.2) and For the invariant function Π cc (s) we may write the spectral representation with one subtraction The factorizable contributions to the form factors H i (k ′2 , k 2 ) are related to f, a 1 , a 2 as follows (5.14) Obviously, H fact V,A,3 (k ′2 , k 2 ) vanish for k ′2 = 0. 3 Therefore, the factorizablecc contribution to the amplitude A (2) and, respectively, to C 7γ , vanish; thecc contribution to A (2) comes exclusively from nonfactorizable gluon exchanges.
The factorizablecc contribution to A (1) can be described as a universal q 2 -addition to the coefficient C 9V :

Adding non-factorizable corrections to the amplitude Hµα
The functions H i (q 2 , 0) may be obtained at q 2 ≪ 4m 2 c using the method of QCD sum rules. The necessary calculation for the B → γl + l − amplitude are not available yet; however, in [42] the functions H i (q 2 , 0) were calculated for the B → K * l + l − amplitude. As already mentioned above, if the charm-loop effects are described as corrections to the Wilson coefficients, the latter are obtained as the ratio of the functions H i (q 2 , 0) and the appropriate B → V form factors (V = K * or V = γ). There are good reasons to expect that the effects related to the difference in the final states cancel to large extent in these ratios. Therefore, we will use the results for ∆C 9V (q 2 ) and ∆C 7γ (q 2 ) obtained at low q 2 from [42] for our analysis.
In [42] nonfactorizable corrections at low k 2 have been calculated using light-cone QCD sum rules. The authors noticed that in distinction to the positive-definite factorizable contributions, nonfactorizable corrections are not positive-definite and therefore different charmonium resonances may in principle appear with different signs. Recall that the absolute values of the amplitudes A i BψK * (cf. Eq. (5.7)) for ψ and ψ ′ are known from the experimental data on B → (ψ, ψ ′ )K * decays, but the phases are unknown.
From our point of view, no conclusion about the relative signs of the resonance contributions may be drawn from the results for ∆C 9V (q 2 ) obtained at q 2 ≤ 4m 2 c where the calculation is trustable. Following [42], we describe h(q 2 ) as an effective heavier resonance of zero width, h(q 2 ) = c/(M 2 R − q 2 ). The unknown parameters are now the subtraction constants a, b, and the parameters of the effective pole c and M 2 R . Fig. 6 presents two different fits to the results of [42] for ∆C 9V (q 2 ) at q 2 < 4 GeV 2 (for our analysis ∆C from [42] are relevant; within errors both are equal to each other so we take the same ∆C 9V in the vector and the axial-vector amplitudes): one fit assumes the standard same positive contributions of ψ and ψ ′ ; another fit assumes an opposite sign for the ψ ′ contribution. Obviously, even the knowledge of ∆C 9V (q 2 ) at q 2 < 4 GeV 2 with the accuracy of a few percent would not allow one to discriminate between the same-phase and the opposite-phase cases. Taking into account the expected uncertainty of about 30-50% of the results from QCD sum rules (see Fig. 5 of [42]), the question of the relative phases between ψ and ψ ′ remains fully open.
We would like to recall that the LHCb collaboration tested the charmonia contributions in the B → Kl + l− decays but did not manage to decide in favour of one or another phase assignments between ψ and ψ ′ . One should take into account, however, that in principle the pattern of ψ and ψ ′ signs may be different in B → Kl + l − and in B → (K * , γ)l + l − decays.
To facilitate calculations at q 2 > 4M 2 D , we take into the contributions of the known broad vector ψ n (n = 3, . . . , 6) resonances (and add a heavier effective pole) to h(q 2 ) in (5.7). This may be done by the following addition to ∆C 9V (q 2 ): .
Here, the factorizable contribution of each ψ n is multiplied by a fudge factor κ n (following the old way of taking into account nonfactorizable corrections [43]) and a q 2 -dependent Γ n (q 2 ) = n is used to enable using this expression also below the DD threshold. It seems reasonable to take κ n ∼ 1.5 − 2.0 for all excited vector charmonia: the experimental data for B → ψK * and B → ψ ′ K * lead to κ J/ψ = 1.6 and κ ψ ′ = 1.9. The subtraction constants a, b and the parameters of the effective pole are again fixed by requiring that ∆C 9V (q 2 ), corresponding to (5.7) with the addition (5.16), reproduces the sum-rule results at 0 < q 2 < 4 GeV 2 . And again, the latter may be easily reproduced with 1% accuracy.   [42] at 0 < q 2 < 4 GeV 2 . In the range 0 < q 2 < 4 GeV 2 both prescriptions for the resonance phases reproduce the results from QCD sum rules with better than 1% accuracy.
To conclude, the results from QCD sum rules available at q 2 ≤ 4m 2 c cannot give a conclusive answer of the relative signs of ψ and ψ ′ resonances (cf. [42]). We shall discuss later some observables which are particularly sensitive to the relative phases of the ψ and ψ ′ and could shed light on this issue if measured experimentally. Fig. 7 shows the typical weak-annihilation (WA) diagrams, which emerge from diagram Fig. 4 after integrating out the W -bosons and taking into account QCD radiative corrections. Diagrams of Fig. 7(a) and similar diagrams with gluon exchanges between quark from the same loop lead to factorizable contributions ∼ f Bs G γγ (k ′2 , k 2 ), where the form factor G γγ (k ′2 , k 2 ) does not depend on the B s -meson structure; the B S -meson contribution is reduced to a single quantity, f Bs . Diagrams of the type Fig. 7(b) are not reduced to f Bs but contain more complicated quantities describing the B s -meson structure. Diagrams with c-quarks replaced by u-quarks should be included but they are CKM suppressed compared to the charm-loop contributions. We denote as A WA the corresponding contribution of these diagrams to the B → γl + l − amplitude, and we take into account both c and u quarks in the loop. The vertex describing thebs →Ū U transition (U = u, c) reads

B. Weak annihilation
with a 1 = C 1 + C 2 /N c , N c number of colors [50]. For N c = 3 one finds a 1 = −0.13. We now have to take TheŪ U contribution to this amplitude can be written as where the form factor G(p 2 , k 2 , q 2 |m 2 U ) is defined as follows [51] γ For massless u-quark in the loop, axial anomaly [52,53] fixes the form factor G γγ (p 2 , k 2 , . For c-quark there is an additional q 2 -dependent contribution given by the amplitude m c γ * (k)γ * (q)|cγ 5 c)|0 ∼ m 2 c /M 2 B , which contains ψ and ψ ′ resonances at q 2 > 0. The latter contribution is numerically negligible compared to contributions discussed in the previous sections for all q 2 in the reaction of interest. Therefore, we have The WA contribution is enhanced at small q 2 , but even here it is suppressed by a power of a heavy quark mass compared to the contributions discussed in the previous sections [55]. 6. THE B → γl + l − DIFFERENTIAL DISTRIBUTION For convenience, we recall here the results from [26] for the differential distributions. The amplitudes discussed in Sections 3(A-C) have the same Lorentz structure, whereas the structure of the Bremsstrahlung amplitude (Section 3 C) is different. Therefore, it is convenient to write the cross-section as the sum of three contributions: square of the amplitude A 1+2+W A which we denote Γ (1) , square of the amplitude A Brems which we denote Γ (2) , and their mixing, denoted as Γ (12) (in this Section F V,A (q 2 ) stands for F A,V (q 2 , 0)): 4 In the above formulas the complex form factorsF T V,T A and defines as follows The expressions (6.2) and (6.3) contain the infrared pole which requires a cut in the energy of the emitted photon. Clearly, the contribution of the pole is proportional to the lepton mass.

FORM FACTORS FROM THE DISPERSION APPROACH BASED OF THE RELATIVISTIC CONSTITUENT QUARK PICTURE
So far our discussion was fully general. The problem we face now is to obtain the necessary form factors describing theB → γ transition. These form factors are very complicated objects which involve theB mesons in the initial state and require a treatment with nonperturbative QCD. In the previous section we obtained the constraints on the form factors coming from the electromagnetic gauge invariance of the amplitude. There is a number of general constraints on the form factors which emerge within large-energy effective theory (LEET) [56,57]. In our previous analysis we have made use of a simple model for the form factors based on LEET and the location of meson singularities in the corresponding channels. The SU (3) breaking effects in the form factors have been neglected in that work.
We now improve our predictions calculating the form factors using the dispersion formulation of the relativistic quark model. In this paper we present a model calculation of all necessary form factor within the relativistic dispersion approach based on the constituent quark picture. This approach has been formulated in detail in [33,34] and applied to the weak decays of heavy mesons in [36].
The pseudoscalar meson is described in the dispersion approach by the vertexq The pseudoscalar-meson wave function φ P is normalized according to The decay constant is represented through φ P (s) by the spectral integral The vector meson is described in the dispersion approach by the vertexq 2  The vector-meson decay constant is represented through φ V (s) by the spectral integral [58] f The photon of virtuality k 2 is described by The P → γ form factors f i (k ′2 , k 2 ) defined in the previous Section are obtained as double spectral representations in terms of the relativistic wave function of the B-meson in the form The double spectral representation (7.6) corresponds to considering the double cut of the triangle diagram in variables p 2 and k 2 , treating the variable k ′2 as the fixed external current virtuality. The double spectral densities ∆ i (s 1 , s 2 , k ′2 |m 2 , m 1 , m 1 ) in variables p 2 and k 2 may be obtained from the known spectral densities of the P → V transition form factors given by Eqs. There is, however, also another possibility to obtain the double spectral representation for f i (k ′2 , k 2 ): one can consider the double cut of the triangle diagram in p 2 and k ′2 , at a fixed value of k 2 .
f i (k ′2 , k 2 ) = ds 1 φ P (s 1 )ds 2 1 s 2 − k ′2∆ i (s 1 , s 2 , k 2 |m 1 , m 1 , m 2 ). (7.7) The double spectral densities∆ i (s 1 , s 2 , k 2 |m 1 , m 1 , m 2 ) differ from ∆ i (s 1 , s 2 , k ′2 |m 2 , m 1 , m 1 ), but the form factors calculated from (7.6) and (7.7) are of course equal to each other. The benefit of using the spectral representations in the form (7.7) shows up when one considers the transition to the real photon, k 2 = 0: in this case, ∆ i (s 1 , s 2 , k 2 → 0|m 1 , m 1 , m 2 ) →ρ i (s 1 |m 1 , m 2 )δ(s 1 − s 2 ), and the double spectral representation (7.7) is reduced to the single dispersion representation; whereas (7.6) remains a double spectral representation also for k 2 = 0. Taking into account this property, we obtain and use the single spectral representations for the form factors F A,V,T A,T V (q 2 , 0), but employ the known double spectral representations in the form (7.6) for F T A,T V (0, q 2 ). We have checked that for F i (q 2 , 0) both representations give the same results.  The form factor F A describing theB in the initial state is given by the diagrams of Fig. 8. Fig. 8a shows F

(b)
A , the contribution to the form factor of the process when the b quark emits the photon; Fig. 8b describes the contribution of the process when the quark d emits the photon while b remains a spectator. It is convenient to change the direction of the quark line in the loop diagram of Fig. 8b. This is done by performing the charge conjugation of the matrix element and leads to a sign change for the γ ν γ 5 vertex. Now, both diagrams in Fig. 8a,b are reduced to the same diagram where quark 1 emits the photon and quark 2 is a spectator: setting For the form factor F A a single dispersion integral was obtained in [59]: where ρ + (q 2 , m 1 , m 2 ) = (m 2 − m 1 ) λ 1/2 (q 2 , m 2 1 , m 2 2 ) s + m 1 log q 2 + m 2 1 − m 2 2 + λ 1/2 (q 2 , m 2 1 , m 2 2 ) q 2 + m 2 1 − m 2 2 − λ 1/2 (q 2 , m 2 1 , m 2 2 ) , (7.10) . (7.11) Notice that this expression differs from the analogous expression from [59]: in the second term we have the factor 1/(s − q 2 ) instead of the factor 1/(M 2 B − q 2 ) in Eq. (3.7) of [59]. This corresponds to a slightly different subtraction prescription in the dispersion integral: the factor 1/(M 2 B − q 2 ) would lead to the appearance of the unphysical pole at q 2 = M 2 B . For the case of a fixed q 2 = M 2 V , considered in [59], both subtraction prescriptions lead to very close numerical results. FV (q 2 , 0) The consideration of the form factor F V is very similar to the form factor F A . F V is determined by the two diagrams shown in Fig. 9: Fig. 9a gives F Again, we change the direction of the quark line in the loop diagram of Fig. 9b by performing the charge conjugation of the matrix element. For the vector current γ ν in the vertex the sign does not change (in contrast to the γ ν γ 5 case considered above). Then the contribution of both diagrams in Fig. 9a,b are given in terms of the form factor F (1)

Form factor FT
The form factor F T A contains two contributions corresponding to the cases when the photon is emitted from b (Fig. 10a) and from d(s) (Fig. 10b) (7.14) For the form factor F (1) T A may be written in the form, see (4.12), For the form factors g 0,2 we obtained the following spectral representations (7.16) with the spectral density (7.17) and ρ k 2 ⊥ given by Eq. (7.11). A subtraction has been performed in g 0 to provide its vanishing at q 2 = M 2 B .

Form factor FT
The form factor F T V also contains two contributions corresponding to the cases when the photon is emitted from b or s quarks of theB s -meson. These contributions are shown in Fig. 11a and Fig. 11b. b.
Changing the direction of the quark line in the loop diagram of Fig. 11b by performing the charge conjugation of the matrix element, we reduce the contribution of each of the diagrams of Fig. 11a,b to the same form factor F (1) T V where quark 1 with mass m 1 emits the photon and quark 2 with mass m 2 remains spectator such that The form factor F (1) T V may be written in the form, see (4.12), (7.19) with the form factors g 0,2 given by Eq. (7.16).
The form factor F T A (0, q 2 ) = F T V (0, q 2 ), which we denote as F T (q 2 ), also has two contributions related to the virtual photon emission by the b and the d(s) valence quark. Recall, that the photon emission from the heavy valence quark is strongly suppressed compared to the photon emission from the light valence quark by a parameter m d,s /m b , where m d,s are the constituent light-quark masses.
The double spectral representations allow us to calculate the form factors at q 2 < 4m 2 q , where m q is the mass of the quark which emits the photon. For the case of the photon emission by the b-quark, the region q 2 < 4m 2 b fully covers the decay region 0 < q 2 < M 2 B . The form factor F (b) (0, q 2 ) in this region is the real function and can be reliably calculated within the dispersion approach.
The situation, however, changes when we consider the process with the virtual photon emission from the light, d-or s-quark: e.g., the physical form factor F (d) (0, q 2 ) has imaginary part at q 2 > 4m 2 π , and contains contribution of light neutral vector-meson resonances (ρ 0 and ω for B-decays and φ for B s -decays) in the physical q 2 -region. Obviously, our calculation based on the quark degrees of freedom is trustable (far) from the hadron thresholds and should not be applied in the region of hadron resonances.
To obtain the form factor F T (0, q 2 ) at 0 < q 2 < M 2 B we therefore proceed as follows: 5 we calculate the form factors F (d,s) (0, q 2 ) using the gauge-invariant version [60] of the vector meson dominance [61][62][63] where M V and Γ V are the mass and the width of the vector meson resonance, g B→V + (0) are the B → V transition form factors, defined according to the relations For the calculation of the B → V transition form factor g B→V + (0), we make use of the same dispersion approach of [33,34]. The e.m. leptonic decay constant of a vector meson is given by V . (7.21)

NUMERICAL RESULTS
A. Calculation of the transition form factors

Parameters of the model
The wave function φ B (s), can be written as with w(k 2 ) normalized as follows The meson weak transition form factors from dispersion approach reproduce correctly the structure of the heavy quark expansion in QCD for heavy-to-heavy and heavy-to-light meson transitions, as well as for the meson-photon 5 Another possibility would be to calculate the form factor F (b) (0, q 2 ) using the dispersion approach and to apply Vector-Meson Dominance (VMD) to F (d,s) (0, q 2 ). However, F (b) (0, q 2 ) is a relatively flat function at 0 < q 2 < M 2 B (it has pole at q 2 = M 2 Υ ), and provides numerically small contribution at the level of about 10% compared to F (d,s) (0, q 2 ). We therefore find eligible to apply the VMD approximation to the full form factor F (d,s) (0, q 2 ). transitions, if the radial wave functions w(k 2 ) are localized in a region of the order of the confinement scale, k 2 ≤ Λ 2 [33,34].
Following [36,54], we make use of a simple Gaussian parameterization of the radial wave function which satisfies the localization requirement for β ≃ Λ QCD and proved to provide a reliable picture of a large family of the transition form factors [36].
In [36] we fixed the parameters of the quark model -constituent quark masses and the wave-function parameters β i of the Gaussian wave functions -by requiring that the dispersion approach reproduces (i) decay constants of pseudoscalar mesons and (ii) some of the well-measured lattice QCD results for the form factors at large q 2 . The analysis of [36] demonstrated that a simple Gaussian Ansatz for the radial wave functions allows one to reach this goal (to great extent due to the fact that the dispersion representations satisfy rigorous constrains from non-perturbative QCD in the heavy-quark limit). With these few model parameters, [36] gave predictions for a great number of weaktransition form factors in the full kinematical q 2 -region of weak decays. However, the analysis of [36] made use of some approximations which need to be updated: namely, the wave-function parameters of ρ 0 and ω mesons, β ρ and β ω , were assumed to be equal to each other, and β φ has been set equal to that of η meson.
In this paper, we make use of the same effective constituent quark masses as obtained in [36]  but make the following updates on the determination of the meson wave-function parameters: First, we make use of the recent results f B = 189(4) MeV and f Bs = 225(4) from lattice QCD [64] to fix the wave-function parameters of B and B s ; this new inputs lead to a more reliable calculation of the form factors F i (q 2 , 0) which involve only B-meson wave functions. Second, for the calculation of the form factors F i (0, q 2 ) based of vector-meson dominance, we need the transition form factor g + (0) for the transition of B s,d meson to neutral vector mesons ρ 0 , ω, φ. As mentioned, in [36] the wave-function parameters of light vector mesons have been set equal to the wave-function parameters of the corresponding pseudoscalar mesons. This turns out to be a rather crude approximation; we now improve the analysis and fix also the vector-meson wave-function parameter β from the reproduction of its decay constant. The new parameters turn out to be in some cases rather different from those used in [36].

Form factors Fi(q 2 , 0)
For the calculation of the form factors F i (q 2 , 0), we need the B (s) wave-function parameters, which we fix by the requirement to reproduce the latest results for the leptonic decay constants f B and f Bs [64]. The obtained wavefunction parameters and the corresponding decay constants of beauty mesons are quoted in Table 1.  With the wave function of B(B s ) meson fixed, we calculate the form factors F 0)) via the spectral representations given in the previous Section. Fig. 12 shows the results of our calculation.
Formally, the spectral representations allow one to calculate the form factors in the region q 2 < (m b + m q ) 2 . However, the results from an approach based on quark degrees of freedom should not be trusted in the region close to hadron thresholds. Therefore, one cannot guarantee our calculations to be reliable at q 2 ≥ 20 GeV 2 . On the other hand, we know that the form factors have poles at q 2 = M 2 R , where M R is the meson with the appropriate quantum  numbers: 1 − for F V and F T V ; 1 + for F A and F T A . Therefore, following [36] we parameterize our results by the "modified" pole function which explicitly takes into account the presence of the pole in the right location. Eq. (8.5) approximates the calculation results with better than 3% accuracy in a broad range 0 < q 2 < 25 GeV 2 . Table 2 gives the corresponding parameters for the B d,s → γ form factors. We also consider a single-pole parametrization of the form factors suggested by large-energy effective theory and used in [25]. At q 2 ≤ 15 GeV 2 both parametrizations agree with better than 10% accuracy, which is expected to be the typical systematic uncertainty of the form factors from the dispersion approach. At larger q 2 , the deviation between the two parametrizations increases, reaching ∼ 20% at q 2 ∼ 25 GeV 2 . We take this difference as a magnitude of the systematic uncertainty of our predictions at large q 2 .
We would like to comment on the systematic uncertainty of our predictions: the constituent quark picture for the form factors is an approximation to a very complicated picture in QCD. It is therefore not possible to provide any rigorous estimate of the systematic uncertainties of the obtained form factors. Based on the comparison of the results from our dispersion approach to the results from QCD sum rules and lattice QCD where such results are available, we expect the systematic uncertainty of our form factors at the level of 10% in the region 0 < q 2 < 15 GeV 2 and about 20% at 15 < q 2 GeV 2 .

Form factors FT
For the form factor F T A (0, q 2 ) = F T V (0, q 2 ), we make use of the VMD formula (7.20) which involves the B → V transition form factor and the decay constants f V . According to adopted procedure, we fix the wave-function parameter β for neutral vector mesons based on their leptonic decay constants, and then, having at hand the wave functions of the B s,d mesons, calculate the B → V transition form factors of interest.
As inputs for the decay constants of the vector mesons, we make use of two sets of values: one set of f V is obtained directly from the experimental data [65] neglecting the mixing effects; the second set makes use of the values obtained in [66] from the experimantal data including also mixing effects. Table 3 gives the wave function parameters of neutral vector mesons corresponding to these two sets of the decay constants. For the parameters β V we make use of the range, the boundary values of which yield the decay constants with/without mixing effects. The electromagnetic decay constants which enter the VMD formula (7.20) are related to f V as follows: f e.m.
Having fixed the wave-function parameters, we calculate the transition form factor g + using the double spectral representation for the P → V form factors given in [35]. The obtained results are summarized in Table 4 and compared with the determinations from other approaches. We point out a sizeable reduction of the B s → φ form factor compared to the result of [36]. This change just reflects the natural sensitivity of the form factor to the shape of the wave functions of the participating mesons. As already pointed out, in [36] β φ was assumed equal to β η , which turns out a rather crude approximation. Fixing β φ from the known value of f φ , as done in this work, is obviously much more reliable. As seen from Table 4, our present form-factor results, corresponding to the wave-function parameters fixed by using the decay constants as inputs, demonstrate the agreement with the latest results from light-cone sum rules, within the expected 10-15% uncertainty.
We take into account the contributions to the form factor F T A,T V (0, q 2 ) from the light ground-state vector mesons, neglecting those from the excited mesons. The latter are difficult to estimate directly, but the experience from the pion elastic and weak transition form factors in the timelike region suggests that the magnitude of the contributions of the excited states is numerically much smaller compared to the ground-state contributions [60]. Then, the contributions of the excited vector mesons to F T A,T V (0, q 2 ) may be neglected compared with other contributions to the B → γl + l − amplitude.  With all form factors known, we are in a position to calculate numerous differential distributions in B → γl + l− decays. The necessary inputs such as the V CKM and the quark masses are taken from [65]; the values of the Wilson coefficient are summarized at the end of Section 2. The only theoretical ingredient which still contains ambiguities is the contribution of charm-loops in the charmonia resonance region. As discussed above in Section 5, one can obtain an excellent description of the QCD sum rule results [42] available at low q 2 with different assignments of the resonance phases. The phase ambiguity cannot be resolved on the basis of the theoretical arguments only, but needs further inputs from the experimental measurements which will become available in the future. Most of the differential distributions presented below are obtained for the standard assignment of positive contributions of all charmonia, that follows the patters of the factorizable contributions. The only exception is the forward-backward asymmetry, A FB , in which case we discuss two different assignments and demonstrate a strong sensitivity of A FB in the q 2 -region between ψ and ψ ′ to the specific choice of the resonance relative signs.

The differential branching ratios
The differential branching ratios are shown in Fig. 13. The results in Fig. 13 correspond to the description of the charm-loop effects according to Eq. (5.7) and adding the contributions of the broad charmonia according to Eq. (5.16), and further assuming that all charmonia contribute with the same positive sign (coinciding with the sign of the factorizable contribution). The subtraction constants a and b in Eq. (5.7) are determined by the requirement to reproduce the known results at q 2 ≤ 4m 2 c , including nonfactorizable corrections calculated in [42]. As discussed in Section 5, this reproduction may be reached with an excellent few percent accuracy for the different assignments of the resonance phases thus leaving the question of the relative resonance phases open. In the region q 2 ≤ 6 GeV 2 , the charming loops provide a mild contribution at the level of a few percent, and therefore the branching fractions in this region may be predicted with a controlled accuracy, mainly limited by the form-factor uncertainty. Our estimates read B(B s → γl + l − )| q 2 ∈[1,6] GeV 2 = (6.01 ± 0.08 ± 0.70)10 −9 B(B d → γl + l − )| q 2 ∈[1,6] GeV 2 = (1.02 ± 0.15 ± 0.05)10 −11 . (8.7) The first uncertainty in these estimates reflects merely the 10% uncertainty in the B → γ transition form factors. The second uncertainty reflects the uncertainty in the contributions of the light vector mesons ρ, ω, and φ. We would like to emphasize that in the case of the B s → γl + l − transition, the dominant contribution is given by the narrow φ-meson pole. In the case of the B d → γl + l − transition, the known contribution of the vector resonances is less important, and as the result, the branching ratio uncertainty reflects to a large extent the form factor uncertainty of 10%. Table 5 gives the branching ratios integrated over several bins below the charmonia region (the charmonia region 0.33 ≤ q 2 /M 2 Bs ≤ 0.55 is normally excluded from the data analysis). Table 6 shows the results of the integration  above the charmonia region, from q 2 = 0.55M 2 Bs to q 2 = M 2 B − 2M B E γ for different values of E γ , the photon energy in the B-meson rest frame: since the Bremsstrahlung contribution leads to a divergence of the branching ratio at large q 2 , a certain cut on the photon energy is required. The values in the range 0.1 GeV < E γ min < 0.5 GeV correspond to the photon selection criteria at the Belle II detector [69], while the interval 0.5 GeV < E γ min < 1 GeV is relevant for those at the LHCb detector [70,71].
In the actual data analysis, the photon-energy cut is applied in the laboratory frame, not in the B-rest frame. However, as can be seen from Tables 5 and 6, the specific value of the energy have a marginal impact on the total branching ratio, since the contribution of the region above charmonia resonances, after applying the cuts, contributes to the branching ratio at the level below 10%.
2. Ratio R µ/e of the muon/electron distributions in B → γl + l − decays Fig. 14 shows the ratio R µ/e of the differential distributionsB → γµ + µ − toB → γe + e − . Such a ratio is a standard observable for probing violations of the lepton universality in rare FCNC semileptonic decays B → (P, V )l + l − ; recently, this ratio was proposed as a useful variable also in radiative leptonic decays [30]. The ratio is found to be close to unity in the region q 2 ≤ 5 GeV 2 . Notice however, that in B → γl + l − decays, unlike the B → (P, V )l + l − processes, the lepton-mass effects, mainly the Bremsstrahlung contribution to the amplitude, come into the game. At large q 2 , the terms proportional to C 7γ can be neglected compared to those containing C 9V and C 10A and one obtains the following expressions for different contribition to the decay rate: For the electrons in the final state, all terms containing m l may be neglected; for the muons Γ 12 gives a negligible contribution compared to Γ 1 and Γ 2 . Finally, for q 2 above the narrow charmonia region, one finds an approximate relation Obviously, the Bremsstrahlung contribution drives the ratio R µ/e far above unity at q 2 → M 2 B . In practice, at q 2 ≥ 15 GeV 2 the ratio R µ/e is sensitive to the details of the q 2 -behaviour of the transition form factors F A,V . So, if the lepton universality has been verified at small q 2 , measuring R µ/e at large q 2 provides a direct access to the q 2 -dependence of the B → γ transition form factors.

Forward-backward asymmetry
The differential forward-backward asymmetry is given by the relation whereŝ = q 2 /M 2 B , θ is the angle between p and p 2 , the momentum of the negative-charge lepton. The results for A F B in theB s → γµµ and in theB d → γµµ are shown in Fig. 15 (a) and (b), respectively. The asymmetries are practically insensitive to the uncertainties in the B → γ transition form factors, as these uncertainties to large extent cancel each other in the asymmetries.
As discussed above, the relative phases of the charmonium resonances cannot be unambiguously determined on the basis of the QCD sum-rule calculation of the nonfactorizable corrections at q 2 ≤ 4m 2 c . The results shown in Fig. 15 (a) and (b), are obtained making use of the QCD sum-rule results for nonfactorizable contributions at small q 2 and employing the conventional assignment of the signs of all charmonia contributions to be positive, following the patter of the factorizable contributions.
It should be emphasized, that A F B in the region between ψ and ψ ′ provides an unambiguous test of the relative signs of ψ and ψ ′ contributions: as displayed in Fig. 15 (c) and (d), A F B being insensitive to the patter of charmonium resonances in C 9V in the region of small q 2 , demonstrates qualitative differences depending on the relative signs between ψ and ψ ′ in the region between the resonances. Therefore, experimental study of the asymmetry in this region will allow one to check the relative signs of the ψ and ψ ′ . We point out that this observation applies not only to B → γl + l − decays, but also to B → V l + l − decays [49], where such measurement seems much more feasible.
One should take into account, however, that the measurement of the forward-backward asymmetry in the B → γl + l − decay seems to be a hard task, because the final state γl + l − does not carry any information about the flavor of the decaying B-meson. In addition, the signs of the asymmetries corresponding to B andB mesons are opposite. In the absence of flavor tagging, the total asymmetry equals zero aside from CP-violating effects. It appears that flavor tagging is impossible at LHCb; however, at Belle II one can use the fact that neutral B-mesons are produced in an entangled state. Thus, if one of the B-mesons decays to a state with a certain flavor, the other B-meson decaying to γl + l − has the opposite flavor. Now, if the interval between the decays is less than half of the oscillation period, one can claim that the flavor of the second B-meson is known with sufficient probability. For each selection procedure one can also account for the oscillations contribution and therefore improve the prediction accuracy. A method for such calculations was developed in [27].  Fig. 15: Forward-backward asymmetry forBs → γµ + µ − (left column) andB d → γµ + µ − (right column) decays. The lower plots show the asymmetries at q 2 < M 2 ψ ′ for two different relatives signs of ψ and ψ ′ contributions.

CONCLUSIONS
The paper presents a detailed analysis of nonperturbative QCD effects in rare FCNC decays B (s,d) → γl + l − decays in the Standard model. It should be emphasized that FCNC radiative leptonic decays exhibit a much more diverse structure of such effects compared to FCNC semileptonic decays.
Our main results may be summarized as follows: 1. We analyzed the B → γ * transition amplitudes induced by the the vector, axial-vector, tensor, and pseudotensor b → (s, d) quark currents. The invariant form factors, which parameterize the corresponding amplitudes, depend on two variable k ′2 and k 2 , where k ′ is the momentum emitted from the FCNC b → (s, d) vertex, and k is the momentum of the virtual photon. We worked out the constraints on these form factors imposed by electromagnetic gauge invariance. We then calculated all the form factors of interest making use of the relativistic dispersion approach based on the constituent quark picture. The appropriate subtractions in the spectral representations for the form factors have been determined from the constraints imposed by gauge invariance.
For those form factors, describing the processes with l + l − pair emitted from the FCNC vertex, the numerical results in the full q 2 -range of the B → γ decay have been obtained entirely within the dispersion approach. The resulting q 2 -dependences of the form factors exhibit the expected properties, in particular, suggest the pole beyond the physical decay region at the right location required by the quantum numbers of the appropriate meson resonances.
For those form factors, describing the l + l − pair emission from the light valence quark of the B-meson, one encounters the light vector-meson resonances in the physical decay region. To obtain the predictions for the form factors in this case, we combined the results from the direct calculation within our dispersion approach with the gauge-invariant version of the vector-meson dominance model. The numerical predictions for the form factors involving vector mesons have been updated by using the new procedure of fixing the wave-function parameters of vector mesons: namely, the wave-function parameters of the vector mesons involved have been fixed by requiring that the decay constants reproduce the latest results for f V known from other theoretical approaches.
The weak form factors of heavy mesons obtained in the dispersion approach satisfy rigorious contraints known from QCD in the heavy-quark limit. Still, the dispersion approach is a phenomenological approach representing a specific formulation of the relativistic quark model. Therefore, essentially the only way to probe its accuracy is a comparison of its predictions with the predictions from direct QCD related approaches. A comparison with the known results from lattice QCD and QCD sum rules, where such results are available, suggest that the uncertainty of the form factors from the dispersion approach does not exceed the level of about 10%. We therefore assign this accuracy to our predictions for the B → γ transition form factors.
2. We performed a detailed study of the charm-quark contributions and worked out the general constraints on these contributions imposed by electromagnetic gauge invariance. Assuming the similarity of the charm-loop contributions to B → γl + l − and to B → V l + l − amplitudes, V the vector meson, we obtained numerical predictions for the charm contributions to the B → γl + l − amplitude making use of the existing estimates of the nonfactorizable effects in the B → V l + l − decays. We demonstrated that the results for the nonfactorizable corrections available at q 2 below the charm threshold do not allow one to resolve the possible ambiguity in the relative charmonium resonance phases. Additional inputs are necessary to determine the phases unambiguously. We have shown that the forward-backward asymmetry in the q 2 -range between ψ and ψ ′ provides an efficient probe of the relative charmonium contributions.
The first error in these predictions reflects the uncertainty of the B → γ transition form factors; the second error reflects the uncertainty in the contributions of the light vector resonances (ρ, ω, φ).