2023 Charming-loop contribution to B s → γγ decay

,

A number of theoretical analyses of nonfactorizable (NF) charming loops in FCNC B-decays has been published: In [12], an effective gluon-photon local operator describing the charm-quark loop has been calculated as an expansion in inverse charm-quark mass m c and applied to inclusive B → X s γ decays (see also [13,14]); in [15], NF corrections in B → K * γ using local operator product expansion (OPE) have been studied; NF corrections induced by the local photon-gluon operator have been calculated in [16,17] in terms of the light-cone (LC) 3-particle antiquarkquark-gluon Bethe-Salpeter amplitude (3BS) of K * -meson [18][19][20] with two field operators having equal coordinates, 0|s(0)G µν (0)u(x)|K * (p) , x 2 = 0. Local OPE for the charm-quark loop in FCNC B decays leads to a power series in Λ QCD m b /m 2 c ≃ 1.To sum up O(Λ QCD m b /m 2 c ) n corrections, Ref. [21] obtained a nonlocal photon-gluon operator describing the charm-quark loop and evaluated its effect making use of 3BS of the B-meson in a collinear LC configuration 0|s(x)G µν (ux)b(0)| Bs (p) , x 2 = 0 [22,23].The same collinear approximation [known to provide the dominant 3BS contribution to meson tree-level form factors [24,25]] was applied also to the analysis of other FCNC B-decays [26].
In later publications [27][28][29][30], it was demonstrated that the dominant contribution to FCNC B-decay amplitudes is actually given by the convolution of a hard kernel with the 3BS in a different configuration -a double-collinear light-cone configuration 0|s(y)G µν (x)b(0)| Bs (p) , y 2 = 0, x 2 = 0, but x y = 0.The corresponding factorization formula was derived in [30].The first application of a double-collinear 3BS to FCNC B-decays was presented in [31].
In this paper, we study NF charming loops in B s → γγ decays making use of the generic 3BS of the B-meson.The main new features of this paper compared to the previous analyses, in particular to [31], are as follows: (i) The generic 3BS of the B-meson contains new Lorentz structures (compared to the collinear and the double-collinear approximations) and new three-particle distribution amplitudes (3DAs) that appear as the coefficients multiplying these Lorentz structures.Analyticity and continuity of the 3BS as the function of its arguments at the point xp = yp = x 2 = y 2 = 0 leads to certain constraints on the 3DAs [30] which we take into account.(ii) We derive the convolution formulas for the B s → γ * γ * form factors involving this generic 3BS, and obtain the corresponding numerical predictions.We check that the deviation between our analysis and the analysis based on double-collinear 3BS differ by O(λ Bs /M B ) terms that in practical calculations give a ∼20% difference.The paper is organized as follows: Section 2 presents general formulas for the top and charm contribution to the B s → γγ amplitude.Section 3 considers the AV V charm-quark triangle and gives a convenient representation for this quantity via the gluon field strength G µν merely (not involving A µ itself).In Section 4, properties of the 3BS of the B-meson in the general noncollinear kinematics are discussed and properly modified 3DAs are constructed.Section 5 presents the numerical results for the form factors and for the nonfactorizable charm-loop correction to the B s → γγ amplitude.Section 6 gives our concluding remarks.Appendix A compares the definition of the amplitude adopted in this paper with the one of [32].Appendix B gives details of the numerical results for the form factors.A standard theoretical framework for treating 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 [33][34][35]: G F is the Fermi constant and V ij are CKM matrix elements.The SM Wilson coefficients relevant for our analysis at the scale µ 0 = 5 GeV have the following values [corresponding to C 2 (M W ) = −1]: C 1 (µ 0 ) = 0.241, C 2 (µ 0 ) = −1.1,C 7 (µ 0 ) = 0.312 [21,[34][35][36][37].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 the case of our interest virtual c quarks -should be calculated and added to the diagrams generated by the effective Hamiltonian.

B. The penguin contribution
The top-quark contribution to B s → γγ decay is generated by penguin operator in (2.1)1 2) The sign of the b → dγ effective Hamiltonian (2.2) 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 corresponding to the definition of the covariant derivative in the form The amplitude of the B → γγ transition is defined according to [32,38]: Here q, q ′ and ε, ε ′ are momenta and polarization vectors of the outgoing real photons, and F T A and F T V are the form factors F T A (q2 = 0, q ′2 = 0) and F T V (q 2 = 0, q ′2 = 0).The latter are defined as [32,38,39]: and satisfy a rigorous constraint F T A (q 2 , 0) = F T V (q 2 , 0).Notice that the strange-quark charge Q s (or Q b in the 1/m b -subleading diagram where the photon is emitted by the b-quark) is included in the form factors F T A and F T V [32].

