# Charming-loop contribution to $B_{s}\to\gamma\gamma$ decay

Ilia Belova, Alexander Berezhnoyb, and Dmitri Melikhovb,c,d aINFN, Sezione di
Genova, via Dodecaneso 33, I-16146 Genova, Italy
bD. V. Skobeltsyn Institute of Nuclear Physics, M. V. Lomonosov Moscow State
University, 119991, Moscow, Russia
cJoint Institute for Nuclear Research, 141980 Dubna, Russia
dFaculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090 Vienna,
Austria

###### Abstract

We present a detailed theoretical study of nonfactorizable contributions of
the charm-quark loop to the amplitude of the $B_{s}\to\gamma\,\gamma$ decay.
This contribution involves the $B$-meson three-particle Bethe-Salpeter
amplitude, $\langle 0|\bar{s}(y)G_{\mu\nu}(x)b(0)|\bar{B}_{s}(p)\rangle$, for
which we take into account constraints from analyticity and continuity. The
charming-loop contribution of interest may be described as a correction to the
Wilson coefficient $C_{7\gamma}$, $C_{7\gamma}\to C_{7\gamma}(1+\delta
C_{7\gamma})$. We calculate an explicit dependence of $\delta C_{7\gamma}$ on
the parameter $\lambda_{B_{s}}$. Taking into account all theoretical
uncertainties, $\delta C_{7\gamma}$ may be predicted with better than 10%
accuracy for any given value of $\lambda_{B_{s}}$. For our benchmark point
$\lambda_{B_{s}}=0.45$ GeV, we obtain $\delta C_{7\gamma}=0.045\pm 0.004$.
Presently, $\lambda_{B_{s}}$ is not known with high accuracy, but its value is
expected to lie in the range $0.3\leq\lambda_{B_{s}}({\rm GeV})\leq 0.6$. The
corresponding range of $\delta C_{7\gamma}$ is found to be $0.02\leq\delta
C_{7\gamma}\leq 0.1$. One therefore expects the correction given by charming
loops at the level of at least a few percent.

##  1 Introduction

Charming loops in rare flavour-changing neutral current (FCNC) decays of the
$B$-meson have impact on the $B$-decay observables neubert  and provides an
unpleasant noise for the studies of possible new physics effects (see, e.g.,
ciuchini2022 ; ciuchini2020 ; ciuchini2021 ; diego2021 ; matias2022 ;
gubernari2022 ; hurth2022 ; stangl2022 ; diego2023a ; diego2023b ).

A number of theoretical analyses of nonfactorizable (NF) charming loops in
FCNC $B$-decays has been published: In voloshin , 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\to
X_{s}\gamma$ decays (see also ligeti ; buchalla ); in khod1997 , NF
corrections in $B\to K^{*}\gamma$ using local operator product expansion (OPE)
have been studied; NF corrections induced by the local photon-gluon operator
have been calculated in zwicky1 ; zwicky2  in terms of the light-cone (LC)
3-particle antiquark-quark-gluon Bethe-Salpeter amplitude (3BS) of
$K^{*}$-meson braun ; ball1 ; ball2  with two field operators having equal
coordinates, $\langle 0|\bar{s}(0)G_{\mu\nu}(0)u(x)|K^{*}(p)\rangle$,
$x^{2}=0$. Local OPE for the charm-quark loop in FCNC $B$ decays leads to a
power series in $\Lambda_{\rm QCD}m_{b}/m_{c}^{2}\simeq 1$. To sum up
$O(\Lambda_{\rm QCD}m_{b}/m_{c}^{2})^{n}$ corrections, Ref. hidr  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
$\langle 0|\bar{s}(x)G_{\mu\nu}(ux)b(0)|\bar{B}_{s}(p)\rangle$, $x^{2}=0$
japan ; braun2017 . The same collinear approximation [known to provide the
dominant 3BS contribution to meson tree-level form factors braun1994 ;
offen2007 ] was applied also to the analysis of other FCNC $B$-decays
gubernari2020 .

In later publications mk2018 ; m2019 ; m2022 ; m2023 , 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 $\langle
0|\bar{s}(y)G_{\mu\nu}(x)b(0)|\bar{B}_{s}(p)\rangle$, $y^{2}=0$, $x^{2}=0$,
but $x\,y\neq 0$. The corresponding factorization formula was derived in m2023
. The first application of a double-collinear 3BS to FCNC $B$-decays was
presented in wang2022 .

In this paper, we study NF charming loops in $B_{s}\to\gamma\gamma$ 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 wang2022 , 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 m2023  which we take into account.

(ii) We derive the convolution formulas for the $B_{s}\to\gamma^{*}\gamma^{*}$
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(\lambda_{B_{s}}/M_{B})$ terms that in practical calculations give a
$\sim$20% difference.

The paper is organized as follows: Section 2 presents general formulas for the
top and charm contribution to the $B_{s}\to\gamma\gamma$ amplitude. Section 3
considers the $\langle AVV\rangle$ charm-quark triangle and gives a convenient
representation for this quantity via the gluon field strength $G_{\mu\nu}$
merely (not involving $A_{\mu}$ 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}\to\gamma\gamma$ amplitude. Section 6 gives our concluding
remarks. Appendix A compares the definition of the amplitude adopted in this
paper with the one of mnk2018 . Appendix B gives details of the numerical
results for the form factors.

##  2 Top and charm contributions to $B_{s}\to\gamma\gamma$

###  2.1 The $b\to d,s$ effective Hamiltonian

A standard theoretical framework for treating FCNC $b\to q$ ($q=s,d$)
transitions is provided by the Wilson OPE: the $b\to q$ effective Hamiltonian
describing dynamics at the scale $\mu$, appropriate for $B$-decays, reads
Grinstein:1988me ; Burasa ; Burasb :

$\displaystyle H_{\rm eff}^{b\to
q}=\frac{G_{F}}{\sqrt{2}}V^{*}_{tq}V_{tb}\sum_{i}C_{i}(\mu){\cal O}_{i}^{b\to
q}(\mu),$ (2.1)

$G_{F}$ is the Fermi constant and $V_{ij}$ are CKM matrix elements. The SM
Wilson coefficients relevant for our analysis at the scale $\mu_{0}=5$ GeV
have the following values [corresponding to $C_{2}(M_{W})=-1$]:
$C_{1}(\mu_{0})=0.241$, $C_{2}(\mu_{0})=-1.1$, $C_{7}(\mu_{0})=0.312$ Burasa ;
Burasb ; Simulaa ; Simulab ; hidr .

The basis operators ${\cal O}_{i}^{b\to q}(\mu)$ 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}(\mu)$. 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.

###  2.2 The penguin contribution

The top-quark contribution to $B_{s}\to\gamma\gamma$ decay is generated by
penguin operator in (2.1)111 Our notations and conventions are:
$\gamma^{5}=i\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}$,
$\sigma_{\mu\nu}=\frac{i}{2}[\gamma_{\mu},\gamma_{\nu}]$,
$\varepsilon^{0123}=-1$,
$\epsilon_{abcd}\equiv\epsilon_{\alpha\beta\mu\nu}a^{\alpha}b^{\beta}c^{\mu}d^{\nu}$,
$e=\sqrt{4\pi\alpha_{\rm em}}$.

$H_{\text{eff}}^{b\to
s\,\gamma}=\frac{G_{F}}{\sqrt{2}}V_{tb}V_{ts}^{*}\,C_{7\gamma}(\mu){\cal
O}_{7\gamma},\qquad{\cal
O}_{7\gamma}=-\frac{e}{8\pi^{2}}m_{b}\cdot\bar{s}\sigma_{\mu\nu}(1+\gamma_{5})b\cdot
F^{\mu\nu}.$ (2.2)

The sign of the $b\to d\gamma$ 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

$iQ_{q}e\bar{q}\gamma_{\mu}q\epsilon^{\mu},$ (2.3)

corresponding to the definition of the covariant derivative in the form
$D_{\mu}=\partial_{\mu}-ie\,Q_{q}A_{\mu}$.

The amplitude of the $B\to\gamma\gamma$ transition is defined according to
mn2004 ; mnk2018 :

$\displaystyle{\cal A}^{(B\to\gamma\gamma)}_{\text{top}}$
$\displaystyle\equiv$
$\displaystyle\langle\gamma(q,\varepsilon),\gamma(q^{\prime},\varepsilon^{\prime})|H^{b\to
s\gamma}_{\rm eff}|\bar{B}_{s}(p)\rangle$ (2.4) $\displaystyle=$
$\displaystyle-2\,\frac{G_{F}}{\sqrt{2}}V_{tb}V_{ts}^{*}\,\frac{e^{2}}{8\pi^{2}}\,2m_{b}C_{7\gamma}(\mu)\Big{[}F_{TV}\,\epsilon_{\alpha\alpha^{\prime}qq^{\prime}}-iF_{TA}\,(g_{\alpha\alpha^{\prime}}q^{\prime}q-q^{\prime}_{\alpha}q_{\alpha^{\prime}})\Big{]}\varepsilon_{\alpha}\varepsilon^{\prime}_{\alpha^{\prime}}.$

Here $q,q^{\prime}$ and $\varepsilon,\varepsilon^{\prime}$ are momenta and
polarization vectors of the outgoing real photons, and $F_{TA}$ and $F_{TV}$
are the form factors $F_{TA}(q^{2}=0,q^{\prime 2}=0)$ and
$F_{TV}(q^{2}=0,q^{\prime 2}=0)$. The latter are defined as mk2003 ; mn2004 ;
mnk2018 :

$\displaystyle\langle\gamma(q^{\prime},\varepsilon^{\prime})|\bar{s}\sigma_{\mu\nu}\gamma_{5}b|\bar{B}_{s}(p)\rangle\,q^{\nu}=e\,\varepsilon^{\prime}_{\alpha}\,\left(g_{\mu\alpha}\,q^{\prime}q-q_{\alpha}q^{\prime}_{\mu}\right)\,F_{TA}(q^{2},q^{\prime
2}),$ (2.5)
$\displaystyle\langle\gamma(q^{\prime},\varepsilon^{\prime})|\bar{s}\sigma_{\mu\nu}b|\bar{B}_{s}(p)\rangle\,q^{\nu}=i\,e\,\varepsilon^{\prime}_{\alpha}\epsilon_{\mu\alpha
qq^{\prime}}F_{TV}(q^{2},q^{\prime 2}),$ (2.6)

