Precision calculations of the double radiative bottom-meson decays in soft-collinear effective theory

Employing the systematic framework of soft-collinear effective theory (SCET) we perform an improved calculation of the leading-power contributions to the double radiative $B_{d, \, s}$-meson decay amplitudes in the heavy quark expansion. We then construct the QCD factorization formulae for the subleading power contributions arising from the energetic photon radiation off the constituent light-flavour quark of the bottom meson at tree level. Furthermore, we explore the factorization properties of the subleading power correction from the effective SCET current $J^{(A2)}$ at ${\cal O} (\alpha_s^0)$. The higher-twist contributions to the $B_{d, \, s} \to \gamma \gamma $ helicity form factors from the two-particle and three-particle bottom-meson distribution amplitudes are evaluated up to the twist-six accuracy. In addition, the subleading power weak-annihilation contributions from both the current-current and QCD penguin operators are taken into account at the one-loop accuracy. We proceed to apply the operator-production-expansion-controlled dispersion relation for estimating the power-suppressed soft contributions to the double radiative $B_{d, \, s}$-meson decay form factors. Phenomenological explorations of the radiative $B_{d, \, s} \to \gamma \, \gamma$ decay observables in the presence of the neutral-meson mixing, including the CP-averaged branching fractions, the polarization fractions and the time-dependent CP asymmetries, are carried out subsequently with an emphasis on the numerical impacts of the newly computed ingredients together with the theory uncertainties from the shape parameters of the HQET bottom-meson distribution amplitudes.

with the complex short-distance functions computed explicitly at the one-loop accuracy.
Previous calculations of the double radiative B s → γγ decay amplitude with the aid of the static bottom-and strange-quark approximation [19][20][21][22][23][24][25] apparently introduce conceptual and technical limitations for the theory description of the corresponding hadronic matrix elements. In addition, the long-distance hadronic contributions from the charm-meson rescattering mechanism [26,27] and from the vector-meson coupling to the on-shell photon [28,29] have been estimated with model-dependent phenomenological approaches, which prevent us from improving the current theory predictions systematically in the long run.
Advancing our understanding of the radiative penguin B d, s → γγ decay amplitudes further by taking advantage of the perturbative factorization approach and the hadronic dispersion relation will be of wide interest for providing the helpful insights of evaluating the subleading power corrections to exclusive heavy-to-light B q -meson decays (with q = d, s) in QCD and for probing the dynamical information of the bottom-meson distribution amplitudes in HQET. The major new ingredients of the present paper at the technical level can be summarized in the following.
• The next-to-leading logarithmic (NLL) resummmation for the enhanced logarithms of m b /Λ QCD entering the leading-power factorization formula of B d,s → γγ will be accomplished by solving the renormalization group equations for the hard matching coefficients and the leading-twist B-meson LCDA at two loops. As a consequence, our result goes beyond the leading logarithmic (LL) approximation implemented in [16]. Furthermore, the complete scale-independent factorization formulae for the double radiative decay amplitudes at O(α 2 s ) will be achieved by including the two-loop b → (s, d) γ matrix elements of the QCD penguin operators.
• Applying the equations of motion for the soft quark field and the effective heavy quark, we will derive the QCD factorization formulae of the subleading power corrections to the energetic photon radiation off the (massless)-light quark at tree level on the basis of the two-particle and three-particle B q -meson distribution amplitudes. Moreover, it will be verified explicitly that the convolution integrals of the hard-collinear matching coefficients and B q -meson LCDAs appearing in the constructed factorization formulae are convergent. By contrast, the naïve soft-collinear factorization formula of the power suppressed contribution from the light-quark mass effect suffers from the rapidity divergences already at tree level.
• We will derive the perturbative factorization formula for the subleading power softcollinear effective theory (SCET) matrix element 0| ξ W c Γ [i D s /(2m b )] h v |B q by employing the well-established relations between the non-local operators at leading order (LO) in α s . In light of the canonical behaviours of the appearing B-meson distribution amplitudes at small quark and gluon momenta, the obtained convolution integrals of the leading-and higher-twist B q -meson distribution amplitudes with the short-distance functions are both convergent.
• Applying the QCD factorization approach, the higher-twist corrections to the double radiative decay form factors from both the two-particle and three-particle B q -meson distribution amplitudes will be evaluate at O(α 0 s ) up to the twist-six accuracy with the systematic parametrization of the HQET matrix element for the three-body light-ray operatorq s (τ 1 n) g s G µν (τ 2 n) Γ h v β (0) [30]. The resulting factorization formulae expressed in terms of the convolution integrals of the twist-three and twist-four LCDAs with the corresponding hard-collinear functions appear to converge for the phenomenologically acceptable models of φ − B (ω, µ) and Φ 4 (ω 1 , ω 2 , µ).
• The subleading power soft contributions to the radiative B d, s → γγ decay amplitudes characterizing the resolved photon corrections will be computed with the dispersion approach motivated from the operator-product-expansion (OPE) technique and the parton-hadron duality ansatz, which has been successfully applied to estimate the non-perturbative correction to the pion-photon transition form factor at experimentally accessible momentum transfer [31][32][33]. To this end, the spectral representations of the soft-collinear factorization formulae for the two transverse B → γ * form factors at NLO in QCD obtained in [34,35] will be in demand.
The remainder of this paper is structured as follows. We set up the computational framework in Section 2 by introducing the effective weak Hamiltonian governing the b → (s, d) γ transitions and by expressing the decay amplitude of B q → γγ to the lowest non-vanishing order in the electromagnetic interactions in terms of the appropriate QCD correlation functions, where the two equivalent bases of the Lorenz-invariant amplitudes are also discussed for the purpose of the phenomenological explorations. We proceed to determine the SCET I representations of the above-mentioned QCD correlation functions at leading power in the heavy quark expansion by integrating out the short-distance modes with virtualities of order m 2 b and then present the analytical expressions of the resulting hard functions for an arbitrary quark-mass ratio m c /m b at O(α s ) in Section 3. Subsequently, the soft-collinear factorization formulae for the yielding SCET I correlation functions will be constructed by integrating out the hard-collinear fluctuations. The radiative jet function from matching SCET I → SCET II at one loop will be also collected here. The NLL resummation improved factorization formulae for the helicity form factors will be further derived with the standard renormalization-group formalism. In Section 4 we turn to compute the subleading power corrections generated by the energetic photon emission from the light soft quark and to derive their tree-level factorization formulae in terms of the HQET B-meson distribution amplitudes. We then evaluate the subleading power matrix element of the effective heavy-to-light current ξ W c Γ [i D s /(2m b )] h v from matching QCD → SCET I at LO in α s . Drawing upon the QCD factorization technique, we further calculate the subleading power contributions from the higher-twist B-meson distribution amplitudes due to both the off light-cone corrections and the higher Fock-state effects. The power suppressed weak-annihilation contributions from all the four-quark operators at one loop and the (anti)-collinear photon radiation off the bottom quark at tree level will be presented here as well by matching QCD → SCET II directly. Section 5 is devoted to an estimate of the long-distance photon contributions to the helicity amplitudes of B d, s → γγ based upon the OPE-controlled dispersion technique at O(α s ). Phenomenological implications of the newly derived expressions for the double radiative decay amplitudes including the renormalization-group improvement for the leading power contributions and the updated analysis of the subleading power corrections will be investigated in Section 6 with an emphasis on the shape-parameter dependence of the achieved theory predictions on the leading-twist B-meson LCDA. Having at our disposal the numerical results for the helicity amplitudes of B d, s → γγ, we predict a number of observables of experimental interest including the CP-averaged branching fractions, the polarization fractions and the time-dependent CP asymmetries in the presence of the neutral-meson mixing. Our main observations and the concluding discussions will be presented in Section 7. We further present the detailed expressions of the complete NLO QCD corrections to the b → sγ matrix elements obtained in [17] and the renormalization-group evolution functions emerged in the NLL resummation improved factorization formulae of the double radiative B q -meson decay amplitudes at leading power in Appendices A and B, respectively.
2 Theory summary of the double radiative B q -meson decays