C. Nonfactorizable charm-quark loop correction to Bs → γγ
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: where differing from each other in the way color indices i, j, k, l are contracted.By Fierz transformation O 2 may be written in the following form (for anticommuting spinor fields): The singlet-singlet operator O 1 and the singlet-singlet part of O 2 at the leading order generate factorizable charm contributions to the B → γγ amplitude.These factorizable contributions vanish for real photons in the final state.A nonzero contribution is induced by the octet-octet part of the operator O 2 and needs the emission of one soft gluon from the charm-quark loop.So relevant for us is the octet-octet operator Therefore, similar to the top contribution, we find (2.12) Here quark fields are understood as Heisenberg field operators with respect to the SM interactions.Expanding them to the second order in electromagnetic interaction and to the first order in strong interaction gives where quark electormagnetic current has the form j e.m. α = e Q i qi γ α q i .Eq. (2.13) may be rewritten as: Other diagrams are those corresponding to an opposite direction in the charm-quark loop and diagrams with the interchanged photons q ↔ q ′ , ε ↔ ε ′ .
A detailed treatment of the operator, describing charm-quark triangle [second line in Eq. (2.14)] is given in the next Section.Here we only notice two important features of this operator: • The cγ µ c part of the V − A weak current does not contribute and one is left with V V A charm-quark triangle.
• The V V A charm-quark triangle contracted with the gluon field B b ν may be written as a gauge-invariant nonlocal operator containing gluon field strength G b να for any gluon momentum (cf.[14]).Making use of the result for the charm-quark V V A triangle from the next Section, we obtain the following expression for the amplitude: A ρη (q, q ′ ) = 1 (2π) 8  dkdye −i(k−q ′ )y dxdκe −iκx Γ µνρ(ab) cc (κ, q) 0|s(y) (2.16) For real photons in the final state, the amplitude A ρη (q, q ′ ) has the same Lorentz structure as the penguin amplitude and contains two form factors H V = H V (q 2 = 0, q ′2 = 0) and H A = H A (q 2 = 0, q ′2 = 0) (Appendix A presents comparison with the form factors defined in [32]): Comparing Eqs.(2.4) and (2.17), and taking into account that V tb V * ts ≃ −V cb V * cs , it is convenient to desctribe the effect of charm as an additions to the Wilson coefficient C 7γ (however, non-universal, i.e. different in the axial and vector Lorentz structures): with Our goal will be to calculate these corrections.The B-meson structure contributes to the A B→γγ charm amplitude via the full set of 3BS 0|s(y with Γ i the appropriate combinations of γ-matrices.This quantity is not gauge invariant, since it contains field operators at different locations.To make it gauge-invariant, one needs to insert Wilson lines between the field operators.To simplify the full consideration, it is convenient to work in a fixed-point gauge, where the Wilson lines reduce to unity factors.As first noticed in [27], the dominant contribution of charm to amplitudes of FCNC B decays comes from the "double collinear" LC configuration [30], where x 2 = 0, y 2 = 0, but xy = 0, i.e. 4-vectors x and y are not collinear.Respectively, we need to parametrize the 3BS in this kinematics; this is discussed in Sect. 4. But before studying 3BS, we present in the next Section a convenient representation for the operator describing the contribution of charm-quark loop.

CHARM-QUARK V V A TRIANGLE
The charm-quark loop contribution is described by the three-point function (see Fig. 2): where q is the momentum of the external virtual photon (vertex containing index ρ) and κ is the gluon momentum (vertex containing index ν).Here t c , c = 1, . . ., 8 are SU c (3) generators normalized as Tr(t a t b ) = 1 2 δ ab .The octet current c(0)γ µ (1 − γ 5 )t a c(0) is a charm-quark part of the octet-octet weak Hamiltonian.Its vector piece does not contribute to Γ µνρ (ab) cc (Furry theorem) and will be omitted.Taking into account vector-current conservation, it is convenient to parametrize Γ µνρ cc (κ, q) as follows [40] Γ The form factors F 0,1,2 are functions of three independent invariant variables q 2 , κ 2 , and κq.The lowest order QCD diagrams describing Γ µνρ cc (κ, q) are shown in Fig. 2. A convenient representation of the form factors has the form [41] This representation may be applied to the physical amplitude in the region of the external momenta far below the thresholds, q 2 , κ 2 , (κ + q) 2 ≪ 4m 2 c .Taking into account the momentum distribution of quarks and gluons inside the B-meson, the dominant contribution of the charm-quark loop to the B-decay amplitude comes from the region κ 2 ∼ Λ 2 QCD , (q + κ) 2 < 0. So, the representation (3.3) is applicable and proves convenient for numerical calculations.As the next step, one takes the convolution of the amplitude (3.2) with the gluon field B ν (x).The Lorentz structures multiplying F 0 and F 1 contain ǫ νκα1α2 with some indices α 1 and α 2 .After multiplying by B ν (x) and performing parts integration, their contribution may be reduced to the convolution with the gluon strength tensor G αν : The Lorentz structure multiplying F 2 at first glance does not have this property.However, using the identity multiplying it by κ α1 κ α2 q α6 , and setting α 3 → µ, α 4 → ν, α 5 → ρ, this Lorentz structure takes the form and may be also reduced to the convolution with G αν .Finally, the operator describing the contribution of the charm-quark loop takes the form

3BS OF THE B-MESON IN A NONCOLLINEAR KINEMATICS
As already mentioned, the contribution of collinear LC configuration dominates the 3BS corrections to the B → π, K form factors.These corrections reflect the following picture: in the rest frame of the B-meson, a fast light quark, produced in weak decay of an almost resting b-quark, emits a soft gluon and continues to move practically in the same direction, before it fragments into the final light meson.
Contributions of charming loops in FCNC B-decay have a qualitatively different picture [29]: In the rest frame of the decaying B-meson, two fast systems produced in the weak decay of an almost resting b-quark move in opposite space directions.Formulated in terms of the LC variable, this means that the s-quark produced in weak decay moves along one of the LC directions, whereas the cc-pair moves along the other LC direction.Introducing vectors n µ and , one finds that the dominant contribution of charming loops to an FCNC B-decay amplitude comes from the double-collinear configuration [27][28][29][30][31] when the coordinates of the field operators in 0|s(y)G µν (x)b(0)| Bs (p) are alligned along the orthogonal light-cone directions x µ ∼ n µ , y µ ∼ n ′ µ .The 3BS amplitude in the collinear and the double-collinear kinematics contain the same Lorentz structures [42] but the distribution amplitudes corresponding to the collinear and the double-collinear kinematics differ from each other.
In this paper we do not consider the double-collinear approximation but make use of the general noncollinear 3BS.This quantity contains new Lorentz structures and new 3DAs.The B s → γγ amplitude calculated using the general noncollinear 3BS differs by terms O(λ Bs /M B ) from the amplitude calculated within the double-collinear approximation.

A. Collinear 3BS of B-meson
We summarize in this Section well-known results concerning the collinear 3BS that will be used for constructing a generalization to a noncollinear kinematics appropriate for charming loops in FCNC B-decays.

The Lorentz structure of the collinear 3BS
We start with the collinear LC 3BS [23], where the arguments of the s-quark field, s(y), and the gluon field G να (x) are collinear to each other, x = uy, u = 0 is a number (in this case x 2 = 0 leads to y 2 = 0): where D(ω, λ) takes into account rigorous constraints on the variables ω and λ: In Eq. (4.1), Γ is an arbitrary combination of Dirac matrices, v µ = p µ /M B , and all 8 DAs (Ψ A , Ψ V , etc) are functions of two dimensionless arguments 0 < λ < 1 and 0 < ω < 1.Here λ refers to the momentum carried by the s-quark, and ω refers to the momentum carried by the gluon.The normalization conditions for Ψ A and Ψ V have the form [23]: Some of the Lorentz structures in (4.1) contain factors x µ /xp or x µ x ν /(xp) 2 .Since 3BS (4.1) is a continuous regular function at x 2 = 0 and xp = 0, the absence of singularities at xp → 0 leads to the following constraints: These constraints are obtained by expanding the exponential in the integral representation (4.1) under the condition x µ = u y µ to the necessary order and requiring that the coefficients multiplying terms singular in xp → 0 vanish.

Twist expansion of the 3DAs
The DAs in (4.1) have no definite twist.According to [23], the distribution amplitudes might be written as an expansion in functions with definite twist as follows: where we keep the contributions up to twist 6 inclusively (the subscript "i" in φ i and ψ i denote the twist value).

Model for DAs entering the collinear 3BS
The powers of ω and λ determine the behaviour at small quark and gluon momenta.This power scaling is related to the conformal spins of the fields and remains the key property of the model.
The starting point of our analysis will be the set of DAs in LD model of [23] for twist 3-and 4, complemented by twist 5 and 6 DAs reconstructed using the constraints (4.4) [43]: ) ) ) ) Dimensionless parameter ω 0 is related to λ B , the inverse moment of the B-meson LC distribution amplitude, as For this model, the integration limits take the following form (2ω 0 < 1): However, as we shall show shortly, certain modifications of the integrated DAs [which emerge when performing parts integrations of the 3BS (4.1)] at large values of ω and λ will be necessary in order to satisfy the continuity of the 3BS considered in a noncollinear kinematics.