and satisfy a rigorous constraint $F_{TA}(q^{2},0)=F_{TV}(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_{TA}$ and $F_{TV}$ mnk2018 .

###  2.3 Nonfactorizable charm-quark loop correction to $B_{s}\to\gamma\gamma$

As already noticed, the light degrees of freedom remain dynamical and their
contributions should be taken into account separately. The relevant terms in
$H_{\rm eff}^{b\to s}$ are those containing four-quark operators:

$H_{\rm eff}^{b\to
s\bar{c}c}=-\frac{G_{F}}{\sqrt{2}}\,V_{cb}V^{*}_{cs}\,\left\\{C_{1}(\mu){\cal
O}_{1}+C_{2}(\mu){\cal O}_{2}\right\\}$ (2.7)

where

$\displaystyle{\cal
O}_{1}=\left(\bar{s}^{i}\gamma_{\mu}(1-\gamma^{5})c^{j}\right)\;\left(\bar{c}^{k}\gamma^{\mu}(1-\gamma^{5})b^{l}\right)\;\delta_{il}\delta_{kj},$
(2.8) $\displaystyle{\cal
O}_{2}=\left(\bar{s}^{i}\gamma_{\mu}(1-\gamma^{5})c^{j}\right)\;\left(\bar{c}^{k}\gamma^{\mu}(1-\gamma^{5})b^{l}\right)\;\delta_{ij}\delta_{kl},$
(2.9)

differing from each other in the way color indices $i,j,k,l$ are contracted.
By Fierz transformation ${\cal O}_{2}$ may be written in the following form
(for anticommuting spinor fields):

${\cal
O}_{2}=\left(\bar{s}^{i}\gamma_{\mu}(1-\gamma^{5})b^{l}\right)\left(\bar{c}^{k}\gamma^{\mu}(1-\gamma^{5})c^{j}\right)\left(2t^{a}_{il}t^{a}_{jk}+\frac{1}{3}\delta_{il}\delta_{jk}\right),$
(2.10)

The singlet-singlet operator ${\cal O}_{1}$ and the singlet-singlet part of
${\cal O}_{2}$ at the leading order generate factorizable charm contributions
to the $B\to\gamma\gamma$ 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 ${\cal O}_{2}$ and needs the emission of one
soft gluon from the charm-quark loop. So relevant for us is the octet-octet
operator

$H_{\text{eff}}^{b\to s\bar{c}c[8\times
8]}=-\frac{G_{F}}{\sqrt{2}}V_{cb}V^{*}_{cs}2C_{2}\left(\overline{s}\gamma_{\mu}(1-\gamma_{5})t^{a}b\right)\left(\overline{c}\gamma^{\mu}(1-\gamma_{5})t^{a}c\right).$
(2.11)

Therefore, similar to the top contribution, we find

$\displaystyle{\cal
A}^{(B\to\gamma\gamma)}_{\text{charm}}=\langle\gamma(q,\varepsilon),\gamma(q^{\prime},\varepsilon^{\prime})|H^{b\to
s\bar{c}c[8\times 8]}_{\rm eff}|\bar{B}_{s}(p)\rangle.$ (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

$\displaystyle{\cal
A}^{(B\to\gamma\gamma)}_{\text{charm}}=\frac{1}{2}i^{3}\langle\gamma(q,\varepsilon),\gamma(q^{\prime},\varepsilon^{\prime})|$
$\displaystyle\hskip 28.45274pt\times T\Big{\\{}H^{b\to s\bar{c}c[8\times
8]}_{\rm eff}(0),\int dzj^{\rm e.m.}_{\rho}(z)A_{\rho}(z),\int dyj^{\rm
e.m.}_{\eta}(y)A_{\eta}(y),\int
dx\bar{c}(x)\gamma_{\nu}t^{b}c(x)g_{s}B_{\nu}^{b}(x)\Big{\\}}|\bar{B}_{s}(p)\rangle$
(2.13)

where quark electormagnetic current has the form $j^{e.m.}_{\alpha}=e\sum
Q_{i}\bar{q}_{i}\gamma_{\alpha}q_{i}$.

Fig. 1: One of the diagrams describing charming loop contribution to
$B_{s}\to\gamma\,\gamma$ decay via nonfactorizable soft gluon exchange. Other
diagrams are those corresponding to an opposite direction in the charm-quark
loop and diagrams with the interchanged photons $q\leftrightarrow
q^{\prime},\varepsilon\leftrightarrow\varepsilon^{\prime}$.

Eq. (2.3) may be rewritten as:

$\displaystyle{\cal A}^{(B\to\gamma\gamma)}_{\text{charm}}$ $\displaystyle=$
$\displaystyle-i^{3}\frac{G_{F}}{\sqrt{2}}2C_{2}V_{cb}V_{cs}^{*}e^{2}Q_{c}Q_{s}\int
dzdxdy\left[\varepsilon_{\rho}e^{iqz}\varepsilon^{\prime}_{\eta}e^{iq^{\prime}y}+(q\leftrightarrow
q^{\prime},\varepsilon\leftrightarrow\varepsilon^{\prime})\right]$ (2.14)
$\displaystyle\times\langle
0|T\Big{\\{}\bar{c}\gamma_{\rho}c(z),\bar{c}(0)t^{a}\gamma_{\mu}(1-\gamma_{5})c(0),\bar{c}(x)t^{b}\gamma_{\nu}c(x)\Big{\\}}|0\rangle$
$\displaystyle\times\langle
0|T\Big{\\{}\bar{s}(y)\gamma_{\eta}s(y),\bar{s}(0)t^{a}\gamma_{\mu}(1-\gamma_{5})b(0)g_{s}B^{b}_{\nu}(x)\Big{\\}}|\bar{B}_{s}(p)\rangle.$

Fig. 1 shows 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.

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 $\bar{c}\gamma_{\mu}c$ part of the $V-A$ weak current does not contribute
and one is left with $\langle VVA\rangle$ charm-quark triangle.

  * •

The $\langle VVA\rangle$ charm-quark triangle contracted with the gluon field
$B^{b}_{\nu}$ may be written as a gauge-invariant nonlocal operator containing
gluon field strength $G^{b}_{\nu\alpha}$ for any gluon momentum (cf. buchalla
).

Making use of the result for the charm-quark $VVA$ triangle from the next
Section, we obtain the following expression for the amplitude:

$\displaystyle A^{B\to\gamma\gamma}_{\rm charm}$ $\displaystyle=$
$\displaystyle\frac{G_{F}}{\sqrt{2}}4C_{2}V_{cb}V^{*}_{cs}e^{2}Q_{s}Q_{c}\,A_{\rho\eta}(q,q^{\prime})\varepsilon_{\rho}\varepsilon^{\prime}_{\eta},$
(2.15) $\displaystyle A_{\rho\eta}(q,q^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{8}}\int dkdye^{-i(k-q^{\prime})y}dxd\kappa
e^{-i\kappa x}\Gamma_{cc}^{\mu\nu\rho(ab)}(\kappa,q)\langle
0|\bar{s}(y)\gamma^{\eta}\frac{\not{k}+m_{s}}{m_{s}^{2}-k^{2}}\gamma^{\mu}(1-\gamma^{5})t^{a}B^{b}_{\nu}(x)b(0)|\bar{B}_{s}(p)\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{4(2\pi)^{8}}\int
dkdye^{-i(k-q^{\prime})y}dxd\kappa e^{-i\kappa
x}\overline{\Gamma}_{cc}^{\mu\nu\rho\alpha}(\kappa,q)\langle
0|\bar{s}(y)\gamma^{\eta}\frac{\not{k}+m_{s}}{m_{s}^{2}-k^{2}}\gamma^{\mu}(1-\gamma^{5})t^{b}G^{b}_{\nu\alpha}(x)b(0)|\bar{B}_{s}(p)\rangle.$

For real photons in the final state, the amplitude
$A_{\rho\eta}(q,q^{\prime})$ has the same Lorentz structure as the penguin
amplitude and contains two form factors $H_{V}=H_{V}(q^{2}=0,q^{\prime 2}=0)$
and $H_{A}=H_{A}(q^{2}=0,q^{\prime 2}=0)$ (Appendix A presents comparison with
the form factors defined in mnk2018 ):

$A_{\rho\eta}(q,q^{\prime})=H_{V}\epsilon_{\rho\eta
qq^{\prime}}-iH_{A}\left(g_{\rho\eta}\,qq^{\prime}-q^{\prime}_{\rho}q_{\eta}\right),$
(2.17)

Comparing Eqs. (2.4) and (2.17), and taking into account that
$V_{tb}V^{*}_{ts}\simeq-V_{cb}V^{*}_{cs}$, it is convenient to desctribe the
effect of charm as an additions to the Wilson coefficient $C_{7\gamma}$
(however, non-universal, i.e. different in the axial and vector Lorentz
structures):

$\displaystyle\epsilon_{\rho\eta qq^{\prime}}:$ $\displaystyle C_{7\gamma}\to
C_{7\gamma}(1+\delta_{V}C_{7\gamma}),$ $\displaystyle
g_{\rho\eta}\,qq^{\prime}-q^{\prime}_{\rho}q_{\eta}:$ $\displaystyle
C_{7\gamma}\to C_{7\gamma}(1+\delta_{A}C_{7\gamma}),$ (2.18)

with

$\displaystyle\delta_{V(A)}C_{7\gamma}=8\pi^{2}\,Q_{s}Q_{c}\frac{C_{2}}{C_{7\gamma}}\frac{H_{V(A)}}{m_{b}\,F_{TV(TA)}}.$
(2.19)

Our goal will be to calculate these corrections.

The $B$-meson structure contributes to the ${\cal A}^{B\to\gamma\gamma}_{\rm
charm}$ amplitude via the full set of 3BS

$\displaystyle\langle
0|\bar{s}(y)\Gamma_{i}t^{a}b(0)G^{a}_{\nu\alpha}(x)|\bar{B}_{s}(p)\rangle,$
(2.20)

with $\Gamma_{i}$ the appropriate combinations of $\gamma$-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 mk2018 , the dominant contribution of charm
to amplitudes of FCNC $B$ decays comes from the “double collinear” LC
configuration m2023 , where $x^{2}=0$, $y^{2}=0$, but $xy\neq 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.

##  3 Charm-quark $\langle VVA\rangle$ triangle

The charm-quark loop contribution is described by the three-point function
(see Fig. 2):

$\Gamma_{cc}^{\mu\nu\rho\,(ab)}\left(\kappa,q\right)=\int
dx^{\prime}dz\>e^{iqz+i\kappa x^{\prime}}\langle
0|T\\{\bar{c}(z)\gamma^{\rho}c(z),\bar{c}(0)\gamma^{\mu}(1-\gamma_{5})t^{a}c(0),\bar{c}(x^{\prime})\gamma^{\nu}t^{b}c(x^{\prime})\\}|0\rangle=\frac{1}{2}\delta^{ab}\,\Gamma_{cc}^{\mu\nu\rho}(\kappa,q),$
(3.21)

where $q$ is the momentum of the external virtual photon (vertex containing
index $\rho$) and $\kappa$ is the gluon momentum (vertex containing index
$\nu$). Here $t^{c}$, $c=1,\dots,8$ are $SU_{c}(3)$ generators normalized as
${\rm Sp}(t^{a}t^{b})=\frac{1}{2}\delta^{ab}$.

Fig. 2: The $\langle VVA\rangle$ triangle one-loop diagrams for
$\Gamma_{cc}^{\mu\nu\rho\,(ab)}$.

The octet current $\bar{c}(0)\gamma^{\mu}(1-\gamma_{5})t^{a}c(0)$ is a charm-
quark part of the octet-octet weak Hamiltonian. Its vector piece does not
contribute to $\Gamma_{cc}^{\mu\nu\rho\,(ab)}$ (Furry theorem) and will be
omitted. Taking into account vector-current conservation, it is convenient to
parametrize $\Gamma_{cc}^{\mu\nu\rho}(\kappa,q)$ as follows lm

$\Gamma_{cc}^{\mu\nu\rho}(\kappa,q)=-i\left(\kappa^{\mu}+q^{\mu}\right)\epsilon^{\nu\rho\kappa
q}\,F_{0}-i\left(q^{2}\epsilon^{\mu\nu\rho\kappa}-q^{\rho}\epsilon^{\mu\nu
q\kappa}\right)\,F_{1}-i\left(\kappa^{2}\epsilon^{\mu\rho\nu
q}-\kappa^{\nu}\epsilon^{\mu\rho\kappa q}\right)\,F_{2}.$ (3.22)

The form factors $F_{0,1,2}$ are functions of three independent invariant
variables $q^{2}$, $\kappa^{2}$, and $\kappa q$. The lowest order QCD diagrams
describing $\Gamma_{cc}^{\mu\nu\rho}(\kappa,q)$ are shown in Fig. 2. A
convenient representation of the form factors has the form max

$\displaystyle F_{i}\left(\kappa^{2},\kappa q,q^{2}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}}\int\limits_{0}^{1}d\xi\int\limits_{0}^{1-\xi}d\eta\frac{\Delta_{i}(\xi,\eta)}{m_{c}^{2}-2\,\xi\eta\>\kappa
q-\xi(1-\xi)q^{2}-\eta(1-\eta)\kappa^{2}},\quad i=0,1,2,$ (3.23)
$\displaystyle\Delta_{0}=-\xi\eta,\quad\Delta_{1}=\xi(1-\eta-\xi),\quad\Delta_{2}=\eta(1-\eta-\xi).$

This representation may be applied to the physical amplitude in the region of
the external momenta far below the thresholds,
$q^{2},\kappa^{2},(\kappa+q)^{2}\ll 4m_{c}^{2}$. 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 $\kappa^{2}\sim\Lambda^{2}_{\rm QCD}$, $(q+\kappa)^{2}<0$. So, the
representation (3.23) is applicable and proves convenient for numerical
calculations.

As the next step, one takes the convolution of the amplitude (3.22) with the
gluon field $B_{\nu}(x)$. The Lorentz structures multiplying $F_{0}$ and
$F_{1}$ contain $\epsilon_{\nu\kappa\alpha_{1}\alpha_{2}}$ with some indices
$\alpha_{1}$ and $\alpha_{2}$. After multiplying by $B_{\nu}(x)$ and
performing parts integration, their contribution may be reduced to the
convolution with the gluon strength tensor $G_{\alpha\nu}$:

$\displaystyle\int e^{-i\kappa
x}\epsilon^{\alpha\nu\rho\omega}\kappa_{\alpha}B_{\nu}(x)dx$ $\displaystyle=$
$\displaystyle i\int\left[\frac{\partial}{\partial x_{\alpha}}e^{-i\kappa
x}\right]\epsilon^{\alpha\nu\rho\omega}B_{\nu}(x)dx=-\frac{i}{2}\int
e^{-i\kappa x}\epsilon^{\alpha\nu\rho\omega}\left[\frac{\partial}{\partial
x_{\alpha}}B_{\nu}(x)-\frac{\partial}{\partial x_{\nu}}B_{\alpha}(x)\right]dx$
(3.24) $\displaystyle=$ $\displaystyle-\frac{i}{2}\int e^{-i\kappa
x}\epsilon^{\alpha\nu\rho\omega}G_{\alpha\nu}(x)dx.$

The Lorentz structure multiplying $F_{2}$ at first glance does not have this
property. However, using the identity

$\displaystyle
g^{\alpha_{1}\alpha_{2}}\epsilon^{\alpha_{3}\alpha_{4}\alpha_{5}\alpha_{6}}-g^{\alpha_{1}\alpha_{3}}\epsilon^{\alpha_{2}\alpha_{4}\alpha_{5}\alpha_{6}}+g^{\alpha_{1}\alpha_{4}}\epsilon^{\alpha_{2}\alpha_{3}\alpha_{5}\alpha_{6}}-g^{\alpha_{1}\alpha_{5}}\epsilon^{\alpha_{2}\alpha_{3}\alpha_{4}\alpha_{6}}+g^{\alpha_{1}\alpha_{6}}\epsilon^{\alpha_{2}\alpha_{3}\alpha_{4}\alpha_{5}}=0,$
(3.25)

multiplying it by $\kappa_{\alpha_{1}}\kappa_{\alpha_{2}}q_{\alpha_{6}}$, and
setting $\alpha_{3}\to\mu$, $\alpha_{4}\to\nu$, $\alpha_{5}\to\rho$, this
Lorentz structure takes the form

$\displaystyle\kappa^{2}\epsilon^{\mu\nu\rho
q}+\kappa^{\nu}\epsilon^{\mu\rho\kappa q}=\kappa^{\mu}\epsilon^{\kappa\nu\rho
q}+\kappa^{\rho}\epsilon^{\kappa\mu\nu q}-\kappa q\epsilon^{\kappa\mu\nu\rho}$
(3.26)

and may be also reduced to the convolution with $G_{\alpha\nu}$. Finally, the
operator describing the contribution of the charm-quark loop takes the form

$\displaystyle\int d\kappa e^{-i\kappa
x~{}}\Gamma_{cc}^{\mu\nu\rho\,(ab)}(\kappa,q)B^{b}_{\nu}(x)dx=\frac{1}{4}\int
d\kappa e^{-i\kappa
x~{}}{\overline{\Gamma}}_{cc}^{\mu\nu\rho\alpha}(\kappa,q)G^{a}_{\nu\alpha}(x)dx$
(3.27)

with

$\displaystyle{\overline{\Gamma}}_{cc}^{\mu\nu\rho\alpha}(\kappa,q)=\left(\kappa^{\mu}+q^{\mu}\right)\epsilon^{\nu\rho\alpha
q}\,F_{0}+\left(q^{\rho}\epsilon^{\mu\nu\alpha
q}+q^{2}\epsilon^{\mu\nu\rho\alpha}\right)\,F_{1}+\left(\kappa^{\mu}\epsilon^{\alpha\nu\rho
q}+\kappa^{\rho}\epsilon^{\alpha\mu\nu q}-\kappa
q\,\epsilon^{\alpha\mu\nu\rho}\right)\,F_{2}.$ (3.28)

${\overline{\Gamma}}_{cc}^{\mu\nu\rho\alpha}(\kappa,q)$ is real in the
Euclidean region and
$\Gamma_{cc}^{\mu\nu\rho}(\kappa,q)=-i\,{\overline{\Gamma}}_{cc}^{\mu\nu\rho\alpha}(\kappa,q)\kappa_{\alpha}$.

##  4 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\to\pi,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 m2022 : 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 $\bar{c}c$-pair moves along the other LC direction.
Introducing vectors $n_{\mu}$ and $n^{\prime}_{\mu}$ such that
$n^{2}=n^{\prime 2}=0$, $n^{\prime}n=2$,
$v_{\mu}=p_{\mu}/M_{B}=\frac{1}{2}(n_{\mu}+n^{\prime}_{\mu})$, one finds that
the dominant contribution of charming loops to an FCNC $B$-decay amplitude
comes from the double-collinear configuration mk2018 ; m2019 ; m2022 ; m2023 ;
wang2022  when the coordinates of the field operators in
$\langle\bar{s}(y)G_{\mu\nu}(x)b(0)|\bar{B}_{s}(p)\rangle$ are alligned along
the orthogonal light-cone directions $x_{\mu}\sim n_{\mu}$, $y_{\mu}\sim
n^{\prime}_{\mu}$. The 3BS amplitude in the collinear and the double-collinear
kinematics contain the same Lorentz structures wang2022b  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}\to\gamma\gamma$ amplitude calculated using
the general noncollinear 3BS differs by terms $O(\lambda_{B_{s}}/M_{B})$ from
the amplitude calculated within the double-collinear approximation.