The effective weak Hamiltonian
Applying the classical equations of motion [36], the effective weak Hamiltonian for b → qγγ can be reduced to the one for the radiative b → qγ transitions at leading order in the Fermi coupling G F [37]. Employing the unitarity relations of the CKM matrix elements, the relevant effective weak Hamiltonian can be written as where the complete set of physical operators are given by [38] P p 1 = (q L γ µ T a p L ) (p L γ µ T a b L ) , P p 2 = (q L γ µ p L ) (p L γ µ b L ) , P 3 = (q L γ µ b L ) q (q γ µ q ) , P 4 = (q L γ µ T a b L ) q (q γ µ T a q ) , P 5 = (q L γ µ 1 γ µ 2 γ µ 3 b L ) q (q γ µ 1 γ µ 2 γ µ 3 q ) , P 6 = (q L γ µ 1 γ µ 2 γ µ 3 T a b L ) q (q γ µ 1 γ µ 2 γ µ 3 T a q ) , The sign convention of the dipole operators P 7,8 corresponds to the definition of the gauge µ with e f = −1 for the lepton fields [1]. In addition, m b (ν) stands for the bottom-quark mass at the renormalization scale ν in the MS scheme. The four-quark operator basis introduced in (2) is more advantageous than the conventional Buchalla-Buras-Lautenbacher scheme [39] due to the disappearance of the closed fermion loop with an odd number of γ 5 matrices in the interaction vertices. The renormalized effective Wilson coefficients C i (ν) at ν m b at NLL in QCD will be achieved by solving the renormalization-group equations including the three-loop mixing of the four-quark operators into P 7,8 and the two-loop self-mixing of the magnetic moment operators [38,40].

The helicity and transversity form factors
The exclusive radiativeB q → γγ decay amplitude can be expressed as where * 1 and * 2 correspond to the polarization vectors of the collinear and anti-collinear photons, respectively. We work in the rest frame of the B q -meson and define the two lightcone momenta n µ andn µ satisfying the relations n 2 =n 2 = 0 and n ·n = 2 by which allows us to write down the four-velocity vector v µ = (p + q) µ /m Bq = (n µ +n µ )/2. In analogy to the radiative leptonic B q → γ ¯ decays [41], the hadronic matrix element (3) at the lowest non-vanishing order in the electromagnetic interactions can be further brought into the following form The primary objective of the present paper consists in performing improved QCD calculations of the hadronic tensors appearing in (5) with the QCD factorization and the dispersion relation techniques. It is straightforward to derive the explicit expressions of T (p) i, αβ in the following where the electromagnetic current of the active quark is given by Applying the transversality conditions for the photon polarization vectors and the QED Ward-Takahashi identities, we can parameterize the obtained hadronic tensors in terms of the Lorenz-invariant helicity form factors where we have introduced the shorthand notations and conventions The transversity decay amplitudes of B q → γγ corresponding to the linearly polarized (anti)collinear photon states can be constructed along the line of the analogous analysis for the charmless non-leptonic B → V V decays [42] F (p) It is evident that the two-photon final states contributing to F i,⊥ are also the CP eigenstates with the eigenvalues +1 and −1, respectively. Furthermore, the left-handedness of the weak interaction Lagrangian and the helicity conservation for the strong interaction at high energy implies the well-known hierarchy structure 3 Factorization of the helicity form factors at leading power In this section we aim at constructing the perturbative factorization formulae for the helicity form factors at leading power in the heavy quark expansion in the SCET framework and accomplishing the complete NLL resummation of the parametrically large logarithms ln (m b /Λ QCD ) due to the renormalization-group evolutions of the hard matching coefficients and the leadingtwist B-meson distribution amplitude. The short-distance Wilson coefficients from integrating out the hard and hard-collinear fluctuations of the QCD correlation functions (6) and (7) will be determined by implementing the two-step matching program QCD → SCET I → SCET II in sequence. We start with constructing the SCET I representations of the hadronic tensors T (p) i, αβ in accordance with the position-space formalism [43,44] and then providing the detailed expressions of the hard matching coefficients at NLO in QCD. In contrast to the SCET factorization for the heavy-to-light B-meson form factors [45,46], only the A0-type effective weak currents O (A0) j will contribute to the double radiative decay amplitudes at leading power in the Λ QCD /m b expansion [47]. The obtained matching equation from QCD → SCET I can be written as with the chiral projection operator P L = (1 − γ 5 )/2. The desired SCET I representation of the electromagnetic current j em β, SCET I of our interest has been constructed in [47] j em β, where the soft and (anti)-collinear Wilson lines are defined in the standard way The position argument of the soft-quark field x + = (n · x/2) n arises from the multipole expansion. Apparently, the representation for an electromagnetic current carrying the collinear momentum can be obtained from (14) with the replacement rules hc ↔ hc and n ↔n.
The hard matching coefficients from the four-quark and choromomagnetic penguin operators at O(α s ) can be extracted from perturbative corrections to the QCD b → sγ matrix elements, which have been computed independently in [17,18] with different techniques for evaluating the non-trivial two-loop diagrams. Furthermore, we need the one-loop hard function from matching the electromagnetic dipole operator onto SCET I at leading power in the heavy quark expansion [45,46] with Here µ stands for the renormalization scale of the A0-type SCET I current, whereas ν characterizes the renormalization scale of the QCD tensor current. We then derive the NLO short-distance matching coefficients appearing in (13) as follows where the effective Wilson coefficients C eff i in the MS renormalization scheme and in the naive dimensional scheme for γ 5 are given by [38] The complete expressions of the two-loop functions F (p) i,7 (i = 1, ..., 6) and the one-loop function F 8,7 will be presented in Appendix A with no series expansion in the quark-mass ratio m c /m b .
We are now in a position to perform the perturbative matching of the SCET I correlation function in (13) onto SCET II by integrating out the hard-collinear fluctuations According to the definition of the HQET B-meson distribution amplitudes [48,49] we can readily derive the soft-collinear factorization formula for the effective non-local matrix element T αβ in momentum space The hard-collinear function from the matching procedure SCET I → SCET II has been already computed with distinct techniques at the one-loop accuracy [47,50] (see also [51] for the newly derived expression at two loops) Moreover, the HQET B-meson decay constantf Bq (µ) can be expressed in terms of the QCD decay constant f Bq by integrating out the strong interaction dynamics at the hard scalẽ The resulting factorization formula for the double radiative B q → γγ decay amplitude at leading power in the large energy expansion then reads Evidently, there is no single choice of the factorization scale µ to get rid of the enhanced logarithms of m b /Λ QCD in the factorized amplitude. Taking the factorization scale of order m b Λ QCD , implementing the resummation of ln(m b /Λ QCD ) for both the hard coefficient functions and the leading-twist B q -meson distribution amplitude to all orders in perturbation theory will be required in the NLL approximation. To achieve this goal, we will first employ the renormalization-group evolution equations for V (p) 7, eff and K −1 in momentum space, whose general solutions can be derived in the following form The explicit expression of the evolution functionÛ 1 can be obtained from U 1 (E γ , µ h , µ) displayed in [52] with the replacement E γ → m Bq /2. Furthermore, the evolution kernelÛ 2 can be deduced fromÛ 1 with the proper substitutions of the anomalous dimensions. The NLL resummation of large logarithms due to the scale evolution of the B-meson LCDA can be accomplished by virtue of the corresponding renormalization group equation at two loops [51,53] dφ The upper incomplete gamma function Γ(s, x) and the perturbative kernel h(x) due to conformal symmetry breaking are given by The perturbative expansions for the cusp anomalous dimension Γ cusp (α s ) and the non-logarithmic term γ η (α s ) are defined by [51,52] where the series coefficients of particular relevance for the NLL resummation are Applying the Mellin transformation with respect to the variable ω with a fixed reference scaleω, the general solution to the non-linear partial differential equation (28) in Mellin representation can be constructed explicitly [54] φ It proves helpful to introduce four dimensionless functions where the definition of the beta function reads [55] β Equivalently, one can perform an alternative integral transformation following [56] to reduce the momentum-space evolution equation (28) to an elegant form 1 Implementing the one-loop approximation for the general eigenfunctions of the Lange-Neubert evolution kernel [57] in the j-space displayed in (3.12) of [56] yields Q λ (j) = Γ(j + 2) exp [−j(λ − γ E + 1)], which precisely correspond to the momentum-space eigenfunctions Q s (ω) = 1 √ ω s J 1 (2 √ ω s) originally derived in [58] (see [59] for an independent construction with the aid of the collinear conformal symmetry) with the desired relation (µ e γ E s) −j = exp [− j (λ − γ E + 1)] at O(α s ) and the integral transformation (36).
The general expression of the renormalization kernel V (j, α s ) is given by where j(µ) is determined by the solution to the following differential equation The coefficient functions γ + (α s ) and ϑ(j) at the NLL accuracy take the form It is straightforward to demonstrate the equivalence of the two renormalization group equations (28) and (37) with the aid of the transformation relation The obtained solution (34) to the Mellin-space evolution equation of the B-meson distribution amplitude allows us to derive the NLL resummation improved expression for the leading-power factorization formula (26) exactlȳ where the operator-valued functionĴ is defined in terms of the hard-collinear function Figure 1: Diagrammatical representations of the QCD correlation functions (6) and (7) contributing to the double radiative B q → γγ decay amplitudes at LO in the strong coupling constant, where the accompanying diagrams due to the exchanges of two energetic photons in the final states are not presented.
Implementing perturbative expansions for the anomalous dimensions as well as the QCD βfunction at NLL leads to the approximate expressions of the renormalization group evolution functions displayed in Appendix B. Consequently, we can readily write down the factorized expressions for the helicity form factors at leading order in the power expansion of Λ QCD /m b where the explicit form of the dimensionful function R(m b , µ 0 , µ) can be found in (163).