B. Generalization to a noncollinear kinematics
When the coordinates x and y are independent variables, the 3BS has the following decomposition that involves more Lorentz structures and more 3DAs compared to the collinear approximation: All invariant amplitudes Φ = Ψ A , Ψ V , . . .are functions of 5 variables, Φ(ω, λ, x2 , y 2 , xy), for which we may write Taylor expansion in x 2 , y 2 , xy.Here we limit our analysis to zero-order terms in this expansion.The corresponding zero-order terms in Φ's are functions of dimensionless arguments λ and ω and are referred to as the DAs.These DAs contain at least the kinematical constraint θ(1 − ω − λ).However, the DAs may have support in more restricted areas: e.g., the DAs of the LD model Eqs. (4.17) The amplitude (4.16) contains two independent kinematical singularities 1/xp and 1/yp in the Lorentz structures.These kinematical singularities of the Lorentz structures should not be the singularities of the amplitude; this requirement leads to certain constraints which we are going to consider now.We shall present these constraints for the case when all DAs contain θ(2ω 0 − ω − λ), 2ω 0 < 1.
For the amplitudes of the type F , the Lorentz structures of which contain first power of 1/xp or 1/yp, the appropriate constraint is obtained by expanding the exponential in (4.16) to zero order and requiring that the singular terms vanish: Let us introduce the primitives which, by virtue of (4.18), vanish at the boundary of the DA support region: For the functions Z and W , the Lorentz structure of which contain 1/(xp) 2 and 1/(yp) 2 , the exponential should be expanded to first order leading to: By introducing primitives and double primitives (ω, λ) (ω, λ) one can show that the requirements (4.22) lead to the vanishing of these functions at the boundaries of the DAs support regions: These relations are verified making use of Eq. (4.22).For instance, Z The analogous relations are valid for W .
For calculating the contribution of (4.16) to the FCNC amplitude, the presence of the 1/xp and 1/yp structures are inconvenient.To facilitate the calculation, one may perform parts integration in ω for the structure containing 1/xp and in λ for the structure containing 1/yp.The conditions (4.21) and (4.24) lead to the absence of the surface terms when performing the parts integrations.So, we can rewrite (4.16) in the convenient form not containing the 1/xp and 1/yp factors: C. Adapting the LD model for the case of noncollinear 3BS First, let us point out that the primitives calculated for the DAs of the LD model (4.6)-(4.13)do not satisfy the constraints (4.21) and (4.24).As already emphasized, the parametrizations (4.16) and (4.26) are then not equivalent and cannot be both correct.So, there are two different possibilities to handle this situation: Since the behaviour of the DAs of the definite twist at small ω and λ are fixed, we choose the following procedure: 1.The functions Ψ A (ω, λ) and Ψ V (ω, λ) remain intact.2. We split each function Φ (Φ = X A , Y A , XA , ỸA , Z, W ) into Φ (x) (ω, λ) and Φ (y) (ω, λ) as follows:3 (4.27) 3. We make use of the DAs of LD model for the functions Φ(ω, λ) and calculate primitives of higher orders.Taking into account that higher primitives obtained in this way do not satisfy the constraints (4.21) and (4.24), we modify them in the way explained below and allow them to depend on one additional parameter a. 4. We still want that after taking the appropriate derivatives of the modified primitives, we reproduce Eq. ( 4.16) with 6 modified functions (X A , Y A , XA , ỸA , Z, W ) which in turn also depend on the additional parameter a.The simplest consistent scheme that fulfills this requirement is the following: • For the functions in Eq. (4.16) containing factors 1/xp and 1/yp (i.e.X = X A , Y A , X A , Y A ) the primitives entering Eq. (4.26) are constructed as follows: • For the functions in Eq. (4.16) containing factors 1/xp and 1/yp and 1/(xp) 2 and 1/(yp) 2 (i.e., Z and W ) the primitives entering Eq. (4.26) are constructed in a different way (the same formulas for W ): The function R(ω, λ, a) is chosen in the form such that for a = 0 the function itself and all its derivatives vanish at the boundary ω + λ = 2ω 0 .Respectively, for a = 0, all modified higher primitives of the necessary order vanish at the boundary 2ω 0 − ω − λ = 0, providing the absence of the 1/xp and 1/yp singularities at xp → 0 and yp → 0 and the continuity of the 3BS at the point x 2 = 0, y 2 = 0, xp = 0 and yp = 0.For small values of the parameter a, our model DAs reproduce well the collinear DAs including their magnitudes and power behaviour at small ω and λ but (strongly) deviate from them near the upper boundary.