###  4.1 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.

####  4.1.1 The Lorentz structure of the collinear 3BS

We start with the collinear LC 3BS braun2017 , where the arguments of the
$s$-quark field, $\bar{s}(y)$, and the gluon field $G_{\nu\alpha}(x)$ are
collinear to each other, $x=uy$, $u\neq 0$ is a number (in this case $x^{2}=0$
leads to $y^{2}=0$):

$\displaystyle\langle
0|\bar{s}(y)G_{\nu\alpha}(u\,y)\Gamma\,b(0)|\bar{B}_{s}(p)\rangle=\frac{f_{B}M_{B}^{3}}{4}\int
D(\omega,\lambda)\,e^{-i\lambda yp-i\omega uyp}{\rm
Sp}\Bigg{\\{}\gamma_{5}\Gamma(1+\not{v})$ $\displaystyle\hskip
14.22636pt\times\bigg{[}(p_{\nu}\gamma_{\alpha}-p_{\alpha}\gamma_{\nu})\frac{1}{M_{B}}[\Psi_{A}-\Psi_{V}]-i\sigma_{\nu\alpha}\Psi_{V}-\frac{(y_{\nu}p_{\alpha}-y_{\alpha}p_{\nu})}{yp}\left(X_{A}+\frac{\not{y}}{yp}M_{B}W\right)$
$\displaystyle\hskip
28.45274pt+\frac{(y_{\nu}\gamma_{\alpha}-y_{\alpha}\gamma_{\nu})}{yp}M_{B}\left(Y_{A}+W+\frac{\not{x}}{yp}M_{B}Z\right)-i\epsilon_{\nu\alpha\mu\beta}\,\frac{y^{\mu}p^{\beta}}{yp}\gamma^{5}\tilde{X}_{A}+i\epsilon_{\nu\alpha\mu\beta}\,\frac{y^{\mu}\gamma^{\beta}}{yp}\gamma^{5}M_{B}\tilde{Y}_{A}\bigg{]}\Bigg{\\}},$
(4.29)

where $D(\omega,\lambda)$ takes into account rigorous constraints on the
variables $\omega$ and $\lambda$:

$\displaystyle
D(\omega,\lambda)=d\omega\,d\lambda\,\theta(\omega)\theta(\lambda)\theta(1-\omega-\lambda).$
(4.30)

In Eq. (4.1.1), $\Gamma$ is an arbitrary combination of Dirac matrices,
$v_{\mu}=p_{\mu}/M_{B}$, and all 8 DAs ($\Psi_{A}$, $\Psi_{V}$, etc) are
functions of two dimensionless arguments $0<\lambda<1$ and $0<\omega<1$. Here
$\lambda$ refers to the momentum carried by the $s$-quark, and $\omega$ refers
to the momentum carried by the gluon.

The normalization conditions for $\Psi_{A}$ and $\Psi_{V}$ have the form
braun2017 :

$\displaystyle\int
D(\omega,\lambda)\,\Psi_{A}\left(\omega,\lambda\right)=\frac{\lambda_{E}^{2}}{3M_{B}^{2}},\qquad\int
D(\omega,\lambda)\,\Psi_{V}\left(\omega,\lambda\right)=\frac{\lambda_{H}^{2}}{3M_{B}^{2}}.$
(4.31)

Some of the Lorentz structures in (4.1.1) contain factors $x_{\mu}/xp$ or
$x_{\mu}x_{\nu}/(xp)^{2}$. Since 3BS (4.1.1) is a continuous regular function
at $x^{2}=0$ and $xp=0$, the absence of singularities at $xp\to 0$ leads to
the following constraints:

$\displaystyle\int
D(\omega,\lambda)\,\big{\\{}X_{A},Y_{A},\tilde{X}_{A},\tilde{Y}_{A},Z,W\big{\\}}=0,$
$\displaystyle\int
D(\omega,\lambda)\,\omega\,\big{\\{}Z,W\big{\\}}=0,\quad\int
D(\omega,\lambda)\,\lambda\,\big{\\{}Z,W\big{\\}}=0.$ (4.32)

These constraints are obtained by expanding the exponential in the integral
representation (4.1.1) under the condition $x_{\mu}=u\,y_{\mu}$ to the
necessary order and requiring that the coefficients multiplying terms singular
in $xp\to 0$ vanish.

####  4.1.2 Twist expansion of the 3DAs

The DAs in (4.1.1) have no definite twist. According to braun2017 , the
distribution amplitudes might be written as an expansion in functions with
definite twist as follows:

$\displaystyle\Psi_{A}(\omega,\lambda)$ $\displaystyle=$
$\displaystyle(\phi_{3}+\phi_{4})/2,$ $\displaystyle\Psi_{V}(\omega,\lambda)$
$\displaystyle=$ $\displaystyle(-\phi_{3}+\phi_{4})/2,$
$\displaystyle{X}_{A}(\omega,\lambda)$ $\displaystyle=$
$\displaystyle(-\phi_{3}-\phi_{4}+2\psi_{4})/2,$
$\displaystyle{Y}_{A}(\omega,\lambda)$ $\displaystyle=$
$\displaystyle(-\phi_{3}-\phi_{4}+\psi_{4}-\psi_{5})/2,$
$\displaystyle\tilde{X}_{A}(\omega,\lambda)$ $\displaystyle=$
$\displaystyle(-\phi_{3}+\phi_{4}-2\widetilde{\psi}_{4})/2,$
$\displaystyle\tilde{Y}_{A}(\omega,\lambda)$ $\displaystyle=$
$\displaystyle(-\phi_{3}+\phi_{4}-\widetilde{\psi}_{4}+\widetilde{\psi}_{5})/2,$
$\displaystyle W(\omega,\lambda)$ $\displaystyle=$
$\displaystyle(\phi_{4}-\psi_{4}-\widetilde{\psi}_{4}+\phi_{5}+\psi_{5}+\widetilde{\psi}_{5})/2,$
$\displaystyle{Z}(\omega,\lambda)$ $\displaystyle=$
$\displaystyle(-\phi_{3}+\phi_{4}-2\widetilde{\psi}_{4}+\phi_{5}+2\widetilde{\psi}_{5}-\phi_{6})/4,$
(4.33)

where we keep the contributions up to twist 6 inclusively (the subscript “$i$”
in $\phi_{i}$ and $\psi_{i}$ denote the twist value).

####  4.1.3 Model for DAs entering the collinear 3BS

The powers of $\omega$ and $\lambda$ 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
braun2017  for twist 3- and 4, complemented by twist 5 and 6 DAs reconstructed
using the constraints (4.1.1) wang2019 :

$\displaystyle\phi_{3}$ $\displaystyle=$
$\displaystyle\frac{105(\lambda_{E}^{2}-\lambda_{H}^{2})}{32\omega_{0}^{7}M_{B}^{2}}\,\lambda\omega^{2}\left(2\omega_{0}-\omega-\lambda\right)^{2}\theta\left(2\omega_{0}-\omega-\lambda\right),$
(4.34) $\displaystyle\phi_{4}$ $\displaystyle=$
$\displaystyle\frac{35(\lambda_{E}^{2}+\lambda_{H}^{2})}{32\omega_{0}^{7}M_{B}^{2}}\,\omega^{2}\left(2\omega_{0}-\omega-\lambda\right)^{3}\theta\left(2\omega_{0}-\omega-\lambda\right),$
(4.35) $\displaystyle\psi_{4}$ $\displaystyle=$
$\displaystyle\frac{35\lambda_{E}^{2}}{16\omega_{0}^{7}M_{B}^{2}}\,\lambda\omega\left(2\omega_{0}-\omega-\lambda\right)^{3}\theta\left(2\omega_{0}-\omega-\lambda\right),$
(4.36) $\displaystyle\widetilde{\psi}_{4}$ $\displaystyle=$
$\displaystyle\frac{35\lambda_{H}^{2}}{16\omega_{0}^{7}M_{B}^{2}}\,\lambda\omega\left(2\omega_{0}-\omega-\lambda\right)^{3}\theta\left(2\omega_{0}-\omega-\lambda\right),$
(4.37) $\displaystyle\phi_{5}$ $\displaystyle=$
$\displaystyle\frac{35\left(\lambda_{E}^{2}+\lambda_{H}^{2}\right)}{64\omega_{0}^{7}M_{B}^{2}}\,\lambda\left(2\omega_{0}-\omega-\lambda\right)^{4}\theta\left(2\omega_{0}-\omega-\lambda\right),$
(4.38) $\displaystyle\psi_{5}$ $\displaystyle=$
$\displaystyle-\frac{35\lambda_{E}^{2}}{64\omega_{0}^{7}M_{B}^{2}}\,\omega\left(2\omega_{0}-\omega-\lambda\right)^{4}\theta\left(2\omega_{0}-\omega-\lambda\right),$
(4.39) $\displaystyle\widetilde{\psi}_{5}$ $\displaystyle=$
$\displaystyle-\frac{35\lambda_{H}^{2}}{64\omega_{0}^{7}M_{B}^{2}}\,\omega\left(2\omega_{0}-\omega-\lambda\right)^{4}\theta\left(2\omega_{0}-\omega-\lambda\right),$
(4.40) $\displaystyle\phi_{6}$ $\displaystyle=$
$\displaystyle\frac{7\left(\lambda_{E}^{2}-\lambda_{H}^{2}\right)}{64\omega_{0}^{7}M_{B}^{2}}\,\left(2\omega_{0}-\omega-\lambda\right)^{5}\theta\left(2\omega_{0}-\omega-\lambda\right).$
(4.41)

Dimensionless parameter $\omega_{0}$ is related to $\lambda_{B}$, the inverse
moment of the $B$-meson LC distribution amplitude, as

$\displaystyle\omega_{0}=\frac{5}{2}\frac{\lambda_{B}}{m_{B}}.$ (4.42)

For this model, the integration limits take the following form
($2\omega_{0}<1$):

$\displaystyle\int
D(\omega,\lambda)\,\theta\left(2\omega_{0}-\omega-\lambda\right)\,{(\dots)}=\int\limits_{0}^{2\omega
0}d\omega\int\limits_{0}^{2\omega_{0}-\omega}d\lambda\,{(\dots)}.$ (4.43)

However, as we shall show shortly, certain modifications of the integrated DAs
[which emerge when performing parts integrations of the 3BS (4.1.1)] at large
values of $\omega$ and $\lambda$ will be necessary in order to satisfy the
continuity of the 3BS considered in a noncollinear kinematics.

###  4.2 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:222We do not include here those
structures that vanish in the collinear limit $x=uy$, such as e.g.
$(x_{\sigma}y_{\nu}-x_{\nu}y_{\sigma})/xy$. We also do not consider structures
of the type $1/(xp\,yp)$ that may emerge when generalizing the Lorentz
structures multiplying $W$ and $Z$ DAs in (4.1.1); according to our analysis
the $W$ and $Z$-structures anyway give a marginal contribution to the FCNC
$B$-decay amplitude.