Factorization of the helicity form factors at subleading power
In this section we proceed to construct the tree-level factorization formulae for the subleading power corrections to the exclusive radiative B q → γγ decay amplitude from distinct sources in the framework of QCD factorization. To achieve this goal, we will categorize the power suppressed mechanisms of interest according to the effective weak operators defining the QCD correlation functions (6) and (7) and further take advantage of the non-trivial relations for the light-cone HQET operators due to the classical equations of motion extensively.

The NLP corrections from the electromagnetic dipole operator
The first type of the next-to-leading power (NLP) correction originates from the subleading terms of the hard-collinear quark propagator due to the (anti)-collinear photon emission off the constituent light-flavour quark of the B q -meson as displayed in Figure 1(a). Such NLP contribution from the electromagnetic penguin operator P 7 can be determined from the QCD correlation function (6) T hc, NLP where the symmetry factor "2" arises from the exchanges of two identical particles in the final states. Applying the light-cone decomposition of the Dirac matrix γ µ [60] and employing the general parametrization of the hadronic matrix element (45), one can immediately demonstrate that the yielding contributions from the second and third terms in (46) vanish in the four-dimensional space. We are therefore led to the reduced formula with the aid of the standard method of the integration by parts (IBP) and the precise relation of the four-vectorsn µ = 2 v µ − n µ . Taking advantage of the well-known operator identity due to the classical equations of motion [61,62] v where the derivative with respect to the total translation is defined by the subleading power HQET matrix element (47) can be further computed as follows where the hadronic parameterΛ characterizes the mass splitting between the heavy quark and the physical heavy-hadron state under discussion [63] Moreover, the coordinate-space factorization formula (50) at LO in the strong coupling has been established by implementing the systematic parametrization for the light-cone hadronic matrix element of the three-body effective operator at the twist-six accuracy [30] 0|q The emerging eight independent invariant functions can be expressed in terms of the HQET distribution amplitudes with the definite collinear twist In the light of the momentum-space representations of the two-particle and three-particle the resulting factorized expression for the NLP amplitude (47) can be written as where the first inverse moment of the leading-twist B q -meson LCDA fulfills the established relation under the renormalization-scale evolution in the LO approximation (see [54] for the NLL improvement) We have introduced the following definition of the inverse-logarithmic moments σ (n) Accordingly, the obtained subleading power corrections to the two helicity form factors reads Apparently, the NLP hard-collinear contribution defined by (45) preserves the large-recoil symmetry F i, ⊥ of the transversality amplitudes, in accordance with the analogous observation for the power suppressed effect in B q → γ ¯ due to the B-type insertion of P 7 [41].
In addition, it would be interesting to investigate the subleading power correction to the hadronic matrix element T 7,αβ from the non-vanishing light-quark mass. A straightforward calculation at tree level leads to the following expression where the convolution integral over the variable ω is unfortunately divergent based upon the asymptotic behaviour of φ − B (ω, µ). Employing the prescription described in [64], we parameterize this contribution by introducing an unknown complex quantity X NLP where an arbitrary perturbative scale Λ UV is introduced to guarantee the ultraviolet finiteness of the newly defined phenomenological parameter. For concreteness, the nonperturbative model for the coefficient X NLP adopted in this work reads (see [65] for further discussions) The achieved results of the helicity form factors are given by which implies the anticipated scaling rule in the heavy quark expansion We now turn to evaluate the NLP correction from the power suppressed term in the SCET I representation of the heavy-to-light transition current appearing in (6) which arises from the HQET expansion of the heavy quark field in QCD [66] b As we will not go beyond the accuracy of Λ QCD /m b in the derivative expansion, the lowestorder equation of motion in HQET (see [67] for the generalization at NLO) (65). The corresponding NLP correction to the correlation function (6) can be derived as follows Implementing further the classical HQET operator identities [61,62] q the subleading power amplitude T A2, NLP 7, αβ can be cast in the following form with Applying the definitions of the two-particle and three-particle B q -meson distribution amplitudes (22) and (52) and carrying out the essential integrations immediately gives rise to We explicitly verify that substituting the covariant derivative − → / D by the transverse derivative − → / D in the tree-level expression (68) indeed yields the identical factorization formula. It is then straightforward to identify the resulting NLP corrections to the helicity form factors We are now in a position to compute the subleading power contributions to the QCD matrix element (6) from the higher-twist B q -meson distribution amplitudes with the perturbative factorization technique. As emphasized in [68,69], the consistent QCD calculations of the higher twist corrections will necessitate taking into account the non-minimal partonic configurations with additional quark-gluon fields and the non-vanishing partonic transverse momenta in the leading Fock state simultaneously. Including the off-light-cone corrections for the non-local two-body HQET operators up to the O(x 2 ) accuracy [30] we can readily derive the factorized expression of the yielding higher-twist contribution Employing the non-trivial relation of the coordinate-space distribution amplitudeŝ we can express the twist-four LCDA g + B (ω, µ) in terms of the "genuine" three-particle distribution amplitude of the same twist and the Wandzura-Wilczek contribution [70] related to the lower twist two-particle distribution amplitudes The manifest expression ofḡ Substituting the established relation (77) into the tree-level factorization formula (75) leads to the equivalent form in the following On the other hand, the subleading twist three-particle corrections to the vacuum-to-B q -meson matrix element (6) can be computed from the light-cone expansion of the quark propagator in the background soft-gluon filed [72] 0|T Applying the general parametrization of the light-cone HQET matrix element (52) allows us to write down the factorized expression of the higher-Fock-state contribution Putting together the achieved two-particle and three-particle higher twist contributions, we can readily derive the factorization formula for the resulting NLP correction . (82) According to the momentum-space relation for the HQET distribution amplitudes [58] we are then able to rewrite the first convolution integral appearing in (82) as follows The newly defined distribution amplitude ∆φ − B (ω, µ) can be further constructed from the "genuine" twist-three nonperturbative function φ − tw3 where the Wandzura-Wilczek term φ − WW