Scenario II
One declares Eq. (4.26) as the starting point but calculates the necessary primitives using the 3DAs from Eq. (4.16).Then, however, Eq. (4.16) itself is not complete and should contain nonzero surface terms obtained by rewriting Eq. (4.26) to the form (4.16).These surface terms guarantee the absence of singularities in (4.16) at xp → 0 and yp → 0. This scenario does not seem to us logically satisfactory: For instance, the double-collinear limit [30,31] that may be readily taken in the 3BS given by Eq. (4.16) is not reproduced by Eq. (4.26) if the surface terms are nonzero!Nevertheless, for comparison we also present the results for the form factors calculated using this prescription referred to as Scenario II.It should be emphasized that for one and the same set of 3DAs, the form factors obtained using Scenario I in the limit a → 0 do not reproduce the form factors obtained using Scenario II: the difference amounts to certain surface terms that arise in the limit a → 0.

RESULTS FOR THE Bs → γγ AMPLITUDE AND δC7γ
We are ready to evaluate the form factors H A,V .We shall also calculate the penguin form factor F T , and in the end, the charming-loop correction to C 7γ .(κ, q) by parts integration in κ.Any factor y α may be represented as y α → ∂ ∂kα e −iky , then moving the derivative onto the strange-quark propagator by parts integration in k.Doing so, we get rid of all powers of x and y, and the integrals over x and y in (2.16) then lead to the δ functions: which allow us to further take the integrals over κ and k.In the end, the invariant functions H i , i = A, V are given by integral representations of the form where h i are linear combinations of the DAs and their primitives entering Eq. (4.26).Eq. (5.3) gives now the form factors H A,V as the convolution integrals of the known expression for the charming loops and the DAs and its primitives from (4.26).We make use of the modified LD model described above.Obviously, the form factors H A and H V depend linearly on λ 2 E and λ 2 H . So, we present the results for the form factors R iE and R iH defined as follows: (5.4) The form factors R iE,iH depend on λ Bs , but do not contain dependence on λ H and λ E .These three parameters, however, are expected to be strongly correlated and we take into account this correlation later when evaluating δC 7γ .
Table 1: Input parameters for the calculation of the form factors RiE,iH defined in Eq. (5.4).
mc(mc) ms(2 GeV) MB s M φ fB s λB s (1 GeV) a wx 1.30 GeV 0.1 GeV 5.3 GeV 1.020 GeV 0.23 GeV 0.45 ± 0.15 GeV 0.05 ± 0.05 0.5 ± 0.5 function R, Eq. (4.30), and the weight factors w x .We present our results for the form factors for w x = 0.5 and reflect the sensitivity to its variations in the ranges 0 < w x < 1 in the final uncertainties.We restrict the parameter a by the requirement that the 3DAs are affected by, say, not more than 10% in the region of "small" λ and ω, e.g., at λ/(2ω 0 ), ω/(2ω 0 ) ≤ 0.1.This leads to the restriction a ≤ 0.1.So we set the benchmark point a = 0.05 and allow its variations in the range 0 ≤ a ≤ 0.1.We shall see that the sensitivity of the calculated form factors to a > 0.1 is rather mild.
Other parameters in Table 1 are taken from [50] and [51].We are interested in obtaining the form factors R iE,iH at q 2 = 0 and q ′2 = 0. Whereas q 2 = 0 may be readily set in the integral representation (5.3), this integral representation is not expected to give reliable predictions at q ′2 = 0: The s-quark propagator at q ′2 = 0 is not sufficiently far off-shell and at q ′2 → 0 one observes a steep rise of R iE,iH related to the nearby quark singularity.This rise is unphysical as the nearest physical singularity in the form factors is located at q ′2 = M 2 φ .To avoid this problem, we choose the following strategy: we take the results of our calculation for R iH,iE (0, q ′2 ) in the interval −5 < q ′2 (GeV 2 ) < −0.7 and extrapolate them numerically to q ′2 = 0 making use of a modified pole formula which takes into account the presence of the physical pole at q ′2 = M 2 φ : . (5.5) In this way we obtain R iE,iH (0, 0). Figure 3 shows the results of our direct calculation and the fits obtained using Eq.(5.5). Figure 3(a,b) gives the form factors evaluated with the modified 3DAs for a = 0.05 (Scenario I). Figure 3(c,d) shows the results for unmodified functions (Scenario II).For comparison, Fig. 3(e,f) shows the contribution of the 3DAs ψ A and ψ V .Good news is that the latter give the dominant contribution (which does not depend on w x and a and thus reduces the uncertainty in the R i related to the values of these parameters).The contribution of other Lorentz 3DAs is however not negligible and depend on the way one handles these 3DAs: e.g. for R AE (0, 0) in Scenario I they give a negative correction of 30%.Appendix B presents details of the contributions of different Lorentz 3DAs to the form factors R iE,iH (0, q ′2 ) at q ′2 = −2 GeV 2 for Scenarios I and II.
Fig. 4 demonstrates the sensitivity of R iE,iH (0, 0) to λ Bs .Since the form factors at q ′2 = 0 are obtained via extrapolation, we take λ Bs from the range 0.3 < λ Bs (GeV) < 0.5, calculate R iE,iH (0, q ′2 ) in the range −5 < q ′2 (GeV 2 ) < −0.7 and then extrapolate to q ′2 = 0 for each value of λ Bs .The same procedure applies to the dependence on parameter a [enters the function R, Eq. (4.30)] shown in Fig. 5.We restirct a in the range 0 ≤ a ≤ 0.1, but, as seen from Fig. 5, the sensitivity to the value of a in the region a > 0.01 is rather weak anyway.
The uncertainties in our theoretical predictions for the form factors thus come from the following sources: (i) The sensitivity to the precise way one handles the 3BS at large values of ω and λ.This is probed by comparing the results obtained with our benchmark Scenario I with those from Scenario II.
(ii) The sensitivity to the precise functional form of the B-meson 3DAs; this is probed by using as our benchmark model (4.6)-(4.13)[Model IIB from [23]] and comparing with a different model [Model IIA given by Eq. (5.23) from [23]].The uncertaintes related to (i) and (ii) are denotes as [3DA].(iii) The sensitivity to the numerical values of the parameters of the 3DAs, mainly parameter a of the function R, Eq. (4.30).Notice that we do not perform any averaging over λ Bs ; we keep the full dependence on this parameter in the form factors and in our final result for the correction δC 7γ .
(iv) The sensitivity to the extrapolation procedure from the interval where the form factors may be calculated by our approach to the physical point q ′2 = 0. We make use of the calculations of the form factors in the interval −5 ≤ q ′2 (GeV 2 ) ≤ −0.7.However, the value at q ′2 = 0 obtained by the extrapolation is sensitive to the precise choice of the upper boundary of this interval.For instance, moving the upper boundary in the range (−1 ÷ −0.5) GeV 2 leads to the variations of the extrapolated value of the form factor at q ′2 = 0 by ±5%.We therefore assign the additional uncertainty of 5% (denoted as [extr]) related to the extrapolation procedure.
Table 2 compares the form factors at q 2 = q ′2 = 0 for two different sets of 3DAs and for two different Scenarios to treat the 3DAs.The individual form factors R iE and R iH demonstrate a sizeabe dependence on the Scenario.However, in Table 2: Our results for the form factors RiE,iH in GeV −1 defined in Eq. (5.4).In addition to the results obtained in Scenario I and Scenario II for the basic set of 3DAs given in (4.6)-(4.13),referred to as (B), we present the results obtained using Scenario I for 3DAs of Model IIA of Eq. (5.23) from [23], referred to as (A). the combinations appropriate for the calculation of δC 7γ , R V E + 2R V H and R AE + 2R AH these uncertainties cancel to great extent and we obtain for the benchmark value λ Bs = 0.45 GeV: and we have slightly increased our final "theoretical uncertainty".In the end, for a given value of the λ Bs , we predict the combination of the form factors (5.6) and (5.7) relevant for the calculation of C 7γ (see Subsection C) with about 8% accuracy.
B. The penguin form factor FT (0, 0) The form factors F T A,T V (q 2 , q ′2 ) may be calculated via the B-meson 2DAs using HQET formula (see e.g.[25]) leading to the following expression for the F T (0, 0) = F T A (0, 0) = F T V (0, 0): where (5.10) The absence of the kinematical singularity at vx → 0 in Eq. (5.8) requires that the primitive Φ(λ) vanishes at the boundaries of the 2DA support region.Then the parts integration in λ, necessary to handle the 1/vx term, does not contain any nonzero surface term.All contributions O(m s /M B ) in the numerator of (5.9) are omitted; keeping such terms will be inconsistent as contributions of this order arise also from the diagrams describing photon emission from the B-quark which are not taken into account.According to the estimates obtained in [32], corrections due to the photon emission from the b-quark amount to 10-20% of the leading contribution (5.9).
For the 2DAs we use the same model from [23] as we do for 3DAs: (5.12) An explicit check shows that Φ(λ = 0) = Φ(2ω 0 ) = 0.The contribution of the term ∼ (λ 2 E − λ 2 H ) to the form factor turns out numerically negligible and may be safely omitted.The form factor FT A(0, q ′2 ), where q ′ is the momentum of the photon emitted from the s-quark, for 2DAs given by (5.11) and (5.12) and our benchmark value λB s = 0.45 GeV.Dotted line -calculation results; solid line -interpolation using the modified pole formula Eq. (5.5).(b) The dependence of FT A(0, 0) on λB s in the range λB s = 0.45 ± 0.15 GeV.
Let us obtain numerical estimates for F T (0, 0).Notice that the numerator of the integrand in Eq. (5.9) contains factor λ, so no singularity at λ → 0 arises in the integrand even in the limit m s → 0. Therefore, the F T A (0, 0) may be calcualted by applying directly the representation for the form factor Eq. (5.9).[Recall that in order to obtain the form factors H A,V (0, 0) we had to perform extrapolation from the region q ′2 −0.5 GeV 2 .] Figure 6(a) shows that the q ′2 -behaviour of F T A (0, q ′2 ) is well compatible with the location of the physical pole at q ′2 = M 2 φ in a broad range of q ′2 < 0, up to q ′2 = 0. Figure 6(b) illustrates the sensitivity of F T (0, 0) to λ Bs .This dependence is not negligible and partly compensates the λ Bs -dependence of H A,V (0, 0), leading to more stable predictions for δC 7γ .
For our further numerical estimates we use F T (0, 0) = 0.155 ± 0.015 for λ Bs = 0.45 GeV, where the uncertainty of ∼10% is assigned on the basis of the size of the neglected 1/M B -suppressed contributions which have been calculated in [32].