$\displaystyle\langle
0|\bar{s}(y)G_{\nu\alpha}(x)\Gamma\,b(0)|\bar{B}_{s}(p)\rangle=\frac{f_{B}M_{B}^{3}}{4}\int
D(\omega,\lambda)\,e^{-i\lambda yp-i\omega xp}\,{\rm
Sp}\Bigg{\\{}\gamma_{5}\Gamma\,(1+\not{v})$
$\displaystyle\times\bigg{[}(p_{\nu}\gamma_{\alpha}-p_{\alpha}\gamma_{\nu})\frac{1}{M_{B}}[\Psi_{A}-\Psi_{V}]-i\sigma_{\nu\alpha}\Psi_{V}$
$\displaystyle-\frac{(x_{\nu}p_{\alpha}-x_{\alpha}p_{\nu})}{xp}\left(X_{A}^{(x)}+\frac{\not{x}}{xp}M_{B}W^{(x)}\right)+\frac{(x_{\nu}\gamma_{\alpha}-x_{\alpha}\gamma_{\nu})}{xp}M_{B}\left(Y^{(x)}_{A}+W^{(x)}+\frac{\not{x}}{xp}M_{B}Z^{(x)}\right)$
$\displaystyle-\frac{(y_{\nu}p_{\alpha}-y_{\alpha}p_{\nu})}{yp}\left(X_{A}^{(y)}+\frac{\not{y}}{yp}M_{B}W^{(y)}\right)+\frac{(y_{\nu}\gamma_{\alpha}-y_{\alpha}\gamma_{\nu})}{yp}M_{B}\left(Y^{(y)}_{A}+W^{(y)}+\frac{\not{y}}{yp}M_{B}Z^{(y)}\right)$
$\displaystyle-i\epsilon_{\nu\alpha\mu\beta}\,\frac{x^{\mu}p^{\beta}}{xp}\gamma^{5}\tilde{X}^{(x)}_{A}+i\epsilon_{\nu\alpha\mu\beta}\,\frac{x^{\mu}\gamma^{\beta}}{xp}\gamma^{5}M_{B}\tilde{Y}^{(x)}_{A}-i\epsilon_{\nu\alpha\mu\beta}\,\frac{x^{\mu}p^{\beta}}{xp}\gamma^{5}\tilde{X}^{(y)}_{A}+i\epsilon_{\nu\alpha\mu\beta}\,\frac{x^{\mu}\gamma^{\beta}}{xp}\gamma^{5}M_{B}\tilde{Y}^{(y)}_{A}\bigg{]}\Bigg{\\}}.$
(4.44)

All invariant amplitudes $\Phi=\Psi_{A},\Psi_{V},\dots$ are functions of 5
variables, $\Phi(\omega,\lambda,x^{2},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 $\Phi$’s are
functions of dimensionless arguments $\lambda$ and $\omega$ and are referred
to as the DAs. These DAs contain at least the kinematical constraint
$\theta(1-\omega-\lambda)$. However, the DAs may have support in more
restricted areas: e.g., the DAs of the LD model Eqs. (4.34)–(4.41) have
support in the region $\theta(2\omega_{0}-\omega-\lambda)$, $2\omega_{0}<1$.

Obviously, the functions $\Psi_{A}$ and $\Psi_{V}$ in (4.1.1) and (4.2) are
the same. Other DAs in (4.1.1) and (4.2) are related to each other as follows:

$\displaystyle X=X^{(x)}+X^{(y)},\quad
X=\\{X_{A},Y_{A},\tilde{X}_{A},\tilde{Y}_{A}\\},\quad W=W^{(x)}+W^{(y)},\quad
Z=Z^{(x)}+Z^{(y)}.$ (4.45)

The amplitude (4.2) 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
$\theta(2\omega_{0}-\omega-\lambda)$, $2\omega_{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.2) to zero order and requiring that the
singular terms vanish:

$\displaystyle\int_{0}^{2\omega_{0}-\lambda}d\omega\,X^{(x)}(\omega,\lambda)=0\qquad\forall\lambda,$
$\displaystyle\int_{0}^{2\omega_{0}-\omega}d\lambda\,X^{(y)}(\omega,\lambda)=0\qquad\forall\omega.$
(4.46)

Let us introduce the primitives

$\displaystyle\overline{X}^{(x)}(\omega,\lambda)=\int_{0}^{\omega}d\omega^{\prime}{X}^{(x)}(\omega^{\prime},\lambda),$
(4.47)
$\displaystyle\overline{X}^{(y)}(\omega,\lambda)=\int_{0}^{\lambda}d\lambda^{\prime}{X}^{(y)}(\omega,\lambda^{\prime}),$
(4.48)

which, by virtue of (4.46), vanish at the boundary of the DA support region:

$\displaystyle\overline{X}^{(x)}(0,\lambda)=\overline{X}^{(x)}(2\omega_{0}-\lambda,\lambda)=0\quad\forall\lambda,$
$\displaystyle\overline{X}^{(y)}(\omega,0)=\overline{X}^{(y)}(\omega,2\omega_{0}-\omega)=0\quad\forall\omega.$
(4.49)

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:

$\displaystyle\int\limits_{0}^{2\omega_{0}-\lambda}d\omega\,\omega^{n}Z^{(x)}(\omega,\lambda)=0\quad\forall\lambda,\quad\int\limits_{0}^{2\omega_{0}-\omega}d\lambda\,\lambda^{n}Z^{(y)}(\omega,\lambda)=0\quad\forall\omega,\quad
n=0,1.$ (4.50)

By introducing primitives and double primitives

$\displaystyle\overline{Z}^{(x)}(\omega,\lambda)=\int\limits_{0}^{\omega}d\omega^{\prime}Z^{(x)}(\omega^{\prime},\lambda),\qquad\overline{\overline{Z}}^{(x)}(\omega,\lambda)=\int\limits_{0}^{\omega}d\omega^{\prime}\overline{Z}^{(x)}(\omega^{\prime},\lambda)$
$\displaystyle\overline{Z}^{(y)}(\omega,\lambda)=\int\limits_{0}^{\lambda}d\lambda^{\prime}Z^{(y)}(\omega,\lambda^{\prime}),\qquad\overline{\overline{Z}}^{(y)}(\omega,\lambda)=\int\limits_{0}^{\lambda}d\omega^{\prime}\overline{Z}^{(y)}(\omega,\lambda^{\prime}),$
(4.51)

one can show that the requirements (4.50) lead to the vanishing of these
functions at the boundaries of the DAs support regions:

$\displaystyle\overline{Z}^{(x)}(0,\lambda)=\overline{Z}^{(x)}(2\omega_{0}-\lambda,\lambda)=0,\qquad\overline{\overline{Z}}^{(x)}(0,\lambda)=\overline{\overline{Z}}^{(x)}(2\omega_{0}-\lambda,\lambda)=0\qquad\forall\lambda,$
$\displaystyle\overline{Z}^{(y)}(\omega,0)=\overline{Z}^{(y)}(\omega,2\omega_{0}-\omega)=0,\qquad\overline{\overline{Z}}^{(y)}(\omega,0)=\overline{\overline{Z}}^{(y)}(\omega,2\omega_{0}-\omega)=0\qquad\forall\omega.$
(4.52)

These relations are verified making use of Eq. (4.50). For instance,

$\displaystyle\overline{\overline{Z}}^{(x)}(2\omega_{0}-\lambda,\lambda)$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{2\omega_{0}-\lambda}d\omega^{\prime}\,\int\limits_{0}^{\omega^{\prime}}d\omega\,Z^{(x)}(\omega,\lambda)\,\theta(\omega<2\omega_{0}-\lambda)=\int\limits_{0}^{2\omega_{0}-\lambda}d\omega\,Z^{(x)}(\omega,\lambda)\int
d\omega^{\prime}\,\theta(0<\omega<\omega^{\prime}<2\omega_{0}-\lambda)$ (4.53)
$\displaystyle=$
$\displaystyle\int\limits_{0}^{2\omega_{0}-\lambda}d\omega\,Z^{(x)}(\omega,\lambda)(2\omega_{0}-\lambda-\omega)=0.$

The analogous relations are valid for $W$.

For calculating the contribution of (4.2) 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 $\omega$ for the structure
containing $1/xp$ and in $\lambda$ for the structure containing $1/yp$. The
conditions (4.2) and (4.2) lead to the absence of the surface terms when
performing the parts integrations. So, we can rewrite (4.2) in the convenient
form not containing the $1/xp$ and $1/yp$ factors:

$\displaystyle\langle|\bar{s}(y)G_{\nu\alpha}(x)\Gamma\,b(0)|\bar{B}_{s}(p)\rangle=\frac{f_{B}M_{B}^{3}}{4}\int
D(\omega,\lambda)\;e^{-i\lambda\,(yp)-i\omega\,(xp)}{\rm
Sp}\Bigg{\\{}\gamma_{5}\Gamma\,(1+\not{v})\bigg{[}$
$\displaystyle\times(p_{\nu}\gamma_{\alpha}-p_{\alpha}\gamma_{\nu})\frac{1}{M_{B}}[\Psi_{A}-\Psi_{V}]-i\sigma_{\nu\alpha}\Psi_{V}$
$\displaystyle-i(x_{\nu}p_{\alpha}-x_{\alpha}p_{\nu})\left(\overline{X}_{A}^{(x)}+i\not{x}M_{B}\overline{\overline{W}}^{(x)}\right)+i(x_{\nu}\gamma_{\alpha}-x_{\alpha}\gamma_{\nu})M_{B}\left(\overline{Y}_{A}^{(x)}+\overline{W}^{(x)}+i\not{x}M_{B}\overline{\overline{Z}}^{(x)}\right)$
$\displaystyle-i(y_{\nu}p_{\alpha}-y_{\alpha}p_{\nu})\left(\overline{X}_{A}^{(y)}+i\not{y}M_{B}\overline{\overline{W}}^{(y)}\right)+i(y_{\nu}\gamma_{\alpha}-y_{\alpha}\gamma_{\nu})M_{B}\left(\overline{Y}_{A}^{(y)}+\overline{W}^{(y)}+i\not{y}M_{B}\overline{\overline{Z}}^{(y)}\right)$
$\displaystyle+\epsilon_{\nu\alpha\mu\beta}\,x^{\mu}p^{\beta}\gamma^{5}\overline{\tilde{X}}_{A}^{(x)}-\epsilon_{\nu\alpha\mu\beta}\,x^{\mu}\gamma^{\beta}\gamma^{5}M_{B}\overline{\tilde{Y}}_{A}^{(x)}+\epsilon_{\nu\alpha\mu\beta}\,y^{\mu}p^{\beta}\gamma^{5}\overline{\tilde{X}}_{A}^{(y)}-\epsilon_{\nu\alpha\mu\beta}\,y^{\mu}\gamma^{\beta}\gamma^{5}M_{B}\overline{\tilde{Y}}_{A}^{(y)}\bigg{]}\Bigg{\\}}.$
(4.54)

We emphasize that if the conditions (4.46) and (4.50) are not fulfilled, then
the parametrizations (4.2) and (4.2) are not equivalent: they differ by a
nonzero surface term. Noteworthy, the requirements of the absence of
singularities at $xp\to 0$ and $yp\to 0$ and the continuity of the 3BS (4.2)
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.2) and (4.2).

###  4.3 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.34)–(4.41) do not satisfy the constraints (4.2) and (4.2). As already
emphasized, the parametrizations (4.2) and (4.2) are then not equivalent and
cannot be both correct. So, there are two different possibilities to handle
this situation:

####  4.3.1 Scenario I

Since the behaviour of the DAs of the definite twist at small $\omega$ and
$\lambda$ are fixed, we choose the following procedure:
1\. The functions $\Psi_{A}(\omega,\lambda)$ and $\Psi_{V}(\omega,\lambda)$
remain intact.
2\. We split each function $\Phi$
$(\Phi=X_{A},Y_{A},\tilde{X}_{A},\tilde{Y}_{A},Z,W)$ into
$\Phi^{(x)}(\omega,\lambda)$ and $\Phi^{(y)}(\omega,\lambda)$ as follows:333In
principle, one can take different $w_{x}$ for each of the 3DAs. We will not
discuss this possibility but will assume that $w_{x}$ is the same for all
3DAs.

$\displaystyle\Phi^{(x)}(\omega,\lambda)=w_{x}\Phi(\omega,\lambda),\quad\Phi^{(y)}(\omega,\lambda)=w_{y}\Phi(\omega,\lambda),\quad
w_{x}+w_{y}=1.$ (4.55)

3\. We make use of the DAs of LD model for the functions
$\Phi(\omega,\lambda)$ and calculate primitives of higher orders. Taking into
account that higher primitives obtained in this way do not satisfy the
constraints (4.2) and (4.2), 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.2) with 6 modified functions
$(X_{A},Y_{A},\tilde{X}_{A},\tilde{Y}_{A},Z,W)$ which in turn also depend on
the additional parameter $a$. The simplest consistent scheme that fulfills
this requirement is the following:

$\bullet$ For the functions in Eq. (4.2) containing factors $1/xp$ and $1/yp$
(i.e. $X=X_{A},Y_{A},\widetilde{X}_{A},\widetilde{Y}_{A}$) the primitives
entering Eq. (4.2) are constructed as follows:

$\displaystyle\overline{X}^{(x)}=\frac{\partial}{\partial\lambda}\bigg{[}R(\omega,\lambda,a)\int\limits_{0}^{\omega}d\omega^{\prime}\int\limits_{0}^{\lambda}d\lambda^{\prime}X(\omega^{\prime},\lambda^{\prime})\bigg{]},$
$\displaystyle\overline{X}^{(y)}=\frac{\partial}{\partial\omega}\bigg{[}R(\omega,\lambda,a)\int\limits_{0}^{\omega}d\omega^{\prime}\int\limits_{0}^{\lambda}d\lambda^{\prime}X(\omega^{\prime},\lambda^{\prime})\bigg{]}.$
(4.56)

$\bullet$ For the functions in Eq. (4.2) 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.2) are constructed in a different way (the same formulas for $W$):