B
(ω, µ) can be expressed in terms of the leading twist The integral representation of the higher twist distribution amplitude (Ψ 4 −Ψ 4 ) motivated by the corresponding renormalization group equation suggests the following decomposition where the twist-three contribution can be deduced from We can then derive an equivalent form of the second convolution integral in (82) which together with the obtained relation (84) leads to Implementing the exact identity for the coordinate-space distribution amplitudes we are ready to establish the momentum-space relation accordingly which implies another useful representation of the higher twist factorization formula Neglecting the nonleading partonic configuration corrections from the four-body light-cone HQET operators of the typesq G G h v andq qq h v facilitates the construction of the approximate expression for the twist-four three-particle distribution amplitudes We are then led to an even more compact factorization formula at tree level in analogy to the factorized expression for the subleading twist corrections to the radiative leptonic B → γ ν decays as displayed in [35]. As a consequence, the resulting higher twist corrections to the helicity form factors can be written as Furthermore, the subleading power contribution to the QCD correlation function (6) due to the electromagnetic current of the bottom quark can be computed as follows which evidently generates the local NLP corrections to the two helicity form factors It is perhaps worth mentioning that the (anti-)collinear photon emission from the heavy quark gives rise to the power suppressed correction preserving the large-recoil symmetry relation of the transversality amplitudes, in contrast to the counterpart contribution to the radiative decay B → γ ν form factors [52].

The NLP corrections from the four-quark operators
Bearing in mind that we only aim at evaluating the subleading power contributions to the double radiative B q -meson decay form factors at O(α 0 s ) accuracy, it is therefore unnecessary to take into account the NLP correction to the QCD correlation function (7) defined by the chromomagnetic dipole operator P 8 . By contrast, there exist two different types of the power suppressed contributions to the hadronic matrix elements of the four-quark operators already at tree level. On the one hand, the factorizable quark-loop diagrams with an insertion of the QCD penguin operators in Figure 1(b) yield the NLP contributions analogous to the discussions presented in Section 4.1 2 , which can be readily achieved by replacing the Wilson coefficient C 7 (µ) multiplying each individual subleading power correction with the effective coefficient C eff 7 (µ). On the other hand, the QCD matrix elements of the four-quark operators further generate the weak-annihilation-type of the subleading power corrections as displayed in Figure 1(c), with no energetic photons coupling to the partonic constituents of the initial B q -meson state directly. The resulting factorization formulae of the weak-annihilation contributions to the two helicity form factors ofB q → γγ can be written as The explicit expressions of the one-loop hard matching coefficient functions with contributions from a complete set of the hadronic operators are given by where we have introduced the following conventions for perturbative loop functions with It is apparent that the subleading power weak annihilation contributions from the QCD penguin operators will spoil the large-recoil symmetry relation between the two transversity amplitudes perturbatively. In particular, the non-trivial strong phase of the perturbative kernel H LL (m c ) arises from the discontinuities in the variable (p + q) 2 , which can be understood from the final-state rescattering mechanismB q → H c H c → γ γ at hadronic level with H c and H c standing for the appropriate (anti-)charm-hadron states [26,27] (see also [73] for further discussions in the context of B → K ). We also mention in passing that the weak annihilation diagrams involving the light-quark loop will not generate the perturbative strong phase at O(α 0 s ) in the limit m q → 0, due to the helicity and angular momentum conservations.

Final factorized expressions for the NLP corrections
We now summarize the tree-level factorization formulae of the various subleading power corrections to the helicity form factors considered so far which together with the NLL resummation improved factorization formulae (44) at leading power in the heavy quark expansion constitute the main technical results of this paper.

The resolved photon corrections from the dispersion approach
The major objective of this section is to compute the power suppressed soft contributions to the double radiative B q -meson decay amplitudes with the OPE-controlled dispersion technique. We will implement the nonperturbative prescription originally constructed for accessing the long-distance photon correction to the pion-photon form factor [31] by investigating the QCD correlation function responsible for the more general transition B q → γ * γ with an off-shell and transversely polarized photoñ It is then straightforward to derive the leading power factorization formula for the generalized form factorF 7,L (n · q, n · q) in the Λ QCD /m b expansion for the hard-collinear q 2 regimẽ where the short-distance matching coefficientJ at the one-loop accuracy is given by [34] Isolating the contributions of the lowest-lying hadronic states and applying the standard definitions of the vector-meson decay constant and the heavy-to-light B q -meson decay form factors we can further derive the hadronic dispersion relation for the hadronic tensor (105) with the aid of an exact QCD identity T where the non-vanishing coefficients a Matching the two different representations of the correlation function (105) with the partonhadron duality ansatz for the physical spectral density and performing the Borel transformation with respect to the variable n · q leads to the desired sum rules for the tensor B q → V decay form factor in QCD The factorized expression (106) allows us to construct the QCD spectral function immediately with the master formulae collected in Appendix A of [34] (see also [35]) with the effective heavy-meson "distribution amplitude" encoding both the soft and hardcollinear strong interaction dynamics The resulting predictions from the SCET sum rules (112)  with the intervals of the Borel mass and the threshold parameter determined in [74] are explicitly verified to be compliable with the numerical results obtained from an alternative LCSR approach [75]. Substituting the obtained sum rules (112) into the hadronic dispersion relation (109) and taking the light-cone limit for the four-momentum q µ →n·q n µ /2 gives rise to the nonperturbatively refined spectral representation where the first and third lines precisely correspond to the SCET factorization formula of the leading power contribution for the correlation function (6). Consequently, we can readily identify the resulting NLP corrections to the helicity form factors of B q → γγ as follows According to the power counting scheme for the effective threshold and the Borel parameter the peculiar NLP mechanism responsible for (116) indeed arises from the long-distance correlation of the electromagnetic current j em β (x) and the dipole operator P 7 (0) with x 2 ∼ 1/Λ 2 QCD , which cannot be probed with the standard light-cone OPE technique without introducing the additional nonperturbative distribution amplitudes. Inspecting the factorization property of the subleading power correction in question with the heavy-meson and (anti)-collinear-photon LCDAs simultaneously indicates that the soft-collinear convolution integrals appearing in the factorized expression for the effective matrix element of the yielding A0-type SCET I current develop rapidity divergences, in analogy to the counterpart for the semileptonic B q → V form factors with a transversely polarized vector meson [49,76]. Under this circumstance, evaluating the NLP "resolved" photon corrections in QCD can be addressed independently by employing the LCSR method with the photon distribution amplitudes [77], as already proposed in [12] in the context of B → γ ν . We are now ready to summarize the power suppressed soft contributions from both the magnetic penguin operators and the four-quark operators where the effective hard function V (p) 7, eff at the accuracy of O(α s ) has been presented in (19). Collecting the different pieces together, the yieldingB q → γγ amplitude can be eventually expressed in terms of two independent helicity amplitudes where the manifest expressions ofĀ L andĀ R can be derived in the following 4 In order to facilitate the phenomenological explorations of the time-dependent decay observables in the presence of the neutral meson-antimeson mixing, we further introduce the transversity amplitudes with definite CP transformation properties which allows for writing the exclusive radiative decay amplitude (119) in an alternative form