C. δC7γ
We are ready to evaluate the relative contribution of the nonfactorizable charming loops given by Eq. (2.19): This expression is useful as it allows us to separate the uncertainties related to the precise values of the Wilson coefficients and m b to be used in the numerical estimates and those related to the description of the B-meson structure.
Let us focus on ρ (i) cc .The parameters λ Bs , λ E and λ H , entering ρ cc , are strongly correlated with each other and should not be treated as independent quantities.As noticed in [23] (see also [52,53]), QCD sum rules suggest an approximate relation λ 2 H = 2λ 2 E .Combining this relation with the constraints from the QCD equations of motion [23], one obtains approximate relations  Fig. 7 shows the dependence of ρ i cc on λ Bs for R iH,iE given in Fig. 4 and F T given in Fig. 6.The results presented in the plot correspond to Scenario I calculated with 3DAs of Model IIB.To a good accuracy, we may set ρ (5.15) The final result for δC 7γ is very sensitive to the precise value of λ Bs .Recall, however, that for a given value of λ Bs , ρ cc may be calculated with an accuracy around 10%.So we prefer to present our results for δC 7γ as the function of λ Bs .For our benchmark point λ 0 Bs = 0.45 GeV, we find δC 7γ (λ 0 Bs ) = 0.045 ± 0.004. (5.16) For λ Bs in the range 0.3 < λ Bs (GeV) < 0.6, the corresponding δC 7γ covers the range δC 7γ = (2 ÷ 10)%. (5.17) We therefore conclude that the effect of nonfactorizable charming loops is expected at the level of a few percent.As soon as the parameter λ Bs is known with a better accuracy, our results allow one to obtain a (relatively) accurate estimate for δC 7γ .