$\displaystyle\overline{\overline{Z}}^{(x)}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}}{\partial\lambda^{2}}\bigg{[}R(\omega,\lambda,a)\int\limits_{0}^{\omega}d\omega^{\prime}\int\limits_{0}^{\lambda}d\lambda^{\prime}\int\limits_{0}^{\omega^{\prime}}d\omega^{\prime\prime}\int\limits_{0}^{\lambda^{\prime}}d\lambda^{\prime\prime}Z(\omega^{\prime\prime},\lambda^{\prime\prime})\bigg{]},$
$\displaystyle\overline{\overline{Z}}^{(y)}$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}}{\partial\omega^{2}}\bigg{[}R(\omega,\lambda,a)\int\limits_{0}^{\omega}d\omega^{\prime}\int\limits_{0}^{\lambda}d\lambda^{\prime}\int\limits_{0}^{\omega^{\prime}}d\omega^{\prime\prime}\int\limits_{0}^{\lambda^{\prime}}d\lambda^{\prime\prime}Z(\omega^{\prime\prime},\lambda^{\prime\prime})\bigg{]},$
$\displaystyle\overline{Z}^{(x)}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\omega}\frac{\partial^{2}}{\partial\lambda^{2}}\bigg{[}R(\omega,\lambda,a)\int\limits_{0}^{\omega}d\omega^{\prime}\int\limits_{0}^{\lambda}d\lambda^{\prime}\int\limits_{0}^{\omega^{\prime}}d\omega^{\prime\prime}\int\limits_{0}^{\lambda^{\prime}}d\lambda^{\prime\prime}Z(\omega^{\prime\prime},\lambda^{\prime\prime})\bigg{]},$
$\displaystyle\overline{Z}^{(y)}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\lambda}\frac{\partial^{2}}{\partial\omega^{2}}\bigg{[}R(\omega,\lambda,a)\int\limits_{0}^{\omega}d\omega^{\prime}\int\limits_{0}^{\lambda}d\lambda^{\prime}\int\limits_{0}^{\omega^{\prime}}d\omega^{\prime\prime}\int\limits_{0}^{\lambda^{\prime}}d\lambda^{\prime\prime}Z(\omega^{\prime\prime},\lambda^{\prime\prime})\bigg{]}.$
(4.57)

The function $R(\omega,\lambda,a)$ is chosen in the form

$\displaystyle
R(\omega,\lambda,a)=\frac{2}{1+\exp\left(\frac{a}{2\omega_{0}-\omega-\lambda}\right)},$
(4.58)

such that for $a\neq 0$ the function itself and all its derivatives vanish at
the boundary $\omega+\lambda=2\omega_{0}$. Respectively, for $a\neq 0$, all
modified higher primitives of the necessary order vanish at the boundary
$2\omega_{0}-\omega-\lambda=0$, providing the absence of the $1/xp$ and $1/yp$
singularities at $xp\to 0$ and $yp\to 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 $\omega$
and $\lambda$ but (strongly) deviate from them near the upper boundary.

####  4.3.2 Scenario II

One declares Eq. (4.2) as the starting point but calculates the necessary
primitives using the 3DAs from Eq. (4.2). Then, however, Eq. (4.2) itself is
not complete and should contain nonzero surface terms obtained by rewriting
Eq. (4.2) to the form (4.2). These surface terms guarantee the absence of
singularities in (4.2) at $xp\to 0$ and $yp\to 0$. This scenario does not seem
to us logically satisfactory: For instance, the double-collinear limit
wang2022 ; m2023  that may be readily taken in the 3BS given by Eq. (4.2) is
not reproduced by Eq. (4.2) 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\to 0$ do not reproduce the form factors obtained
using Scenario II: the difference amounts to certain surface terms that arise
in the limit $a\to 0$.

##  5 Results for the $B_{s}\to\gamma\gamma$ amplitude and $\delta
C_{7\gamma}$

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\gamma}$.

###  5.1 The charming-loop form factors $H_{A,V}(0,0)$

Using Eqs. (2.15) and (LABEL:Acharm3b) and calculating the trace in (4.2) for

$\displaystyle\Gamma=\gamma^{\eta}(\not{k}+m_{s})\gamma^{\mu}(1-\gamma_{5}),$
(5.59)

leads to the $B_{s}\to\gamma\gamma$ amplitude expressed via the DAs and their
primitives.

The expression (4.2) contains powers of $x$ and/or $y$, but these are easy to
handle. Let us go back to Eq. (LABEL:Acharm3b): Any factor $x_{\alpha}$ may be
represented as
$x_{\alpha}\to\frac{\partial}{\partial\kappa_{\alpha}}e^{-i\kappa x}$ then
taking the $\kappa$-derivative over to
$\overline{\Gamma}_{cc}^{\mu\nu\rho\alpha}(\kappa,q)$ by parts integration in
$\kappa$. Any factor $y_{\alpha}$ may be represented as
$y_{\alpha}\to\frac{\partial}{\partial k_{\alpha}}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 (LABEL:Acharm3b) then lead to the $\delta$ functions:

$\displaystyle\int dx\rightarrow(2\pi)^{4}\delta\left(\kappa+\omega
p\right),\qquad\int dy\rightarrow(2\pi)^{4}\delta\left(k-q^{\prime}+\lambda
p\right),$ (5.60)

which allow us to further take the integrals over $\kappa$ and $k$. In the
end, the invariant functions $H_{i}$, $i={A,V}$ are given by integral
representations of the form

$\displaystyle H_{i}(q^{2},q^{\prime
2})=\int\limits_{0}^{2\omega_{0}}d\omega\int\limits_{0}^{2\omega_{0}-\omega}d\lambda\int\limits_{0}^{1}d\xi\int\limits_{0}^{1-\xi}d\eta\;h_{i}(\omega,\lambda,\xi,\eta\,|\,q^{2},q^{\prime
2})$ (5.61)

where $h_{i}$ are linear combinations of the DAs and their primitives entering
Eq. (4.2). Eq. (5.61) 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.2). We make use of the modified LD model described above.

|
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
Fig. 3:  The contributions $R_{iE}$ and $R_{iH}$ to the form factors $H_{i}$
[$i=A,V$] defined according to Eq. (5.62) for different scenarios of treating
the $B$-meson 3DA. (a,b): The appropriate modifications of the 3DAs $X_{A}$,
$Y_{A}$, etc. are taken into account (Scenario I). (c,d): The form factors are
calculated using Eq. (4.2) with the primitives obtained from $X_{A}$, $Y_{A}$,
etc. without modifications (Scenario II). (e,f): Only the contributions of
$\Psi_{A}$ and $\Psi_{V}$ are taken into account. Dashed lines show the
calculation results for $\lambda_{B_{s}}=0.45$ GeV, $a=0.05$,
$w_{x}=w_{y}=0.5$ and other inputs from Table 1, and solid lines are the fits
using Eq. (5.63).

Obviously, the form factors $H_{A}$ and $H_{V}$ depend linearly on
$\lambda_{E}^{2}$ and $\lambda_{H}^{2}$. So, we present the results for the
form factors $R_{iE}$ and $R_{iH}$ defined as follows:

$\displaystyle H_{i}(0,q^{\prime 2})$ $\displaystyle=$ $\displaystyle
R_{iE}(0,q^{\prime 2})\lambda_{E}^{2}+R_{iH}(0,q^{\prime
2})\lambda_{H}^{2},\qquad i=A,V.$ (5.62)

The form factors $R_{iE,iH}$ depend on $\lambda_{B_{s}}$, but do not contain
dependence on $\lambda_{H}$ and $\lambda_{E}$. These three parameters,
however, are expected to be strongly correlated and we take into account this
correlation later when evaluating $\delta C_{7\gamma}$.

All numerical inputs necessary for calculating $R_{i}$ are summarized in Table
1.

Table 1: Input parameters for the calculation of the form factors $R_{iE,iH}$ defined in Eq. (5.62).  $\overline{m}_{c}(\overline{m}_{c})$ | $\overline{m}_{s}(2\,{\rm GeV})$ | $M_{B_{s}}$ | $M_{\phi}$ | $f_{B_{s}}$ | $\lambda_{B_{s}}(1\,{\rm GeV})$ | $a$ | $w_{x}$
---|---|---|---|---|---|---|---
$1.30$ GeV | 0.1 GeV | 5.3 GeV | 1.020 GeV | 0.23 GeV |  $0.45\pm 0.15$ GeV | $0.05\pm 0.05$ | $0.5\pm 0.5$

Making use of the estimate for the ratio $\lambda_{B_{s}}/\lambda_{B}=1.19\pm
0.04$ khodjamirian2020  and available theoretical results for $\lambda_{B}$
kou ; braun2004 ; beneke2011 ; zwicky2021 ; im2022 , we consider the range
$\lambda_{B_{s}}=0.45\pm 0.15$ GeV with $\lambda^{0}_{B_{s}}=0.45$ GeV as the
benchmark value.

Presently, the theoretical knowledge of noncollinear 3DAs is poor to give good
arguments for the choice of the function $R$, Eq. (4.58), 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 not more than 10% in the region of “small” $\lambda$ and
$\omega$, e.g., at $\lambda/(2\omega_{0}),\omega/(2\omega_{0})\leq 0.1$. This
leads to the restriction $a\leq 0.1$. So we set the benchmark point $a=0.05$
and allow its variations in the range $0\leq a\leq 0.1$.

Other parameters in Table 1 are taken from pdg  and lattice .

We are interested in obtaining the form factors $R_{iE,iH}$ at $q^{2}=0$ and
$q^{\prime 2}=0$. Whereas $q^{2}=0$ may be readily set in the integral
representation (5.61), this integral representation is not expected to give
reliable predictions at $q^{\prime 2}=0$: The $s$-quark propagator at
$q^{\prime 2}=0$ is not sufficiently far off-shell and at $q^{\prime 2}\to 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^{\prime 2}=M_{\phi}^{2}$. To avoid this
problem, we choose the following strategy: we take the results of our
calculation for $R_{iH,iE}(0,q^{\prime 2})$ in the interval $-5<q^{\prime
2}({\rm GeV}^{2})<-0.7$ and extrapolate them numerically to $q^{\prime 2}=0$
making use of a modified pole formula which takes into account the presence of
the physical pole at $q^{\prime 2}=M_{\phi}^{2}$:

$F(q^{2})=\frac{F(0)}{\left(1-q^{2}/M_{R}^{2}\right)\left(1-\sigma_{1}(q^{2}/M_{R}^{2})-\sigma_{2}(q^{2}/M_{R}^{2})^{2}\right)}.$
(5.63)

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.63). 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 $\psi_{A}$ and $\psi_{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 II 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^{\prime 2})$ at $q^{\prime 2}=-2$ GeV2 for Scenarios I
and II.

Fig. 4 demonstrates the sensitivity of $R_{iE,iH}(0,0)$ to $\lambda_{B_{s}}$.
Since the form factors at $q^{\prime 2}=0$ are obtained via extrapolation, we
take $\lambda_{B_{s}}$ from the range $0.3<\lambda_{B_{s}}({\rm GeV})<0.5$,
calculate $R_{iE,iH}(0,q^{\prime 2})$ in the range $-5<q^{\prime 2}({\rm
GeV}^{2})<-0.7$ and then extrapolate to $q^{\prime 2}=0$ for each value of
$\lambda_{B_{s}}$.

Fig. 4:  The dependence on parameter $\lambda_{B_{s}}$ of the form factors
obtained in Scenario I: (a) $R_{iE,iH}(0,0)$, $i=A,V$; (b) Linear combination
$R_{iE}(0,0)+2R_{iH}(0,0)$ that determines $\delta C_{7\gamma}$ taking into
account approximate relation $\lambda_{H}^{2}\simeq 2\lambda_{E}^{2}$.

The same procedure applies to the dependence on parameter $a$ [enters the
function $R$, Eq. (4.58)] shown in Fig. 5. We restirct $a$ in the range $0\leq
a\leq 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.

Fig. 5:  The dependence on parameter $a$ of the form factors obtained in
Scenario I: (a) $R_{iE,iH}(0,0)$, $i=A,V$; (b) Linear combination
$R_{iE}(0,0)+2R_{iH}(0,0)$ that determines $\delta C_{7\gamma}$ taking into
account approximate relation $\lambda_{H}^{2}\simeq 2\lambda_{E}^{2}$.

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
$\omega$ and $\lambda$. 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.34)-(4.41) [Model IIB from braun2017
] and comparing with a different model [Model IIA given by Eq. (5.23) from
braun2017 ]. 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.58). Notice that we do not perform
any averaging over $\lambda_{B_{s}}$; we keep the full dependence on this
parameter in the form factors and in our final result for the correction
$\delta C_{7\gamma}$.

  * (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^{\prime 2}=0$. We make use of the calculations of the form factors in the
interval $-5\leq q^{\prime 2}({\rm GeV}^{2})\leq-0.7$. However, the value at
$q^{\prime 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\div-0.5)$ GeV2 leads to the variations of the
extrapolated value of the form factor at $q^{\prime 2}=0$ by $\pm 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^{\prime 2}=0$ for two different
sets of 3DAs and for two different Scenarios to treat the 3DAs.

Table 2: Our results for the form factors $R_{iE,iH}$ in GeV-1 defined in Eq. (5.62). In addition to the results obtained in Scenario I and Scenario II for the basic set of 3DAs given in (4.34)–(4.41), referred to as (B), we present the results obtained using Scenario I for 3DAs of Model IIA of Eq. (5.23) from braun2017 , referred to as (A).  | $10^{5}R_{VE}$ | $10^{5}R_{VH}$ | $10^{5}(R_{VE}+2R_{VH})$ | $10^{5}R_{AE}$ | $10^{5}R_{AH}$ | $10^{5}(R_{AE}+2R_{AH})$
---|---|---|---|---|---|---
Scenario I, 3DA (B) | 138.1 | 128.1 | 394.3 | 134.7 | 131.0 | 396.7
Scenario I, 3DA (A) | 146.4 | 133.7 | 413.6 | 142.7 | 137.0 | 425.4
Scenario II, 3DA (B) | 180.2 | 116.7 | 413.8 | 182.6 | 121.4 | 416.7

The individual form factors $R_{iE}$ and $R_{iH}$ demonstrate a sizeabe
dependence on the Scenario. However, in the combinations appropriate for the
calculation of $\delta C_{7\gamma}$, $R_{VE}+2R_{VH}$ and $R_{AE}+2R_{AH}$
these uncertainties cancel to great extent and we obtain for the benchmark
value $\lambda_{B_{s}}=0.45$ GeV:

$\displaystyle 10^{5}(R_{VE}+2R_{VH})\,[{\rm GeV}^{-1}]$ $\displaystyle=$
$\displaystyle 410\pm 10\,[{\rm 3DA}]\pm 15\,[a]\pm 20\,[{\rm extr}]\to 410\pm
30,$ (5.64) $\displaystyle 10^{5}(R_{AE}+2R_{AH})\,[{\rm GeV}^{-1}]$
$\displaystyle=$ $\displaystyle 415\pm 15\,[{\rm 3DA}]\pm 15\,[a]\pm 20\,[{\rm
extr}]\to 415\pm 30,$ (5.65)

and we have slightly increased our final “theoretical uncertainty”. In the
end, for a given value of the $\lambda_{B_{s}}$, we predict the combination of
the form factors (5.64) and (5.65) relevant for the calculation of
$C_{7\gamma}$ (see Subsection C) with about 8% accuracy.

###  5.2 The penguin form factor $F_{T}(0,0)$

The form factors $F_{TA,TV}(q^{2},q^{\prime 2})$ may be calculated via the
$B$-meson 2DAs using HQET formula (see e.g. offen2007 )

$\displaystyle\langle 0|\bar{s}(x)\Gamma
b(0)|\bar{B}_{s}(p)\rangle=-\frac{if_{B}M_{B}}{4}\int d\lambda e^{-i\lambda
px}\,{\rm
Sp}\bigg{\\{}\gamma_{5}\Gamma(1+\not{v})\bigg{[}\phi_{+}(\lambda)-\frac{\phi_{+}(\lambda)-\phi_{-}(\lambda)}{2vx}\not{x}\bigg{]}\bigg{\\}}$
(5.66)

leading to the following expression for the
$F_{T}(0,0)=F_{TA}(0,0)=F_{TV}(0,0)$:

$\displaystyle F_{T}(0,0)=-Q_{s}f_{B}M_{B}\int
d\lambda\frac{\phi_{+}(\lambda)(1-\lambda)+\bar{\Phi}(\lambda)}{m_{s}^{2}+\lambda(1-\lambda)M_{B}^{2}}$
(5.67)

where

$\displaystyle\bar{\Phi}(\lambda)=\int\limits_{0}^{\lambda}d\lambda^{\prime}\big{[}\phi_{+}(\lambda^{\prime})-\phi_{-}(\lambda^{\prime})\big{]}.$
(5.68)

The absence of the kinematical singularity at $vx\to 0$ in Eq. (5.66) requires
that the primitive $\bar{\Phi}(\lambda)$ vanishes at the boundaries of the 2DA
support region. Then the parts integration in $\lambda$, 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.67) 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 mnk2018 ,
corrections due to the photon emission from the $b$-quark amount to 10-20% of
the leading contribution (5.67).