Numerical analysis
Having at our disposal the improved results of the helicity form factors including both the NLL resummation of the leading power contributions and the newly obtained NLP corrections, we are now ready to investigate their numerical implications on a number of observables for the exclusive radiative B q → γγ decay processes of experimental interest. To this end, we will proceed by specifying the different types of theory inputs (Wilson coefficients, quark masses, B q -meson LCDAs, CKM parameters), which will be essential to determine the two independent helicity amplitudesĀ L andĀ R numerically.

Theory inputs
For definiteness we will employ the three-loop running of the strong coupling constant in the MS scheme by adopting the interval of α The Wilson coefficients C i (ν) in the weak effective Hamiltonian (1) for the hard scale ν m b at the NLL accuracy will be evaluated from the NLO matching expressions at the initial scale 4 For later convenience we further introduce the shorthand notations for the separate amplitudes   which are in excellent agreement with the numerical values displayed in [2] (see also [41]) with slightly different choices of the SM input parameters. The MS mass of the bottom quark in the magnetic dipole operators P 7, 8 collected in Table 1 is determined from the lattice QCD simulations with N f = 2+1+1 flavours of sea quarks. By contrast, the bottom-quark mass appearing in the perturbative hard and hard-collinear functions from QCD → SCET I → SCET II matching is usually understood to be the pole mass due to on-shell kinematics [1,2,45]. However, converting the precisely known MS-scheme heavy quark mass to the counterpart pole-scheme mass will lead to the numerical results sensitive to the loop order of the matching relation in the corresponding perturbative expansion. We will therefore use the potential-subtracted (PS) renormalization scheme [85] (see [86] for an alternative short-distance quark-mass scheme) for both the bottom-and charm-quark masses consistently in the short-distance functions displayed in (17), (19), (24), (25), (101) and (114) following the prescription of [1,41]. We now turn to discuss the nonperturbative hadronic inputs entering the factorized expressions of the B q → γγ helicity amplitudes (120). The leptonic decay constants of the pseudoscalar B d -and B s -mesons in QCD are taken from the lattice averages of N f = 2 + 1 + 1 results in the isospin-symmetry limit [80] (see [87] for further discussions on the strong-isospin breaking effect due to m u = m d and [88,89] for the systematic strategies to take into account the electromagnetic corrections). In addition, the two-particle and three-particle B q -meson distribution amplitudes in HQET apparently serve as the fundamental ingredients of the derived factorization formulae, encoding the strong interaction dynamics from the soft-scale fluctuation, for both the LO and NLO terms in the heavy quark expansion. Following [35] we will introduce the general three-parameter ansatz for the leading-twist LCDA φ + B (ω, µ 0 ) which can be analytically evolved into the hard-collinear scale µ hc under the renormalizationgroup flows at the LL accuracy in terms of the generalized hypergeometric 2 F 2 (a, b; c, d; x) functions. The shape parameter λ Bq defined in (56) and the associated inverse-logarithmic moments σ (n) Bq (58) can be determined straightforwardly with the following identities In spite of the intensive investigations of understanding the dynamical aspects of the twist-two B q -meson distribution amplitude (as well as the counterpart transverse-momentum-dependent wavefunction) with distinct techniques and strategies [12,34,35,52,[92][93][94][95][96], the current theory constraints on λ Bq (µ 0 ) still turn out to be far from satisfactory due to the uncontrollable systematic uncertainties. Consequently, we will vary the input value of λ B d (µ 0 ) in the conservative interval as presented in Table 1 and further assign approximately O(15 %) SU(3)-flavour symmetry breaking effect 5 for the ratio λ Bs (µ 0 )/λ B d (µ 0 ) as recently advocated in [41]. Moreover, we will employ the phenomenological models for the subleading-twist heavymeson distribution amplitudes satisfying the asymptotic behaviours at small quark and gluon momenta and the classical equations of motion constraints [35] (see [30,71] for the particular two-parameter models of the following ansatz) where the normalization condition and the first three moments of the function f (ω) read The HQET parameters λ 2 E and λ 2 H are defined by the hadronic matrix element of the threebody effective local operator [30,48] Applying their renormalization-group equations at the one-loop order [90,91] we can readily write down the corresponding solutions in the LL approximation where the matrixV diagonalize the evolution kernel γ EH with the yielding eigenvalues as irrational numbers remarkably The available QCD sum rule predictions for these two dimensionful quantities {(0.03 ± 0.02) GeV 2 , (0.06 ± 0.03) GeV 2 } , [91] (133) are observed to differ significantly in magnitude due to the sizeable radiative correction to the quark-gluon condensate contribution and the yet higher power correction from the dimensionsix operator qq 2 included in the improved analysis [91]. Taking advantage of the default ansatz (127) for the subleading twist B q -meson distribution amplitudes allows for the enormous simplification of the factorized expressions for the power suppressed contributions at tree level derived in Section 4.1 with the interesting relations which are apparently independent of the precise shape of the profile function f (ω). Keeping in mind that the predicted ratio λ 2 E (µ 0 )/λ 2 H (µ 0 ) shown in (133) is insensitive to the higher-order corrections in the sum rule calculations and the normalization constant κ(µ 0 ) is actually determined by the particular combination 2 λ 2 E (µ 0 ) + λ 2 H (µ 0 ), it proves to be most advantageous to take the aforementioned two quantities as the independent hadronic inputs parametrizing the chromoelectric and chromomagnetic matrix elements rather than λ 2 E (µ 0 ) and λ 2 H (µ 0 ) themselves as already discussed in [35]. The renormalization-scale dependence of the first logarithmic moment σ (1) Bq (µ) entering the expression of I mq NLP can be readily determined from the evolution equation of the HQET distribution amplitude φ + B (ω, µ) [57] σ (1) which can be also derived from the one-loop evolution equation of the moment σ Bq (µ) presented in [58,93] alternatively. Following [89], we will set the "effective mass" of the B q -meson state appearing in the subleading power factorization formulae (59) and (73) asΛ = m Bq − m b with m b = (4.8 ± 0.1) GeV numerically (see [98,99] for further discussions on the scheme dependence of this HQET quantity).
In addition, we will vary the hard scale µ h of the short-distance matching coefficient V

Theory predictions for the helicity amplitudes
We now proceed to explore the phenomenological impacts of the newly derived NLP corrections at O(α 0 s ) and the NLL resummation improved leading power contributions on the exclusive radiative B q → γγ helicity amplitudes. To develop a transparent understanding of the higher order corrections from distinct dynamic mechanisms, we display in Figure 2 the λ Bq dependence of the leading power contributions at the LL and NLL accuracy, the combined local and non-local NLP corrections (104) from the QCD factorization approach and the subleading power soft contributions estimated with the dispersion technique for both b → d γ and b → s γ transitions, including the theory uncertainties from the variations of the hard scales µ h and ν as well as the hard-collinear scale µ. Evidently, the radiative corrections with renormalization-group improvement can significantly reduce the perturbative scale uncertainties for the LL predictions of the helicity amplitudeĀ L (B q → γγ) at leading power in the Λ QCD /m b expansion. In general, the perturbative QCD correction to ReĀ L (B q → γγ) (with q = d, s) at O(α s ) can shift the corresponding LL prediction by an amount of approximately (2 ∼ 10)% for the allowed interval of λ Bq 6 . Moreover, the yielding LL and NLL scale uncertainty bands only overlap marginally for ImĀ L (B d → γγ) and even turn out to be well separated for ImĀ L (B s → γγ), while adopting the central values for all the remaining input parameters. The latter can be attributed to the fact that the negligible LL contribution to ImĀ L (B s → γγ) only arises from the heavily suppressed CKM factors Im [V ub V * us ] and Im [V cb V * cs ] in comparison with the counterpart amplitude ReĀ L (B s → γγ), while the dominant NLL correction to ImĀ L (B s → γγ) originating from the emerged strong phase of the two-loop contribution of the current-current operator P c 2 is apparently free of such CKM suppression. In addition, the power suppressed soft corrections to both ReĀ L (B q → γγ) and ImĀ L (B q → γγ) will gives rise to the destructive interferences with the appropriate leading power contributions for the central input parameters, in analogy to the previous observation forB → γ ν [34,35]. By contrast, the factorizable NLP correction to the amplitude ReĀ L (B q → γγ) can flip the sign numerically with the evolution of the inverse moment λ Bq due to the nontrivial competing mechanisms between the two distinct sectors of {Ā A2, NLP In addition, it has already become evident from the factorization formula (104) that the right-handed helicity amplitudeĀ R (B q → γγ) results from the local subleading power con- completely, in contrast to the large-recoil symmetry breaking effects forB q → γ ¯ [41]. The sizeable uncertainties of the yielding predictions for ReĀ R (B q → γγ) and ImĀ R (B q → γγ) reflect the uncancelled renormalization-scale dependence of the Wilson coefficients C i (ν) for the tree-level calculation described in Section 4.2, which is expected to be compensated by the unknown ν-dependent hard functions at O(α s ) requiring the explicit evaluations of the two-loop one-particle irreducible diagrams depicted in Figure 4 of [14]. It needs to be pointed out further that the obtained results of ReĀ R (B d → γγ) and ImĀ R (B d → γγ) appear to own different signs from the counterpart predictions forB s → γγ due to the contributing CKM matrix elements Re [V cb V * cd ] < 0 and Re [V cb V * cs ] > 0. Inspecting the various NLP corrections to the factorizableB q → γγ amplitudes displayed in Figure 3 numerically reveals the notable sensitivity of the real parts forĀ hc, NLP L ,Ā mq, NLP L andĀ A2, NLP L on the actual value of λ Bq , which together with the remaining NLP corrections independent of the inverse moment characterizes the diverse facets of the strong interaction dynamics encoded in the subleading power heavy-hadron decay amplitudes. It is also worth mentioning that the non-local strange-quark mass correction to ReĀ L (B s → γγ) is of mithe perturbative effects entering the two short-distance functions C eff 7 C   for the radiativē B d → γγ decay process together will fall short of the magnitude of the counterpart weakannihilation contribution ImĀ WA, NLP L consistently. We further present in Figure 4 the final theory predictions for the helicity amplitudes of B q → γγ, as the analytical functions of the inverse moment λ Bq , with the combined uncertainties by adding in quadrature all the separate uncertainties from varying the input parameters discussed in Section 6.1. It is straightforward to identify the dominant theory uncertainties for the left-handed helicity amplitudes arising from the hard and hard-collinear scale fluctuations and the inverse-logarithmic moments σ (1) Bq and σ (2) Bq , in addition to the particular sensitivity to λ Bq 7 . In contrast, the yielding theory uncertainties for the right-handed helicity amplitudes are dominated by the hard-scale variations of the Wilson coefficients C i (ν) entering the tree-level factorization formulae of the weak-annihilation corrections. Moreover, our limited knowledge of the nonperturbative coefficient X NLP parameterizing the rapidity divergence of the NLP correction (63) will not generate the sizeable uncertainty due to the smallness of the ratio m q /λ Bq on top of the λ Bq /m Bq suppression in comparison with the counterpart charmless two-body B q -meson decays in QCD factorization [64]. will only bring about the subdominant uncertainties, since the numerically important contributions stem from the constant factors and the "twist-two" terms proportional toΛ/λ Bq in the curly brackets of (59), (73) and (97).

Phenomenological observables for B q → γγ
Prior to addressing the phenomenological implications of the improved calculations of the double radiative decay amplitudes on the experimental observables, we will construct the systematic expressions for investigating the CP-averaged branching fractions, the polarization fractions and the time-dependent CP asymmetries taking into account the neutral-meson mixing based upon the standard formalism presented in [78,102]. In doing so, we are required to derive the decay amplitude for the counterpart channel B q → γγ by applying the CPtransformation for the amplitude (123) where we have adopted the phase convention CP |B q = −|B q and the two transversity amplitudes A and A ⊥ can be determined from the corresponding "barred" amplitudes with all the weak phases conjugated. Solving the Schrödinger equation governing the time evolution of the B q −B q system immediately leads to the time-dependent decay rates [103,104] where N f stands for the time-independent normalization factor arising from the phase-space integration and the index χ is introduced to characterize the parallel and perpendicular polarization configurations of the two-photon states. Apparently, the two decay amplitudes A (B q → γγ) and A ⊥ (B q → γγ) for the final states possessing the definite CP-parities can be obtained from (140) by keeping the individual terms proportional to A and A ⊥ , respectively. The average width Γ q of the two bottom-meson mass eigenstates depends upon the diagonal matrix elements of the decay matrix with Γ q 11 = Γ q 22 thanks to the CPT invariance of the effective Hamiltonian dictating the neutral-meson mixing pattern. The mass and width splittings between the light and heavy mass eigenstates can be expressed in terms of the dispersive and absorptive matrix elements M q 12 and Γ q 12 . The B q -mixing observables ∆m q and ∆Γ q further facilitate the construction of the two important quantities whose numerical results with uncertainties from the world average values of the experimental measurements [78] (with the combination procedures described in [105]) 8 have been already summarized in Table 1. The small quantity a q characterizing the CP violation in mixing can be defined by the nondiagonal matrix elements in the following form The SM predictions for such flavour specific asymmetries [107,108] imply that we can safely drop out the tiny corrections of O(a n q ) (with n ≥ 1) for the radiative B q → γγ decay phenomenology. The phase-convention independent combination λ f and the appearing three CP violation observables A dir, χ CP , A mix, χ CP and A χ ∆Γ are further defined as which evidently reveals an exact relation Expanding the obtained expression (143) for the phase-convention dependent quantity q/p in terms of the small parameter |Γ q 12 /M q 12 | leads to an approximate solution [106] The SM contributions of the off-diagonal matrix element M q 12 turn out to be dominated by the ∆B = 2 box diagrams with internal top and anti-top quarks due to the very strong GIM cancellation [109] for the yielding two terms proportional to λ 2 c and λ c λ t (with λ p = V pb V * pq ), which can be traced back to the peculiar analytical behaviours of the associated perturbative loop function F m 2 107,110]. We are then led to conclude that the imaginary part of M q 12 originates to a rather good approximation from the CKM factor V tb V * tq uniquely with the standard OPE technique, which further enables us to write down an even illuminating relation for the B q -meson mixing relation y q : x q ≈ y s : x s [78,106], which holds to the first order in |Γ q 12 /M q 11 | by stimulatingly switching off the SU(3)-flavour symmetry violating effects of the "bag" parameters and the subleading power corrections to the relevant ∆B = 2 matrix elements in the heavy quark expansion (see [107] for an overview of the current theory status).
According to (141), we can immediately derive the CP-averaged branching fraction of B q → γγ for the flavour-tagged measurement at an e + e − collider with equal numbers of the produced bottom and anti-bottom mesons [111,112] Interestingly, the resulting correction factor of the time-integrated branching ratio due to the neutral-meson mixing is independent of the polarization index χ for the coherent bottommeson pair production at the SuperKEKB accelerator. The CP-averaged branching fraction in the limit ∆Γ q → 0 can be readily computed from the corresponding transversity amplitudes Consequently, we can proceed to construct a typical set of phenomenological observables in analogy to the two-body hadronic B → V V decays [42] BR where the two polarization fractions are apparently not linearly independent due to the normalization condition f + f ⊥ = 1. Applying the definition for the time-dependent CP asymmetry for the neutral B q -meson decaying into CP eigenstates yields [102] A χ CP (t) =Γ The quantity A dir, χ CP describes the direct CP-violating effect as a consequence of Ā χ : A χ = 1, which requires the appearance of at least two individual terms with distinct strong and weak phases simultaneously contributing to the transversity amplitude A χ . By contrast, A mix, χ CP encodes the CP-violating information due to interference between decays with and without mixing for the final states common to both the bottom and anti-bottom meson decays (i.e., by virtue of the two different decay mechanisms of B q → γγ and B q →B q → γγ). The observable A χ ∆Γ 9 takes the extreme values ±1 for the so-called "golden modes" with CP-eigenstate final states, whose decay amplitudes contain only one CKM structure such that |λ χ f | = 1. We now present the improved theory predictions for the CP-averaged branching fractions of B q → γγ in the presence of the neutral-meson mixing in Figure 5, which summarizes our main phenomenological results for confronting the forthcoming precision Belle II measurements with an ultimate integrated luminosity of 50 ab −1 [15]. The factorizable NLP correction to BR(B q → γγ) appears to reduce the counterpart leading power contribution by approximately an amount of (15 − 20) % at λ B d = 200 MeV (λ Bs = 250 MeV). However, such destructive interference mechanism will gradually evolve into the constructive intervention pattern with the increase of λ Bq (as large as 35 % enhancement at the maximal value of the inverse moment) due to the emergence of the zero-crossing regime for the factorizable NLP correction to the amplitude ReĀ L (B q → γγ) as discussed in Section 6.2 together with the hierarchy structure where the quoted uncertainties arise from the variations of λ Bq merely. The NLP resolved photon contribution to the branching fraction of B q → γγ will result in the sizeable reduction of the corresponding "LP + fac, NLP" prediction with the default ansatz of the B q -meson distribution amplitude (125) at α = β (namely, the exponential model [30,48,71]): numerically (6 − 28) % and (15 − 43) % corrections for b → d γ and b → s γ transitions, respectively. 9 As pointed out in [113], the mass-eigenstate rate asymmetry A χ ∆Γ is also essential to determine the effective lifetime for the untagged B q → f decay.  Inspecting the distinct sources of the yielding theory uncertainties as collected in Tables 2 and  3, the currently poor understanding towards the two-particle and three-particle B q -meson distribution amplitudes in HQET remains the most significant ingredient of preventing us from achieving the phenomenologically encouraging precision up to date, which can be competitive with the experimental prospects of the anticipated Belle II measurements (with the estimated precision of 9.6 % and 23 % for BR(B d → γγ) and BR(B s → γγ) approximately).

Central Value Total Error
Remarkably, including the subleading power soft corrections to the time-integrated decay rates will enlarge the combined theory uncertainties by almost a factor of two for both exclusive radiative decay processes with the central values of λ Bq , which can be attributed to the even more pronounced sensitivity to the shape parameters σ (1) B d and σ (2) B d under this circumstance. We can readily understand such interesting observation from the fact that the two inverselogarithmic moments will start to emerge in the leading power factorization formulae for the helicity amplitudes by means of either the perturbative correction to the hard-collinear matching coefficient or the renormalization-group evolution of the inverse moment λ Bq , both of which will lead to an effective suppression factor α s (µ) C F /(4 π) ≈ 0.04 at µ = 1.5 GeV numerically. On the contrary, the NLP resolved photon contribution (116) estimated with the dispersion approach depends on the functional form (particularly the small-ω behaviour) of the leading twist B q -meson distribution amplitude φ + B (ω, µ) already at tree level (see [93,114] Tables 2 and 3 are in order.
• It is instructive to derive the transparent expression of the ratio for the two branching fractions of B q → γγ in the leading-order approximation of the double expansion in powers of Λ/m b and α s We have verified explicitly that the λ Bq -scaling violation effect at leading power in the heavy quark expansion appears to be insignificant numerically (at the level of O(1 %)); while the subleading power factorizable and soft contributions can result in approximately (10 − 20) % corrections to the scaling behaviour (154) for the default intervals of the two inverse moments. It is therefore of high interest to carry out the precision measurements for such "golden observable" BR(B s → γγ) : BR(B d → γγ) at the Belle II experiment, which allow for an excellent determination of the ratio of the two in-verse moments with the systematic theory uncertainties at the level of (5 − 10) % in the leading-power approximation.
• In comparison with the previous predictions of the exclusive radiative B q → γγ decay rates from the QCD factorization approach [16] our numerical results for the corresponding leading power contributions with the central values of the input parameters turn out to be smaller by a factor of 2 − 3 approximately. The observed discrepancies in large part stem from approximating the MS-mass of the b-quark by the bottom-meson mass for the resulting hard-scattering kernel (19) in [16], from adopting the higher value of the effective Wilson coefficient C eff 7 (4.8 GeV) = −0.390 and from taking the lower value of the inverse moment for the B s -meson distribution amplitude λ Bs (µ 0 ) = 350 MeV implemented in the phenomenological analysis of [16].
• Adding the factorizable NLP weak-annihilation correction on the top of the leadingpower contribution from the magnetic penguin operator P 7 and adopting further the central input parameters gives rise to [14] BR which differ prominently from our numerical predictions with the same theory accuracy The dominating mechanisms accounting for the notable differences between (156) and (157) lie in the replacement of the MS-scheme bottom-quark mass in the effective weak operator P 7 by the counterpart heavy-meson mass m Bq for deriving the leading-power transversity amplitudes (3.2) in [14], which invokes additionally the approximate SU(3)flavour symmetry relation λ Bs = λ B d in the numerical extrapolation.
We proceed to present our predictions for the polarization fractions of the double radiative B q -meson decays within the default intervals of the inverse moments in Figure 6. Apparently, the yielding uncertainties for such ratio observables are substantially improved, when compared with the numerical results of the CP-averaged branching fractions displayed in (5), thanks to the striking cancellation of the errors from varying the shape parameters of the bottom-meson distribution amplitudes. The most significant uncertainties for the polarization fractions arise from the variations of QCD renormalization scale ν as shown in Tables  2 and 3. It is appropriate to emphasize that the deviation of the quantity f : f ⊥ from one evidently characterizes the typical size of the subleading power correction in the Λ QCD /m b expansion. In addition, the SU(3)-flavour symmetry violation effects for the polarization fractions are estimated to be extraordinary small, around (4 − 7) % numerically, on account of an extra suppression factor from the heavy quark mass expansion.
Finally we turn to display in Figure 7 the theory predictions for the three observables determining the time-dependent CP asymmetries of B q → γγ with the nonvanishing width splittings ∆Γ q . The apparent hierarchy structure for the CKM matrix elements contributing to the B s → γ γ helicity amplitudes indicates the model-independent results of the mass-eigenstate rate asymmetries A ∆Γ −1 and A ⊥ ∆Γ +1 by dropping out the negligible V ub V * us terms in (119), which are further supported by the complete NLL QCD calculations at leading power in Λ QCD /m b with the various NLP corrections at tree level as summarized in Table 3 Tables 2 and 3 appear to be more uncertain, mainly from varying the inverse moment λ Bq as well as the perturbative hard scale ν. According to the definition of A mix, ⊥ CP in (145), we can readily identify the emerged considerable uncertainties from the strong sensitivity of the Wilson-coefficient combination [C F C 1 (ν) + C 2 (ν)] entering the factorizable weak-annihilation amplitude F (p), WA V in (101) on the renormalization scale ν. As can be understood from the factorization formula (42), the two direct CP-violating quantities A dir, CP and A dir, ⊥ CP are expected to be identical in the leading-power approximation. However, the obtained numerical results displayed in Tables 2 and 3 tend to suggest that the subleading power contributions to the radiative B q -meson decay amplitudes in the heavy quark expansion can numerically generate the substantial corrections to the above-mentioned asymptotic scaling. This pattern can be attributed to the fact that the higher order terms in the heavy quark mass expansion will generally lead to the enormous impacts on the perturbative strong phases for the exclusive B q -decay amplitudes due to the absence of the α s suppression.

Conclusions
In the present paper we have carried out the improved calculations of the double radiative B q -meson decay amplitudes beyond the leading-power accuracy by taking advantage of the perturbative factorization technique and the dispersion approach. The NLO hard matching coefficients from the two-loop QCD matrix elements of q γ|H eff |b including the QCD penguin operators P 3,...,6 were taken into account for the first time based upon the analytical calculations performed in [17], which are apparently of importance for exploring the phenomenological observables of B d → γ γ. A complete NLL resummation for the enhanced logarithms of m b /Λ QCD appearing in the SCET factorization formula for the radiative B q → γ γ decay amplitude has been accomplished with the renormalization group formalism in momentum space. Subsequently, we constructed the factorized expressions for the NLP contributions from the non-local hard-collinear interaction due to the energetic photon radiation off the light-flavour quark of the bottom meson, from the nonvanishing light-quark mass, from the subleading power heavy-to-light current J (A2) in SCET, from the two-particle and three-particle B q -meson distribution amplitudes up to the twist-six accuracy, from the (anti)-collinear photon emission from the bottom quark and from the one-particle irreducible weak-annihilation diagrams at tree level in detail in Section 4. It is remarkable to mention that all the above-mentioned nonlocal NLP corrections other than the light-quark mass effect (60) can be expressed in terms of the convolution integrals of the hard-collinear functions and the HQET distribution amplitudes free of the potential rapidity divergences, which would invalidate the resulting soft-collinear factorization formulae otherwise. We further evaluated the NLP resolved photon contribution to the two helicity form factors with the OPE-controlled dispersion technique at O(α s ). The large-recoil symmetry violating effect between the two transversity amplitudes has proven to be solely generated by the local NLP contribution by means of the weak-annihilation mechanism.
Adopting the three-parameter ansatz for the essential B q -meson distribution amplitudes [35] fulfilling the classical equation-of-motion constraints and the appropriate asymptotic behaviors at small quark and gluon momenta, we have investigated the numerical implications of the newly derived theory ingredients on the B q → γγ decay phenomenology systematically in Section 6. The attractive advantage of implementing this nonperturbative LCDA model for the nonlocal NLP contributions consists in the complete independence of the yielding factorized results on the very profile function f (ω) parameterizing the leading-twist B q -meson distribution amplitude as presented in (134) In particular, the factorizable NLP corrections to the imaginary part of the helicity amplitude ImĀ L (B d → γγ) turned out to be even more pronounced: as large as (26 − 35) % numerically. Moreover, it has been demonstrated that the subleading power nonfactorizable contribution from the "hadronic" photon component develops the strong sensitivity on the precise shape of the leading twist B q -meson distribution amplitude.
Having at our disposal the desired expressions for the transversity amplitudesĀ and A ⊥ , we further provided the theory predictions for the CP-averaged branching fractions, the polarization fractions and the CP-violating observables for B q → γγ in the presence of the neutral-meson mixing. The major phenomenological results for confronting the forthcoming Belle II measurements have been collected in Tables 2 and 3 Bq , σ (2) Bq and the QCD renormalization scale ν. It is worth noticing that precision measurements of the ratio for the two branching fractions BR(B s → γγ) : BR(B d → γγ) have been proven to allow for determining the quantity λ B d : λ Bs with the estimated systematic uncertainties not worse than (5 − 10)% (see also [97] for an interesting calculation from the HQET sum rules). It has been also verified that the resulting uncertainties for the CP violation observables A Further developments of the precision QCD calculations for the double radiative B q -meson decays beyond the current work can be pursued forward in distinct directions. First, implementing the systematic QCD → SCET I → SCET II matching for the correlation functions (6) and (7) at NLP in the heavy quark mass expansion and constructing the appropriate SCET factorization formulae beyond the O(α 0 s ) accuracy will be in high demand for enhancing the predictive power of the perturbative factorization method and for promoting our understanding towards the emerged B-physics anomalies. To this end, it will be of utmost importance to develop the technical framework for tackling the soft-collinear convolution integrals in the presence of the rapidity divergences with effective filed theories (see [115,116] for further discussions in the different context). Second, improving the current theory constraints on the shape parameters of the twist-two HQET distribution amplitude will be indispensable for pinning down the yet considerable uncertainties particularly for the CP-averaged branching fractions. In this respect, it remains to be demonstrated evidently that whether the newly proposed method of accessing the nonperturbative light-front quantities with the corresponding equal-time correlation functions calculable on a Euclidean lattice and the standard OPE technique will be able to provide us with phenomenological encouraging results [117]. Third, extending the established strategies of evaluating the factorizable NLP corrections to the factorization analysis for a wide range of exclusive B-meson decays B → γ ν , B → ¯ ν [118] and B → M ν (with M = π, ρ, ω, K ( * ) , D ( * ) ) will be of high interest for probing the delicate strong interaction mechanisms governing a variety of heavy quark decays.

A The NLO coefficient functions F
In this appendix, we will collect the explicit expressions of the perturbative coefficients F (p) i,7 entering the NLO hard functions (19) by taking advantage of the obtained results of [17].

B Renormalization-group evolution functions
Here we will summarize the expanded expressions for the evolution functions in the SCET factorization formula (42) at the NLL accuracy in perturbation theory. As already explained in Section 3, the approximate expressions for the evolution kernelsÛ 1 (m b , µ h , µ) andÛ 2 (m b , µ h , µ) at O(α s ) can be derived from the expanded solution of U 1 (E γ , µ h , µ) as displayed in [52].
For completeness, we further collect the expansion coefficients in demand In addition, the newly defined function D(ξ, µ 0 , µ) is given by with the customary definition for the logarithmic derivative of the Gamma function