DISCUSSION AND CONCLUSIONS
We presented a detailed analysis of NF charming loops in FCNC B s → γγ decay and reported the following results: (i) We derived and made use of the expression for the V V A charm-quark loop that is fully given in terms of the gluon field strength G µν (x).This has an advantage that no explicit use of any specific gauge for the gluon field is necessary.Still, nonlocal operator describing the charm-loop contribution to the amplitude of FCNC B-decay contains field operators at different coordinates, s(y)G µν (x)b(0), and thus needs the Wilson lines joining the field operators at different points.These Wilson lines are reduced to unity operators in the Fock-Schwinger gauge x µ A µ (x) = 0, so this gauge is used implicitly.
(ii) We studied the generic non-collinear 3BS of the B-meson; this quantity contains new Lorentz structures and new 3DAs compared to collinear and double-collinear 3BS.We took into account constraints on the non-collinear 3BS coming from the requirement of analyticity and continuity [30] and implemented proper modifications of the corresponding 3DAs X i (ω, λ) at large values of their arguments.
(iii) We calculated the form factors H i (q 2 = 0, q ′2 = 0), i = A, V , describing the B → γγ amplitude.Whereas q 2 = 0 [q the momentum emitted from the charm-quark loop] may be set directly in the analytic formulas, the physical point q ′2 = 0 [q ′ the momentum emitted by the light s-quark] was reached by an extrapolation from the spacelike region −5 ≤ q ′2 (GeV 2 ) ≤ −(0.5 ÷ 1).An explicit dependence of the form factors and the correction δC 7γ on the parameter λ Bs was calculated.Taking into account all uncertainties (excluding that of λ Bs ), we found that for any specific value of λ Bs , δC 7γ may be obtained with better than 8% accuracy.For our benchmark point λ 0 Bs = 0.45 GeV, we found δC 7γ (λ 0 Bs ) = 0.045 ± 0.004.
For λ Bs in the range 0.3 < λ Bs (GeV) < 0.6, δC 7γ covers the range 2-10%.Thus, one should expect the NF charm-loop correction to B → γγ decays at the level of a few %.As soon as a more accurate value of λ Bs is available, our results allow one to obtain the corresponding δC 7γ .
(iv) In the double-collinear kinematics that dominates the NF charm-loop contribution to FCNC B-decay amplitudes in the HQ limit [30], the leading contribution to the amplitude is given by the convolution of the form factor F 0 describing the charm-quark loop Eq. (3.2) and the following combination of the 3DAs [31,42]: [The explicit expressions for these 3DAs show that this combination is proportional to (λ 2 E + λ 2 H ).] Eq. (6.1) implies that in the HQ limit, R V E = R V H and R AE = R AH .Our results presented in Fig. 3 show that these relations are sizeably violated by O(λ Bs /M B ) corrections that come into the game via other Lorentz structures describing the charm-quark loop (3.2) and other 3DAs and worth to be taken into account.Numerically, these corrections are around 20%.
It might be useful to notice that the B s → γγ decay amplitude receives contributions from the weak-annihilation type diagrams [54][55][56].The weak-annihilation mechanism differs very much from the mechanism discussed in this paper and is therefore beyond the scope of our interest here.However, weak-annihilation diagrams should be taken into account in a complete analysis of B s → γγ decays.
In conclusion, we emphasize that the approach of this paper may be readily applied to the analysis of nonfactorizable charming loops in other FCNC B-decays and looks promising for treating NF effects in nonleptonic B-decays (see e.g.[57]).
Appendix B: Anathomy of the form factors HA,V (0, q ′2 0 ) at q ′2 0 = −2 GeV 2 .This Appendix presents detailed numerical results for the form factors at q 2 = 0 and q ′2 = q ′2 0 = −2 GeV 2 , giving separately contributions of different Lorentz 3DAs for the coefficients R iE and R iH defined in Eq. (5.4).Table 3 gives the results obtained with Scenario I of Section 4 C for 3DAs given by Eqs.(4.6)-(4.13)[Model IIB of Eq. (5.27) from [23]].Table 4 presents the results obtained with Scenario II and the same set of 3DAs.Table 5 gives the results corresponding to Scenario I but using different 3DAs, those of Model IIA of Eq. (5.23) from [23].Separate contributions coming from different Lorentz 3DAs are isolated keeping the explicit dependence on w x and w y , w x + w y = 1.The results correspond to λ Bs = 0.45 GeV and the central values of m c and m s .For 3DAs modified according to Scenario I, a = 0.05.