For the 2DAs we use the same model from braun2017  as we do for 3DAs:

$\displaystyle\phi_{+}(\lambda)$ $\displaystyle=$
$\displaystyle\frac{5}{8\omega_{0}^{5}}\lambda(2\omega_{0}-\lambda)^{3}\theta(2\omega_{0}-\lambda),$
(5.69) $\displaystyle\phi_{-}(\lambda)$ $\displaystyle=$
$\displaystyle\frac{5}{192\omega_{0}^{5}}(2\omega_{0}-\lambda)^{2}\bigg{[}6(2\omega_{0}-\lambda)^{2}-\frac{7(\lambda_{E}^{2}-\lambda_{H}^{2})}{\omega_{0}^{2}}(15\lambda^{2}-20\omega_{0}\lambda+4\omega_{0}^{2})\bigg{]}\theta(2\omega_{0}-\lambda).$
(5.70)

An explicit check shows that
$\bar{\Phi}(\lambda=0)=\bar{\Phi}(2\omega_{0})=0$. The contribution of the
term $\sim(\lambda_{E}^{2}-\lambda_{H}^{2})$ to the form factor turns out
numerically negligible and may be safely omitted.

Fig. 6:  (a) The form factor $F_{TA}(0,q^{\prime 2})$, where $q^{\prime}$ is
the momentum of the photon emitted from the $s$-quark, for 2DAs given by
(5.69) and (5.70) and our benchmark value $\lambda_{B_{s}}=0.45$ GeV. Dotted
line - calculation results; solid line - interpolation using the modified pole
formula Eq. (5.63). (b) The dependence of $F_{TA}(0,0)$ on $\lambda_{B_{s}}$
in the range $\lambda_{B_{s}}=0.45\pm 0.15$ GeV.

Let us obtain numerical estimates for $F_{T}(0,0)$. Notice that the numerator
of the integrand in Eq. (5.67) contains factor $\lambda$, so no singularity at
$\lambda\to 0$ arises in the integrand even in the limit $m_{s}\to 0$.
Therefore, the $F_{TA}(0,0)$ may be calcualted by applying directly the
representation for the form factor Eq. (5.67). [Recall that in order to obtain
the form factors $H_{A,V}(0,0)$ we had to perform extrapolation from the
region $q^{\prime 2}\lesssim-0.5$ GeV2.] Figure 6(a) shows that the $q^{\prime
2}$-behaviour of $F_{TA}(0,q^{\prime 2})$ is well compatible with the location
of the physical pole at $q^{\prime 2}=M_{\phi}^{2}$ in a broad range of
$q^{\prime 2}<0$, up to $q^{\prime 2}=0$. Figure 6(b) illustrates the
sensitivity of $F_{T}(0,0)$ to $\lambda_{B_{s}}$. This dependence is not
negligible and partly compensates the $\lambda_{B_{s}}$-dependence of
$H_{A,V}(0,0)$, leading to more stable predictions for $\delta C_{7\gamma}$.

For our further numerical estimates we use $F_{T}(0,0)=0.155\pm 0.015$ for
$\lambda_{B_{s}}=0.45$ GeV, where the uncertainty of $\sim$10% is assigned on
the basis of the size of the neglected $1/M_{B}$-suppressed contributions
which have been calculated in mnk2018 .

###  5.3 $\delta C_{7\gamma}$

We are ready to evaluate the relative contribution of the nonfactorizable
charming loops given by Eq. (2.19):

$\displaystyle\delta_{i}C_{7\gamma}=8\pi^{2}\,Q_{s}Q_{c}\frac{C_{2}}{C_{7\gamma}m_{b}}\rho^{(i)}_{cc},\qquad\rho^{(i)}_{cc}=\frac{R_{iE}\lambda_{E}^{2}+R_{iH}\lambda_{H}^{2}}{F_{T}},\qquad
i=A,V.$ (5.71)

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 $\rho^{(i)}_{cc}$. The parameters $\lambda_{B_{s}}$,
$\lambda_{E}$ and $\lambda_{H}$, entering $\rho^{(i)}_{cc}$, are strongly
correlated with each other and should not be treated as independent
quantities. As noticed in braun2017  (see also grozin1997 ; nishikawa2014 ),
QCD sum rules suggest an approximate relation
$\lambda_{H}^{2}=2\lambda_{E}^{2}$. Combining this relation with the
constraints from the QCD equations of motion braun2017 , one obtains
approximate relations

$\displaystyle\lambda_{E}^{2}\simeq
0.6\lambda_{B_{s}}^{2},\qquad\lambda_{H}^{2}\simeq 1.2\lambda_{B_{s}}^{2}.$
(5.72)

As the result, $\rho^{i}_{cc}$ turns out to be the function of one variable,
$\lambda_{B_{s}}$, and the appropriate combinations of the form factors that
determines the correction $\delta C_{7\gamma}$ are $R_{VE}+2R_{VH}$ and
$R_{AE}+2R_{AH}$, presented in Table 2.

Fig. 7:  The functions $\rho^{(i)}_{cc},i=A,V$ [Eq. (5.71)] vs
$\lambda_{B_{s}}$.

Fig. 7 shows the dependence of $\rho^{i}_{cc}$ on $\lambda_{B_{s}}$ 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 $\rho^{(A)}_{cc}=\rho^{(V)}_{cc}$.

Using the values $C_{2}=-1.1$, $C_{7\gamma}=-0.31$ and taking
$\overline{m}_{b}(\overline{m}_{b})=4.2$ GeV for $m_{b}$, Eq. (5.71) yields:

$\displaystyle\delta C_{7\gamma}(\lambda_{B_{s}})=(0.15\,{\rm
GeV^{-1}})\,10^{2}\rho_{cc}(\lambda_{B_{s}}).$ (5.73)

The final result for $\delta C_{7\gamma}$ is very sensitive to the precise
value of $\lambda_{B_{s}}$. Recall, however, that for a given value of
$\lambda_{B_{s}}$, $\rho_{cc}$ may be calculated with an accuracy of better
than 10%. So we prefer to present our results for $\delta C_{7\gamma}$ as the
function of $\lambda_{B_{s}}$. For our benchmark point
$\lambda_{B_{s}}^{0}=0.45$ GeV, we find

$\displaystyle\delta C_{7\gamma}(\lambda^{0}_{B_{s}})=0.045\pm 0.004.$ (5.74)

For $\lambda_{B_{s}}$ in the range $0.3<\lambda_{B_{s}}({\rm GeV})<0.6$, the
corresponding $\delta C_{7\gamma}$ covers the range

$\displaystyle\delta C_{7\gamma}=(2\div 10)\%.$ (5.75)

We therefore conclude that the effect of nonfactorizable charming loops is
expected at the level of a few percent. As soon as the parameter
$\lambda_{B_{s}}$ is known with a better accuracy, our results allow one to
obtain a (relatively) accurate estimate for $\delta C_{7\gamma}$.

##  6  Discussion and conclusions

We presented a detailed analysis of NF charming loops in FCNC
$B_{s}\to\gamma\gamma$ decay and reported the following results:

(i) We derived and made use of the expression for the $\langle VVA\rangle$
charm-quark loop that is fully given in terms of the gluon field strength
$G_{\mu\nu}(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, $\bar{s}(y)G_{\mu\nu}(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_{\mu}A_{\mu}(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 m2023  and
implemented proper modifications of the corresponding 3DAs
$X_{i}(\omega,\lambda)$ at large values of their arguments.

(iii) We calculated the form factors $H_{i}(q^{2}=0,q^{\prime 2}=0)$, $i=A,V$,
describing the $B\to\gamma\gamma$ 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^{\prime 2}=0$ [$q^{\prime}$ the
momentum emitted by the light $s$-quark] was reached by an extrapolation from
the spacelike region $-5\leq q^{\prime 2}({\rm GeV}^{2})\leq-(0.5\div 1)$. An
explicit dependence of the form factors and the correction $\delta
C_{7\gamma}$ on the parameter $\lambda_{B_{s}}$ was calculated. Taking into
account all uncertainties (excluding that of $\lambda_{B_{s}}$), we found that
for any specific value of $\lambda_{B_{s}}$, $\delta C_{7\gamma}$ may be
obtained with better than 8% accuracy. For our benchmark point
$\lambda^{0}_{B_{s}}=0.45$ GeV, we found

$\displaystyle\delta C_{7\gamma}(\lambda^{0}_{B_{s}})=0.045\pm 0.004.$

For $\lambda_{B_{s}}$ in the range $0.3<\lambda_{B_{s}}({\rm GeV})<0.6$,
$\delta C_{7\gamma}$ covers the range 2-10%. Thus, one should expect the NF
charm-loop correction to $B\to\gamma\gamma$ decays at the level of a few %. As
soon as a more accurate value of $\lambda_{B_{s}}$ is available, our results
allow one to obtain the corresponding $\delta C_{7\gamma}$.

(iv) In the double-collinear kinematics that dominates the NF charm-loop
contribution to FCNC $B$-decay amplitudes in the HQ limit m2023 , the leading
contribution to the amplitude is given by the convolution of the form factor
$F_{0}$ describing the charm-quark loop Eq. (3.22) and the following
combination of the 3DAs wang2022 ; wang2022b :

$\displaystyle\Psi_{A}+\Psi_{V}+2(W^{(y)}+Y_{A}^{(y)}-\tilde{Y}_{A}^{(y)})\sim(\lambda_{E}^{2}+\lambda_{H}^{2}).$
(6.76)

[The explicit expressions for these 3DAs show that this combination is
proportional to $(\lambda_{E}^{2}+\lambda_{H}^{2})$.] Eq. (6.76) implies that
in the HQ limit, $R_{VE}=R_{VH}$ and $R_{AE}=R_{AH}$. Our results presented in
Fig. 3 show that these relations are sizeably violated by
$O(\lambda_{B_{s}}/M_{B})$ corrections that come into the game via other
Lorentz structures describing the charm-quark loop (3.22) and other 3DAs and
worth to be taken into account. Numerically, these corrections are around 20%.

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. piscopo2023 ).

Acknowledgments. We are pleased to express our gratitude to M. Ferre, E. Kou,
O. Nachtmann, and H. Sazdjian for a fruitful collaboration at the early stage
of this project, and to Q. Qin, Y.-L. Shen, C. Wang and Y.-M. Wang for
interesting discussions. Our special thanks and due to Yu-Ming Wang for his
illuminating remarks and comments. D. M. gratefully acknowledges participation
at the Erwin Schrödinger Institute (ESI) thematic program “Quantum Field
Theory at the Frontiers of the Strong Interaction” which promoted a deeper
understanding of the problems discussed and influenced much the final form of
the paper. The research was carried out within the framework of the program
“Particle Physics and Cosmology” of the National Center for Physics and
Mathematics.

##  Appendix A Comparison with the charm-loop contribution from mnk2018

Here we compare the definition of the $B\to\gamma\gamma$ amplitude (2.3) with
the definition adopted in mnk2018 . The starting point in mnk2018  is the
matrix element

$\displaystyle H_{\rho\eta}(q,q^{\prime})=i\int dze^{iqz}\langle
0|T\\{\bar{c}\gamma_{\rho}c(z),j^{\rm e.m.}_{\eta}(0)|B_{s}(p)\rangle,\qquad
p=q+q^{\prime},$ (A.77)

where quark operators are Heisenberg operators in the SM, i.e. the
corresponding $S$-matrix includes weak and strong interactions. For real
photons in the final state, at least one soft gluon should be emitted from the
charm-quark loop, we expand the $S$-matrix to first order in $G_{F}$ and first
order in $g_{s}$:

$\displaystyle H_{\rho\eta}(q^{\prime},q)=i\int dze^{iqz}\langle
0|T\\{\bar{c}(z)\gamma_{\rho}c(z),i\int dy\,L_{\text{weak}}(y),i\int
dx\,L_{Gcc}(x),eQ_{s}\,\bar{s}(0)\gamma_{\eta}s(0)\\}|B_{s}(p)\rangle.$ (A.78)

For comparison with Eq. (2.3), we place $L_{\rm weak}$ at $y=0$ by shifting
coordinates of all operators through the translation ${\cal
O}(x)=e^{i\hat{P}y}{\cal O}(x-y)e^{-iy\hat{P}}$. Using the relations $\langle
0|e^{i(\hat{P}y)}=\langle 0|$ and
$e^{-i(y\hat{P})}|B_{s}(p)\rangle=e^{-i(py)}|B_{s}(p)\rangle$, and changing
the variables $x-y\to x$, $z-y\to z$, $y\to-y$, we find

$\displaystyle H_{\rho\eta}(q,q^{\prime})=i^{3}\,e\,Q_{s}\int
dxdydz\>e^{iqz+iq^{\prime}y}\langle
0|T\\{\bar{c}(z)\gamma_{\rho}c(z),L_{\text{weak}}(0),L_{Gcc}(x),\bar{s}(y)\gamma_{\eta}s(y)\\}|B_{s}(p)\rangle.$
(A.79)

Taking into account that $L_{\rm weak}=-H^{b\to s\bar{c}c[8\times 8]}_{\rm
weak}$ [the latter given by Eq. (2.11)], we obtain

$\displaystyle H_{\rho\eta}(q,q^{\prime})$ $\displaystyle=$ $\displaystyle
2C_{2}\,\frac{G_{F}}{\sqrt{2}}V_{cb}V^{*}_{cs}e\,Q_{s}\,i\int
dxdydz\>e^{iqz+iq^{\prime}y}\langle
0|T\\{\bar{c}(z)\gamma_{\rho}c(z),\bar{c}(0)\gamma_{\mu}(1-\gamma_{5})t^{a}c(0),\bar{c}(x)\gamma_{\nu}t^{b}c(x)|0\rangle$
(A.80) $\displaystyle\hskip 113.81102pt\times\langle
0|T\left\\{\bar{s}(y)\gamma_{\eta}s(y),\bar{s}(0)\gamma_{\mu}(1-\gamma_{5})t^{a}b(0),B^{b}_{\nu}(x)\right\\}|B_{s}(p)\rangle$

It is convenient to insert under the integral (A.80) the identity

$\displaystyle B^{b}_{\nu}(x)=\frac{1}{(2\pi)^{4}}\int d\kappa
dx^{\prime}B^{b}_{\nu}(x^{\prime})e^{i\kappa(x-x^{\prime})}.$ (A.81)

This allows us to isolate the contribution of the charm-quark loop
$\Gamma^{\mu\nu\rho(ab)}_{cc}(\kappa,q)$:

$\displaystyle H_{\rho\eta}(q,q^{\prime})$ $\displaystyle=$ $\displaystyle
2C_{2}\,\frac{G_{F}}{\sqrt{2}}V_{cb}V^{*}_{cs}eQ_{s}\frac{i}{(2\pi)^{4}}\int
dy\,e^{iq^{\prime}y}\;d\kappa e^{-i\kappa x^{\prime}}\langle
0|T\\{\bar{s}(0)\gamma_{\mu}(1-\gamma_{5})b(0),B^{b}_{\nu}(x^{\prime}),\bar{s}(y)\gamma_{\eta}s(y)\\}|B_{s}(p)\rangle$
(A.82) $\displaystyle\times\underbrace{\int dxdz\,e^{iqz+i\kappa x}\langle
0|T\\{\bar{c}(z)\gamma_{\rho}c(z),\bar{c}(0)\gamma_{\mu}(1-\gamma_{5})t^{a}c(0),\bar{c}(x)\gamma_{\nu}t^{b}c(x)\\}|0\rangle}_{\Gamma^{\mu\nu\rho(ab)}_{cc}(\kappa,q)}.$

Using momentum representation for the $s$-quark propagator

$\displaystyle\langle
0|T\\{s(y)\bar{s}(0)\\}|0\rangle=\frac{1}{(2\pi)^{4}i}\int
dk\,e^{-iky}\frac{\not{k}+m_{s}}{m_{s}^{2}-k^{2}-i0},$ (A.83)

we obtain

$\displaystyle H_{\rho\eta}(q,q^{\prime})$ $\displaystyle=$ $\displaystyle
2C_{2}\,\frac{G_{F}}{\sqrt{2}}V_{cb}V^{*}_{cs}eQ_{s}$
$\displaystyle\times\underbrace{\frac{1}{(2\pi)^{8}}\int{dk}dye^{-i(k-q^{\prime})y}dxd\kappa
e^{-i\kappa x}\Gamma_{cc}^{\mu\nu\rho(ab)}(\kappa,q)\langle
0|\bar{s}(y)\gamma^{\eta}\frac{\not{k}+m_{s}}{m_{s}^{2}-k^{2}}\gamma^{\mu}(1-\gamma^{5})t^{a}B^{b}_{\nu}(x)b(0)|B_{s}(p)\rangle}_{A_{\rho\eta}}.$

Comparing this expression with Eq. (2.15) gives a useful relation

$\displaystyle A^{B\to\gamma\gamma}_{\rm
charm}=2eQ_{c}\,H_{\rho\eta}\,\varepsilon_{\rho}\varepsilon^{\prime}_{\eta}.$
(A.85)

The amplitude $H_{\rho\eta}$ may be decomposed via the form factors
$\tilde{H}_{A,V}$ (denoted as $H_{A,V}$ in Eq. (5.2) from mnk2018 ):

$\displaystyle
H_{\rho\eta}(q,q^{\prime})=-\frac{G_{F}}{\sqrt{2}}V_{cb}V^{*}_{cs}e\Big{[}\epsilon_{\rho\eta
qq^{\prime}}\tilde{H}_{V}-i\left(g_{\rho\eta}q^{\prime}q-q^{\prime}_{\rho}q_{\eta}\right)\tilde{H}_{V}\Big{]}.$
(A.86)

Comparing with (LABEL:Acharm3b) we find

$\displaystyle\tilde{H}_{V(A)}=-2C_{2}Q_{s}H_{V(A)}.$ (A.87)

##  Appendix B Anathomy of the form factors $H_{A,V}(0,q_{0}^{\prime 2})$ at
$q_{0}^{\prime 2}=-2$ GeV2.

This Appendix presents detailed numerical results for the form factors at
$q^{2}=0$ and $q^{\prime 2}=q^{\prime 2}_{0}=-2$ GeV2, giving separately
contributions of different Lorentz 3DAs for the coefficients $R_{iE}$ and
$R_{iH}$ defined in Eq. (5.62). Table 3 gives the results obtained with
Scenario I of Section 4.3 for 3DAs given by Eqs. (4.34)–(4.41) [Model IIB of
Eq. (5.27) from braun2017 ]. 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
braun2017 . 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 $\lambda_{B_{s}}=0.45$ GeV and the
central values of $m_{c}$ and $m_{s}$. For 3DAs modified according to Scenario
I, $a=0.05$.

Table 3:  $H_{V}(0,q_{0}^{\prime 2})$ and $H_{A}(0,q_{0}^{\prime 2})$ at $q_{0}^{\prime 2}=-2$ GeV2 calculated for Scenario I and 3DAs of Eqs. (4.34)–(4.41) for $\lambda_{B_{s}}=0.45$ GeV and $a=0.05$. $H_{i}(0,q^{\prime 2})=R_{iE}\lambda_{E}^{2}+R_{iH}\lambda_{H}^{2},i=A,V$. | $10^{5}H_{V}(q^{2},q_{0}^{\prime 2})$ | $10^{5}H_{A}(q^{2},q_{0}^{\prime 2})$
---|---|---
Lorentz 3DA |  $10^{5}R_{VE}$ [GeV-1] |  $10^{5}R_{VH}$ [GeV-1] |  $10^{5}R_{AE}$ [GeV-1] |  $10^{5}R_{AH}$ [GeV-1]
$\psi_{A}$ | 45.84 | 7.54 | 43.44 | 7.67
$\psi_{V}$ | 4.62 | 25.2 | 4.61 | 26.7
${X}_{A}$ | $-0.28\,w_{x}+0.49\,w_{y}$ | $-1.64\,w_{x}+0.99\,w_{y}$ | $-0.39\,w_{x}+0.24\,w_{y}$ | $-1.63\,w_{x}+0.81\,w_{y}$
${Y}_{A}$ | $-8.37\,w_{x}+1.52\,w_{y}$ | $8.33\,w_{x}+0.06\,w_{y}$ | $-8.52\,w_{x}+0.70\,w_{y}$ | $8.48\,w_{x}+0.63\,w_{y}$
${\tilde{X}}_{A}$ | $-0.65\,w_{x}-0.05\,w_{y}$ | $-0.11\,w_{x}+0.11\,w_{y}$ | $-0.70\,w_{x}+0.22\,w_{y}$ | $-0.07\,w_{x}-0.19\,w_{y}$
${\tilde{Y}}_{A}$ | $1.65\,w_{x}-0.01\,w_{y}$ | $-2.45\,w_{x}-1.51\,w_{y}$ | $1.70\,w_{x}-0.59\,w_{y}$ | $-2.51\,w_{x}-0.68\,w_{y}$
${W}$ | $-0.61\,w_{x}-0.39\,w_{y}$ | $-0.61\,w_{x}-0.39\,w_{y}$ | $-0.63\,w_{x}-0.35\,w_{y}$ | $-0.63\,w_{x}-0.35\,w_{y}$
${Z}$ | $0.93\,w_{x}+0.27\,w_{y}$ | $-0.31\,w_{x}-0.19\,w_{y}$ | $0.95\,w_{x}+0.03\,w_{y}$ | $-0.31\,w_{x}-0.02\,w_{y}$

Setting $w_{x}=\frac{1}{2}(1+\delta w)$, $w_{y}=\frac{1}{2}(1-\delta w)$ and
summing the contribuons of all 3DAs in Table 3, we obtain,

$\displaystyle H_{V}(0,q_{0}^{\prime 2})=(47.72-4.58\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{E}^{2}+(33.88+2.07\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{H}^{2}.$ (B.88)

For the expected relationship $\lambda_{H}^{2}=2\lambda_{E}^{2}$, the
dependence on $\delta w$ is further suppressed leading to extremely stable
predicitons with respect to $\delta_{w}$:

$\displaystyle H_{V}(0,q_{0}^{\prime 2})=(115.48-0.44\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{E}^{2}.$ (B.89)

Table 4: $H_{V}(0,q_{0}^{\prime 2})$ and $H_{A}(0,q_{0}^{\prime 2})$ at $q_{0}^{\prime 2}=-2$ GeV2 calculated for Scenario II and 3DAs of Eqs. (4.34)–(4.41); $\lambda_{B_{s}}=0.45$ GeV.  | $10^{5}H_{V}(q^{2},q_{0}^{\prime 2})$ | $10^{5}H_{A}(q^{2},q_{0}^{\prime 2})$
---|---|---
Lorentz 3DA |  $10^{5}R_{VE}$ [GeV-1] |  $10^{5}R_{VH}$ [GeV-1] |  $10^{5}R_{AE}$ [GeV-1] |  $10^{5}R_{AH}$ [GeV-1]
$\psi_{A}$ | 45.84 | 7.54 | 43.44 | 7.67
$\psi_{V}$ | 4.62 | 25.2 | 4.61 | 26.7
${X}_{A}$ | $-1.54\,w_{x}-1.17\,w_{y}$ | $-4.08\,w_{x}-1.71\,w_{y}$ | $-1.02\,w_{x}-1.40\,w_{y}$ | $-4.73\,w_{x}+0.11\,w_{y}$
${Y}_{A}$ | $-17.4\,w_{x}-5.93\,w_{y}$ | $25.3\,w_{x}-11.4\,w_{y}$ | $-17.6\,w_{x}+5.36\,w_{y}$ | $25.4\,w_{x}-12.0\,w_{y}$
${\tilde{X}}_{A}$ | $-1.73\,w_{x}-0.80\,w_{y}$ | $-0.12\,w_{x}+0.32\,w_{y}$ | $-1.62\,w_{x}+0.59\,w_{y}$ | $-0.42\,w_{x}-0.65\,w_{y}$
${\tilde{Y}}_{A}$ | $3.53\,w_{x}+10.9\,w_{y}$ | $-4.89\,w_{x}-5.69\,w_{y}$ | $3.51\,w_{x}+11.5\,w_{y}$ | $-4.96\,w_{x}-5.11\,w_{y}$
${W}$ | $2.27\,w_{x}-0.52\,w_{y}$ | $2.27\,w_{x}-0.52\,w_{y}$ | $2.38\,w_{x}-0.38\,w_{y}$ | $2.38\,w_{x}-0.38\,w_{y}$
${Z}$ | $6.45\,w_{x}+0.75\,w_{y}$ | $-0.96\,w_{x}-0.50\,w_{y}$ | $6.28\,w_{x}+0.05\,w_{y}$ | $-0.93\,w_{x}-0.02\,w_{y}$

For Scenario II ( Table 4), one finds

$\displaystyle H_{V}(0,q_{0}^{\prime 2})=(53.81-11.76\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{E}^{2}+(31.75+18.51\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{H}^{2}.$ (B.90)

For $\lambda_{H}^{2}=2\lambda_{E}^{2}$, this gives

$\displaystyle H_{V}(0,q_{0}^{\prime 2})=(117.31+25.25\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{E}^{2}.$ (B.91)

For $\delta_{w}=0$, both Scenario I and II give close results but the
sensitivity of the form factor to $\delta_{w}$ in Scenario II is strong.

Table 5:  $H_{V}(0,q_{0}^{\prime 2})$ and $H_{A}(0,q_{0}^{\prime 2})$ at $q_{0}^{\prime 2}=-2$ GeV2 calculated for Scenario I and 3DAs of Model IIA given by Eq. (5.23) of braun2017  for $\lambda_{B_{s}}=0.45$ GeV and $a=0.05$.  | $10^{5}H_{V}(q^{2},q_{0}^{\prime 2})$ | $10^{5}H_{A}(q^{2},q_{0}^{\prime 2})$
---|---|---
Lorentz 3DA |  $10^{5}R_{VE}$ [GeV-1] |  $10^{5}R_{VH}$ [GeV-1] |  $10^{5}R_{AE}$ [GeV-1] |  $10^{5}R_{AH}$ [GeV-1]
$\psi_{A}$ | 46.7 | 6.97 | 44.8 | 7.05
$\psi_{V}$ | 4.74 | 28.9 | 4.74 | 30.1
${X}_{A}$ | $0.08\,w_{x}+0.00\,w_{y}$ | $-0.39\,w_{x}+0.17\,w_{y}$ | $0.07\,w_{x}-0.04\,w_{y}$ | $-0.39\,w_{x}+0.18\,w_{y}$
${Y}_{A}$ | $-3.50\,w_{x}+0.96\,w_{y}$ | $1.79\,w_{x}-0.31\,w_{y}$ | $-3.56\,w_{x}+0.70\,w_{y}$ | $1.82\,w_{x}-0.21\,w_{y}$
${\tilde{X}}_{A}$ | $-0.17\,w_{x}-0.02\,w_{y}$ | $0.04\,w_{x}+0.03\,w_{y}$ | $-0.18\,w_{x}+0.06\,w_{y}$ | $0.05\,w_{x}-0.05\,w_{y}$
$\tilde{Y}_{A}$ | $0.52\,w_{x}+0.31\,w_{y}$ | $-1.18\,w_{x}-0.94\,w_{y}$ | $0.54\,w_{x}+0.20\,w_{y}$ | $-1.21\,w_{x}-0.68\,w_{y}$
${W}$ | $-0.04\,w_{x}-0.04\,w_{y}$ | $-0.04\,w_{x}-0.04\,w_{y}$ | $-0.04\,w_{x}-0.04\,w_{y}$ | $-0.04\,w_{x}-0.04\,w_{y}$
${Z}$ | $0.13\,w_{x}+0.05\,w_{y}$ | $-0.08\,w_{x}-0.05\,w_{y}$ | $0.14\,w_{x}+0.01\,w_{y}$ | $-0.08\,w_{x}-0.01\,w_{y}$

For Scenario I but using different 3DAs, those of Model IIA of Eq. (5.23) from
braun2017 , one finds (Table 5):

$\displaystyle H_{V}(0,q_{0}^{\prime 2})=(50.58-2.12\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{E}^{2}+(35.37+0.64\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{H}^{2}.$ (B.92)

For $\lambda_{H}^{2}=2\lambda_{E}^{2}$, this yields:

$\displaystyle H_{V}(0,q_{0}^{\prime 2})=(121.32-0.84\,\delta_{w})\mbox{
[GeV${}^{-1}$]}\lambda_{E}^{2}.$ (B.93)

## References

  * (1) M. Beneke, G. Buchalla, M. Neubert, and C. T. Sachrajda, Penguins with Charm and Quark-Hadron Duality, Eur. Phys. J. C61, 439 (2009).
  * (2) M. Ciuchini, M. Fedele, E. Franco, A. Paul, L. Silvestrini, Lessons from the $B^{0,+}\to K^{*0,+}\mu^{+}\mu^{-}$ angular analyses, Phys. Rev. D 103, 015030 (2021).
  * (3) M. Ciuchini, M. Fedele, E. Franco, A. Paul, L. Silvestrini, New Physics without bias: Charming Penguins and Lepton Universality Violation in $b\to s\ell^{+}\ell^{-}$ decays, Eur. Phys. J. C 83, 64 (2023).
  * (4) D. Guadagnoli, B Discrepancies Hold Their Ground, Symmetry 13, 1999 (2021).
  * (5) M. Algueró, B. Capdevila, A. Crivellin, J. Matias, Disentangling Lepton Flavour Universal and Lepton Flavour Universality Violating Effects in $b\to s\ell^{+}\ell^{-}$ Transitions, Phys. Rev. D 105, 113007 (2022).
  * (6) N. Gubernari, M. Reboud, D. van Dyk, and J. Virto, Improved theory predictions and global analysis of exclusive $b\to s\mu^{+}\mu^{-}$ processes, JHEP 09, 133 (2022).
  * (7) T. Hurth, F. Mahmoudi, D. Martinez Santos, and S. Neshatpour, Neutral current B-decay anomalies, Springer Proc. Phys. 292, 11 (2023).
  * (8) A. Greljo, J. Salko, A. Smolkovic, and P. Stangl, Rare b decays meet high-mass Drell-Yan, JHEP 05, 087 (2023).
  * (9) M. Ciuchini, M. Fedele, E. Franco, A. Paul, L. Silvestrini, Constraints on Lepton Universality Violation from Rare $B$ Decays, Phys. Rev. D 107, 055036 (2023).
  * (10) D. Guadagnoli, C. Normand, S. Simula, L. Vittorio, From $D_{s}\to\gamma$ in lattice QCD to $B_{s}\to\mu\mu\gamma$ at high $q^{2}$, JHEP 07, 112 (2023).
  * (11) D. Guadagnoli, C. Normand, S. Simula, L. Vittorio, Insights on the current semi-leptonic $B$-decay discrepancies – and how $B_{s}\to\mu^{+}\mu^{-}\gamma$ can help, e-Print: 2308.00034 [hep-ph].
  * (12) M. B. Voloshin, Large $O(m_{c}^{-2})$ nonperturbative correction to the inclusive rate of the decay $B\to X_{s}\gamma$, Phys. Lett. B397, 275 (1997).
  * (13) Z. Ligeti, L. Randall, and M. B. Wise, Comment on nonperturbative effects in $\bar{B}\to X_{s}\gamma$, Phys. Lett. B402, 178 (1997).
  * (14) G. Buchalla, G. Isidori, and S. J. Rey, Corrections of order $\Lambda_{\rm QCD}^{2}/m_{c}^{2}$ to inclusive rare $B$ decays, Nucl. Phys. B511, 594 (1998).
  * (15) A. Khodjamirian, R. Ruckl, G. Stoll, and D. Wyler, QCD estimate of the long distance effect in $B\to K^{*}\gamma$, Phys. Lett. B402, 167 (1997).
  * (16) P. Ball and R. Zwicky, Time-dependent CP Asymmetry in $B\to K^{*}\gamma$ as a (Quasi) Null Test of the Standard Model, Phys. Lett. B642, 478 (2006).
  * (17) P. Ball, G. W. Jones, and R. Zwicky, $B\to V\gamma$ beyond QCD factorisation, Phys. Rev. D75, 054004 (2007).
  * (18) I. I. Balitsky, V. M. Braun, A. V. Kolesnichenko, Radiative Decay $\Sigma^{+}\to p\gamma$ in Quantum Chromodynamics, Nucl. Phys. B312, 509 (1989).
  * (19) P. Ball and V. Braun, Higher twist distribution amplitudes of vector mesons in QCD: Twist - 4 distributions and meson mass corrections, Nucl. Phys. B543, 201 (1999).
  * (20) P. Ball, Theoretical update of pseudoscalar meson distribution amplitudes of higher twist: The Nonsinglet case, JHEP 9901, 010 (1999).
  * (21) A. Khodjamirian, T. Mannel, A. Pivovarov, and Y.-M. Wang, Charm-loop effect in $B\to K^{(*)}l^{+}l^{-}$ and $B\to K^{*}\gamma$, JHEP 09, 089 (2010).
  * (22) H. Kawamura, J. Kodaira, C.-F. Qiao, and K. Tanaka, B-meson light cone distribution amplitudes in the heavy quark limit, Phys. Lett. B523, 111 (2001), Erratum: Phys. Lett. B536, 344 (2002).
  * (23) V. Braun, Y. Ji, and A. Manashov, Higher-twist B-meson Distribution Amplitudes in HQET”, JHEP 05, 022 (2017).
  * (24) V. M. Braun and I. Halperin, Soft contribution to the pion form-factor from light cone QCD sum rules, Phys. Lett. B328, 457 (1994).
  * (25) A. Khodjamirian, T. Mannel, and N. Offen, Form-factors from light-cone sum rules with B-meson distribution amplitudes, Phys. Rev. D75, 054013 (2007).
  * (26) N. Gubernari, D. van Dyk, J. Virto, Non-local matrix elements in $B_{(s)}\to\\{K^{(*)},\phi\\}\ell^{+}\ell^{-}$, JHEP 02, 088 (2021).
  * (27) A. Kozachuk and D. Melikhov, Revisiting nonfactorizable charm-loop effects in exclusive FCNC $B$ decays, Phys. Lett. B786, 378 (2018).
  * (28) D. Melikhov, Charming loops in exclusive rare FCNC $B$-decays, EPJ Web Conf. 222, 01007 (2019).
  * (29) D. Melikhov, Nonfactorizable charming loops in FCNC B decay versus B-decay semileptonic form factors, Phys. Rev. D106, 054022 (2022).
  * (30) D. Melikhov, Three-particle distribution in the B meson and charm-quark loops in FCNC B decays, Phys. Rev. D108, 034007 (2023).
  * (31) Q. Qin, Yue-Long Shen, Chao Wang and Yu-Ming Wang, Deciphering the long-distance penguin contribution to $\bar{B}_{d,s}\to\gamma\gamma$ decays, e-Print: 2207.02691 [hep-ph].
  * (32) A. Kozachuk, D. Melikhov, and N. Nikitin, Rare FCNC radiative leptonic $B_{s,d}\to\gamma l^{+}l^{-}$ decays in the Standard Model, Phys. Rev. D97, 053007 (2018).
  * (33) B. Grinstein, M. J. Savage and M. B. Wise, $B\to X(s)e^{+}e^{-}$ in the Six Quark Model, Nucl. Phys. B 319 271 (1989).
  * (34) A. J. Buras and M. Munz, Effective Hamiltonian for $B\to X(s)e^{+}e^{-}$ beyond leading logarithms in the NDR and HV schemes, Phys. Rev. D 52, 186 (1995) [hep-ph/9501281].
  * (35) G. Buchalla, A. J. Buras, and M. E. Lautenbacher, Weak decays beyond leading logarithms, Rev. Mod. Phys. 68, 1125 (1996).
  * (36) D. Melikhov, N. Nikitin and S. Simula, Lepton asymmetries in exclusive $b\to sl^{+}l^{-}$ decays as a test of the standard model, Phys. Lett. B 430, 332 (1998).
  * (37) D. Melikhov, N. Nikitin and S. Simula, Rare exclusive semileptonic $b\to s$ transitions in the standard model, Phys. Rev. D 57, 6814 (1998).
  * (38) D. Melikhov and N. Nikitin, Rare radiative leptonic decays $B_{d,s}\to\gamma l^{+}l^{-}$, Phys. Rev. D70, 114028 (2004).
  * (39) F. Kruger and D. Melikhov, Gauge invariance and form-factors for the decay $B\to\gamma l^{+}l^{-}$, Phys. Rev. D 67, 034002 (2003).
  * (40) W. Lucha and D. Melikhov, The puzzle of the $\pi\to\gamma\gamma^{*}$ transition form factor, J. Phys. G39, 045003 (2012).
  * (41) Max Ferre, private communication.
  * (42) Yu-Ming Wang, private communication.
  * (43) C.-D. Lü, Y.-L. Shen, Y.-M. Wang, and Y.-B. Wei, QCD calculations of $B\to\pi,K$ form factors with higher-twist corrections, JHEP 01, 024 (2019).
  * (44) A. Khodjamirian, R. Mandal, and T. Mannel, Inverse moment of the $B_{s}$-meson distribution amplitude from QCD sum rule, JHEP 10, 043 (2020).
  * (45) P. Ball and E. Kou, $B\to\gamma e\nu$ transitions from QCD sum rules on the light cone, JHEP 0304, 029 (2003).
  * (46) V. M. Braun, D. Yu. Ivanov, and G. P. Korchemsky. The B meson distribution amplitude in QCD, Phys. Rev. D69 034014 (2004).
  * (47) M. Beneke and J. Rohrwild, B meson distribution amplitude from $B\to\gamma l\nu$, Eur. Phys. J. C71, 1818 (2011).
  * (48) T. Janowski, B. Pullin, and R. Zwicky, Charged and neutral $\bar{B}_{u,d,s}\to\gamma$ form factors from light cone sum rules at NLO, JHEP 12, 008 (2021).
  * (49) M. A. Ivanov and D. Melikhov, Theoretical analysis of the leptonic decay $B\to lll^{\prime}\nu^{\prime}$, Phys. Rev. D105, 014028 (2022); Phys. Rev. D106, 119901(E) (2022).
  * (50) Particle Data Group (R. L. Workman et al), Review of Particle Physics. Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * (51) FLAG Collaboration (Y. Aoki et al.), FLAG Review 2021, Eur. Phys. J. C 82, 869 (2022).
  * (52) A. G. Grozin and M. Neubert, Asymptotics of heavy meson form-factors, Phys. Rev. D55, 272 (1997).
  * (53) T. Nishikawa and K. Tanaka, QCD Sum Rules for Quark-Gluon Three-Body Components in the B Meson, Nucl. Phys. B879, 110 (2014).
  * (54) M. L. Piscopo and A. Rusov, Non-factorisable effects in the decays $\bar{B}^{0}_{s}\to D^{+}_{s}\pi^{-}$ and $\bar{B}^{0}\to D^{+}K^{-}$ from LCSR, e-Print: 2307.07594 [hep-ph].