2 .
TOP AND CHARM CONTRIBUTIONS TO Bs → γγ A. The b → d, s effective Hamiltonian

. 14 )Fig. 1 Fig. 1 :
Fig.1shows one of the corresponding diagrams when the photon is emitted by the B-meson valence s-quark.We will neglect the 1/m b -suppressed contribution when the photon is emitted by the valence b-quark.

. 26 )
We emphasize that if the conditions (4.18) and(4.22)are not fulfilled, then the parametrizations (4.16) and (4.26) are not equivalent: they differ by a nonzero surface term.Noteworthy, the requirements of the absence of singularities at xp → 0 and yp → 0 and the continuity of the 3BS (4.16) at x 2 = 0, y 2 = 0, xp = 0 and yp = 0 lead to a number of constraints on the DAs which guarantee the equivalence of the forms (4.16) and (4.26).

A.
The charming-loop form factors HA,V (0, 0) Using Eqs.(2.15) and (2.16) and calculating the trace in (4.26) for Γ = γ η (/ k + m s )γ µ (1 − γ 5 ), (5.1) leads to the B s → γγ amplitude expressed via the DAs and their primitives.The expression (4.26) contains powers of x and/or y, but these are easy to handle.Let us go back to Eq. (2.16): Any factor x α may be represented as x α → ∂ ∂κα e −iκx then taking the κ-derivative over to Γ µνρα cc

Fig. 3 :
Fig. 3: The contributions RiE and RiH to the form factors Hi [i = A, V ] defined according to Eq. (5.4) for different scenarios of treating the B-meson 3DA.(a,b): The appropriate modifications of the 3DAs XA, YA, etc. are taken into account (Scenario I).(c,d): The form factors are calculated using Eq.(4.26) with the primitives obtained from XA, YA, etc. without modifications (Scenario II).(e,f): Only the contributions of ΨA and ΨV are taken into account.Dashed lines show the calculation results for λB s = 0.45 GeV, a = 0.05, wx = wy = 0.5 and other inputs from Table1, and solid lines are the fits using Eq.(5.5).

Fig. 6 :
Fig.6: (a) The form factor FT A(0, q ′2 ), where q ′ is the momentum of the photon emitted from the s-quark, for 2DAs given by (5.11) and (5.12) and our benchmark value λB s = 0.45 GeV.Dotted line -calculation results; solid line -interpolation using the modified pole formula Eq. (5.5).(b) The dependence of FT A(0, 0) on λB s in the range λB s = 0.45 ± 0.15 GeV.

10 5
RV E 10 5 RV H 10 5 (RV E + 2RV H ) 10 5 RAE 10 5 RAH 10 5 (RAE + 2RAH ) As the result, ρ i cc turns out to be the function of one variable, λ Bs , and the appropriate combinations of the form factors that determines the correction δC 7γ are R V E + 2R V H and R AE + 2R AH , presented in Table2.