Form factors of B, Bs ->eta, eta' and D, Ds ->eta, eta' transitions from QCD light-cone sum rules

In the framework of the QCD light-cone sum rules (LCSRs) we present the analysis of all $B, B_{s}\to \eta^{(\prime)}$ and $D, D_{s}\to \eta^{(\prime)}$ form factors ($f^+, f^0$ and $f^T$) by including $m_{\eta^{(\prime)}}^2$ corrections in the leading (up to the twist-four) and next-to-leading order (up to the twist-three) in QCD, and two-gluon contributions to the form factors at the leading twist. The SU(3)-flavour breaking corrections and the axial anomaly contributions to the distribution amplitudes are also consistently taken into account. The complete results for the $f^0$ and $f^T$ form factors of $B,B_s \to \eta^{(\prime)}$ and $D, D_{s} \to \eta^{(\prime)}$ relevant for processes like $B \to \eta^{(\prime)} \tau \nu_{\tau}$ or $B_{s} \to \eta^{(\prime)} l^+ l^-$ are given for the first time, as well as the two-gluon contribution to the tensor form factors. The values obtained for the $f^+$ form factors are as follows: $f^+_{B\eta}(0)= 0.168^{+0.042}_{-0.047}$, $|f^+_{B_s\eta}(0)|= 0.212^{+0.015}_{-0.013}$, $f^+_{B\eta^\prime}(0)= 0.130^{+0.036}_{-0.032}$, $f^+_{B_s\eta^\prime}(0)= 0.252^{+0.023}_{-0.020}$ and $f^+_{D\eta}(0)= 0.429^{+0.165}_{-0.141}$, $|f^+_{D_s\eta}(0)|= 0.495^{+0.030}_{-0.029}$, $f^+_{D\eta^\prime}(0)= 0.292^{+0.113}_{-0.104}$, $f^+_{D_s\eta^\prime}(0)= 0.558^{+0.047}_{-0.045}$. Also phenomenological predictions for semileptonic $B, B_{s}\to \eta^{(\prime)}$ and $D, D_{s}\to \eta^{(\prime)}$ decay modes are given.


Introduction
In the view of the numerous precise new measurements of two-body nonleptonic and semileptonic B, B s and D, D s decays to η ( ) performed by BaBar and Belle recently [1] and the upcoming experimental precision in the next-generation experiments it is timely to provide precise predictions for B, B s → η ( ) and D, D s → η ( ) form factors for analysis of these decays. The form factors parametrize hadronic matrix elements of quark currents and describe the long-distance QCD effects in semi-leptonic and non-leptonic decays.
All those decays are important for testing and understanding the Standard Model flavour interactions, in particular for our understanding of the QCD dynamics in the flavour physics as well as the flavour mixing given by the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix. The B, B s and D, D s decays to η, η pseudoscalar mesons can be used to shed some light on both of these phenomena. paper.

Mixing
To analyze η and η states, we have to deal with several definitions of matrix elements of the flavour-diagonal axial vector and pseudoscalar current: where q = u, d and the isospin limit is taken, m q = 1 2 (m u + m d ). There is also a U (1) A anomaly, which is connected with derivatives of the currents through the equation of motion as and included in h q,s P as In the exact SU(3) flavour-symmetry limit a P = 0. It is known that the SU(3) breaking corrections for η and η are large and that η and η mix since they are not pure a flavour-octet and a flavour-singlet states, respectively.
The mixing of η and η mesons is established in two mixing schemes: the singlet-octet (SO) and the quark-flavour (QF) scheme. Each of the schemes has some advantages and some disadvantages.
In the SO scheme the mixing occurs among SU (3) F singlet |η 1 = 1/ √ 3|uū + dd + ss and octet |η 8 = 1/ √ 6|uū + dd − 2ss components. By defining the coupling of the axialcurrents to η and η mesons as Since only singlet component mixes with the gluonic contributions, the renormalization scale dependence of parameters is diagonalized in the SO scheme and therefore is suitable for the analysis of the gluon distribution amplitudes [14]. Moreover, f 8 is scale independent and f 1 renormalizes multiplicatively: where µ 0 = 1 GeV, the scale at which the values of the mixing parameters are determined [15]. The simpler mixing scheme is QF scheme. There the basic components are |η q = 1/ √ 2|uū + dd and |η s = |ss states and the decay constants are defined as 0|J r 5ν |η ( ) (p) = if r η ( ) p µ , (r = q, s) . (2.8) Their mixing with the decay constants of pure (hypothetical) non-strange and strange states, f q and f s respectively, is given by The main advantage of this scheme is that the mixing is not governed by the (large, 10 − 20%) SU (3) F breaking effects as in the SO scheme, but by the OZI-rule violating contributions which have be proven to be small [15]. Therefore it is possible to parametrize the mixing just with one angle φ and the matrix U (φ) given as which leads to the following expressions The parameters have been determined by fits in [15] as and will be also used in this paper 1 . These values give for the parameters of the SO basis the following: (2.14) Due to the mixing of the flavour-singlet and gluonic components, in the QF scheme both η q and η s will get gluonic contributions and therefore also the physical η and η states. The flavour states in QF scheme and in the approximation above can be written as [14] |η q = f q 2 √ 2N c ψ q (x, µ F )|qq + 2/3ψ g (x, µ F )|gg (2.15) 2N c ψ s (x, µ F )|ss + 1/3ψ g (x, µ F )|gg (2.16) where |qq = (uū + dd)/ √ 2 and ψ q = 1/3(ψ 8 + 2ψ 1 ) and ψ s = 1/3(2ψ 8 + ψ 1 ). By combing above information about the nature of η and η states one can expect that gluonic contributions |gg will be larger for η mesons, which is confirmed by the final results.
Until now there is no available QCD sum rule or lattice QCD calculations of B s to η ( ) transition form factors f +,0,T Bsη ( ) . Since these transitions probe only the |ss content, one can use the approximation in the quark flavour scheme which neglects completely the gluonic contribution. The calculation presented in this paper will check for the SU (3) F breaking effects in the above relations.

Distribution amplitudes
The light-cone distribution amplitudes (DA), giving the momentum fraction distribution of valence quarks of η and η are defined analogously to other meson light-cone DAs, by expanding the non-local operators on the light-cone in terms of increasing twist, but paying attention to the specific flavour structure of η ( ) mesons. The twist 2 two-quark DAs φ i 2,P of P = η ( ) mesons are defined as where as usual z µ is the light-like vector and [z, −z] is the path-ordered gauge connection and u is the momentum fraction of a valence quark. In the SO basis one will have C 1 = 1/ √ 3 and The twist-2 two-quark DAs of η ( ) are symmetric in their argument and therefore they can be expanded in terms of Gegenbauer polynomials as usually: The coefficients a P,i n are the Gegenbauer moments of the quark DA. The gluonic twist-2 DA φ g 2,P of P = η ( ) mesons are defined by the following matrix element (for detailed discussion on the derivation of gluonic DA and its mixing with the quark states in mesons see for example [19]): It is antisymmetric and therefore and it is expanded in terms of C 5/2 n Gegenbauer polynomials where the coefficients b P,g n are the Gegenbauer moments of the gluon DA and we take b η,g n = b η ,g n and keep only the first term in the sum, n = 2. Although b η,g n and b η ,g n could differ, this approximation is justified since their values are subject of large uncertainties.
In the calculation we use the following matrix element of the η ( ) over two gluon fields . (2.23) With the above normalization of the DA, the renormalization mixing of twist-2 quark and gluonic distribution amplitudes is given as and it is numerically small. But, the mixing is important for p 2 = 0 case, since it verifies the collinear 'factorization formula' for the form factors 25) and proves that the separation of the transition form factors in perturbatively calculable hard-scattering T H part and a nonperturbative DA is essentially independent on the factorization scale µ IR [20]. This is an essential step of calculation which is going to be proved for each of the F -correlation function at the order of twist 2, see discussion in the next section. The explicit solutions of (2.24) can be find in [13] and in the appendix B of [12]. In the asymptotic case, when Q 2 = −q 2 → ∞ the twist-2 quark and gluon DAs evolve to their asymptotic forms In that case, there is no gluonic contribution at the twist-2 level to the form factors, and the residual µ IR dependence in the twist-2 NLO quark contribution integrates with φ i 2,P (u) |asym to zero, which again confirms the µ IR independence of the complete result.
To include SU(3) flavour-breaking corrections consistently we keep not only m 2 η ( ) corrections and quark masses in the hard-scattering amplitudes, but also in the distribution amplitudes. Therefore we do not use the approximations in the twist-3 and twist-4 contributions employed in the literature where the following replacements are used in DAs: for M → η q and M → η s decays respectively. Instead we are going to use (in the QF scheme): Although the above quantities, especially h q , are weakly constrained due to the numerical cancellations, we use them for the consistency of our calculation. Actually, we will see later that the approximation in (2.27) for h q is quite bad and causes somewhat large values of form factors of D, B → η ( ) . Distribution amplitudes of higher twist are defined following [12] and [21]. Their parameter evolutions and definitions include now the anomaly contribution a P with the following expressions [22]: (2.30) Therefore in [12] the normalizations of two-particle twist-3 DAs φ p,σ 3 differ from those in [21]. In [12] one can find a consistent treatment of m s corrections up to twist-4 and of anomalous contributions to DA and we take definitions and expressions given there. Then, where r = q, s and 2m r 0|r(z 2 n)σ µν γ 5 r(z 1 n)|P (p) = i(z 1 − z 2 ) 6 (p µ n ν − p ν n µ ) (2.32) The normalization is then and By calculating the mixing of twist-4 DAs, some approximations in the twist-3 DA are made in [12] when compared to the expressions in [21], to keep the same order of calculation in the conformal spin and the quark masses.
The expressions for the two-particle twist-3 DAs used (contributions of higher conformal spin and O(m 2 s ) corrections are neglected; see also [21], Eqs. (3.25-26)) are (2.36) The three-particle quark-gluon-antiquark DA is defined as usual [21] 0|r(z)σ µν γ 5 gG αβ (vz)r(−z) There are two two-particle twist-4 DAs ψ      4P (α). All details and subtleties in derivation of these improved twist-4 DAs with the corrected mass corrections and inclusion of the anomalous contribution can be found in Appendix A of [12]. Here we just quote the expressions: and Φ (r) (2.41) The parameters which appear here are parametrization of various local matrix elements and their values are taken from [21] and listed in Appendix B. The above twist-4 expressions are valid for flavour-octet contributions where there is no mixing with the gluonic twist-4 DA. For the flavour-singlet case one has to take this mixing into account. In the approximation taken in [12] the twist-4 DAs are given by the replacement everywhere at the twist-4 level where the mass m 2 P occurs. As it was discussed in [12] this substitution ensures for the given accuracy the consistent normalization of the twist-3 and twist-4 DA and ensures that the same mixing FKS scheme applies also for the higher-twist contributions.
For the values of parameters involved we will use crude estimates in terms of the pion and kaon DA parameters derived from the sum rules [21], see Appendix: while the corresponding η, η parameters will be given through the mixing as   Figure 1. Diagrams corresponding to the leading-order terms in the hard-scattering amplitudes involving the two-particle (left) and three-particle (right) η ( ) DA's shown by ovals. Solid, curly and wave lines represent quarks, gluons, and external currents, respectively. For B s → η ( ) transition, u is replaced by s. In the case of D → η ( ) transitions, u → d and b → c and correspondingly d is exchanged by s for D s → η ( ) . u → s in (3.1) and j Bs = m bb iγ 5 s interpolating current. Again, D s case is then obtained trivially by replacing b-quark with the c-quark.
Since we want to explore also the SU(3) symmetry breaking, we will keep the η ( ) masses (p 2 = m 2 η ( ) ) in (3.1). The light quark masses will be systematically neglected, except when they occur in ratios in the distribution amplitudes.
The method of the LCSR is very well know and we will here just briefly outline the procedure in order to properly define all ingredients necessary for calculating the form factors. For the large virtualities of the currents above, the correlation function is dominated by the distances x 2 = 0 near the light-cone, and factorizes to the convolution of the nonperturbative, universal part (the light cone distribution amplitude (DA)) and the perturbative, short-distance part, the hard scattering amplitude, as a sum of contributions of increasing twist.
We calculate here contributions up to the twist-4 in the leading order, O(α 0 s ), and up to the twist-3 in NLO, neglecting the three-particle contributions at this level. Schematically, the contributions are shown in Fig.1 and Fig.2. Due to the specific properties of η and η Figure 3. Diagrams contributing to the gluonic hard-scattering amplitudes at O(α s ).The first diagram is IR divergent and its divergence will be absorbed by the evolution of the gluon DA. See text.
mesons discussed above, there are additional gluonic diagrams contributing to M → η ( ) form factors shown in Fig.3. These contribution has only been calculated for f + Bη ( ) form factor at twist-2 level in [13] and for m 2 η ( ) = 0. Here we are going to calculate these contributions for other form factors f 0 M η ( ) and f T M η ( ) by neglecting O(α s m η ( ) ) effects in both DAs and the hard-scattering part. This approximation is justified having in mind that parameters of DA for the gluonic DA of η and η are badly known, see the values of b η ( ) ,g 2 parameter below. By using hadronic dispersion relation in the virtuality (p + q) 2 of the current in the B channel, we can relate the correlation function (3.1) to the B → η ( ) matrix elements, and extract the form factors. In the literature it sometimes appears that the form factors are defined as above by divided by a factor √ 2 to match the transition form factors of η, η with those of a pion when there is no η − η mixing and in the limit of the conserved SU(3)-flavour symmetry [13].
Inserting hadronic states with the B-meson quantum numbers between the currents in (3.1), and isolating the ground-state B-meson contributions for all three invariant amplitudes F (q 2 , (p + q) 2 ), F (q 2 , (p + q) 2 ) and F T (q 2 , (p + q) 2 ) and using (3.2) and (3.3) obtains: The scalar B → η ( ) form factor is then a combination of the vector form factor (3.4) and the form factor from (3.5), and is only present in the semileptonic B (s) , D (s) → η ( ) lν decays when the lepton mass is not neglected and the rare B s , D s → η ( ) l + l − decays. In above, F 0(1) and F 0(1) represent the LO (NLO) contributions and are leading order twist-2 two-gluon contributions calculated explicitly in the paper. At the leading twist-2 level there is no gluonic contribution in (3.5). However, note from (3.7) that this does not mean that twist-2 two-gluon contributions will not appear in the scalar f 0 M η ( ) form factors (3.7). As usual, the quark-hadron duality is used to approximate heavier state contribution by introducing the effective threshold parameter s B 0 and the ground state contribution of B meson is enhanced by the Borel-transformation in the variable (p + q) 2 → M 2 . Completely analogous relations are valid for B s → η ( ) form factors, with the replacement u → s in ). In addition, in the derivation of above expressions for B s , one has to take into account that B s |biγ 5 s|0 /m 2 The same is valid for D, D s form factors with the replacement m b → m c and the appropriate exchanges described before.
The calculation will be performed in M S scheme. The B, B s and D, D s decay constants will be calculated in the M S scheme using the sum rule expressions from [23] with O(α s , m 2 s ) accuracy. In that way we achieve the consistency of the calculation and the cancellation of uncertainties in the sum rule parameters.
Each form factor can be written in the form of the dispersion relation: The leading order parts of the LCSR for Up to now, SU(3)-violating effects for f D (s) η ( ) , f B (s) η ( ) form factors were not systematically studied, since the effects of inclusion of m 2 η ( ) effects complicate the calculation, especially at NLO in the hard-scattering amplitudes. However, while the complete SU(3)symmetry breaking corrections in η ( ) DAs of twist-3 and twist-4 are now known [12], it is worth to have a consistent picture of all SU(3)-breaking corrections and we will include complete SU(3)-breaking effects in both DAs, as well as in the hard scattering amplitudes at LO. At NLO in the hard-scattering amplitudes, for the cases when the mass of a light quark cannot be neglected, as for m s , the inclusion of m s and m 2 η ( ) effects complicate the calculation. As already known from the analysis of B (s) → K from factors done in [11], inclusion of quark mass effects leads to the mixing between different twists and the fully consistent calculation with m s included in the quark propagators is not possible, see discussion in [11]. However, here we have η as a finite-state particle which mass is much larger than m s and therefore, in the NLO quark and gluonic amplitudes we set m s = 0 and p 2 = m 2 η = 0. Each form factor can be expressed as and and explicitly where Fq q 0 and Fs s 0 (Fq q 1 and Fs s 1 ) are LO (NLO) contributions from quark hard-scattering amplitudes for each of the form factors and F gg 1 is the NLO gluonic contribution proportional to the singlet-flavour decay constants The f (r) η ( ) decay constants are given in (2.12). Analogous expressions are valid for D (s) → η ( ) decays.
Obviously, for B, D → η ( ) transitions the main contribution comes from η q meson states and η s contributes only through suppressed gluonic contributions, while for B s , D s → η ( ) transitions the leading η s meson state contribution will receive, through the gluonic diagrams, a small mixture with η q state. Also, implicitly there will be mixing with among twist-2 quark and gluonic distribution amplitudes Eq.(2.24), which will bring b η ( ),g 2 dependence in the twist-2 quark LO (Fq q 0 and Fs s 0 ) and NLO contributions (Fq q 1 and Fs s 1 ) and a η ( ),1 2 dependence to the gluonic contributions F gg 1 . Since η and η are mixtures of the |qq , |ss in the calculation of the quark contributions we will use (with appropriate substitutions) our NLO results for the hard-scattering part for B → π [10] and B (s) → K form factors [11] with the p 2 effects included at the LO (up to twist-4) and NLO level (up to twist-3) and will imply recently derived DAs of η and η with the SU(3)-breaking effects and the axial anomaly contributions included up to twist-4. The gluonic contributions, which are already NLO effect, will be calculated for p 2 = m 2 η ( ) = 0.

LCSR for gluonic contributions to the form factors and consistent treatment of the IR divergences appearing
The gluonic contributions at the O(α s ) to the B(D) → η ( ) and B s (D s ) → η ( ) form factors come from the diagrams in Fig.3. The results for the form factors f +,0,T M η ( ) are presented in subsection 4.1. They are added to the quark contributions (3.11) to get the complete result at the order O(α s ).
The first diagram is Fig.3 is IR (collinear) divergent. This divergence has to disappear for the general collinear factorization formula used here be valid. The scale µ F is the factorization scale. At the twist n = 2 level, as already mentioned, there will be mixing of quark and gluonic contributions and the hard-scattering (perturbative part) T H and the distribution amplitude Φ 2 can be represented as In order to consistently treat this mixing we have to examine the evolution of the DAs, at the same O(α s ) as the calculation of the perturbative part T H . Due to the mixing the standard Brodsky-Lepage (BL) evolution equation [24] µ will be a matrix equation now, where V (u, v, µ F ) is the perturbatively calculable evolution kernel with the LO kernel of our interest of the form and V ij are well-know evolution kernels [14,19,25] which we cite here for convenience: (4.6) These evolution kernels are exactly those which govern the renormalization of the DAs The connection between Z and the evolution kernel V is given as (D = 4 − 2 ). On the other hand, by calculating the hard-scattering part T H , owing to the fact that final-state quarks are taken to be massless and on-shell (for the case p 2 = 0), the amplitude contains collinear singularities. Since T H is a finite quantity by definition, collinear singularities have to be subtracted. Therefore, T factorizes as with collinear singularities being subtracted at the scale µ F and absorbed into the constant Z T,col . As usual The UV singularities are removed by the renormalization of the fields and by the coupling-constant renormalization at the (renormalization) scale µ R . Now, in order that the factorization formula is valid, the following has to be satisfied The divergences of T (u, Q 2 ) and Φ(u) in (4.1) then cancel and at the end we are left with the finite perturbative expressions for all form transition factors It is worth pointing out that the scale µ F representing the boundary between the low-and high-energy parts in (4.1) plays the role of the separation scale for collinear singularities in T (u, Q 2 ), on the one hand, and of the renormalization scale for UV singularities appearing in the perturbatively calculable part of the distribution amplitude Φ(u), on the other hand.
The general discussion and all details of the proof of the cancellation of the factorization scale dependence in the collinear factorization formula (4.1) at all orders of calculation can be found in [20,26]. In our case of calculating the heavy-to-light transition form factors f +,0,T we face the following situation. The hard-scattering, perturbatively calculable pieces coming from the diagrams from Fig.2 have UV and infra-red singularities at O(α s ). We have already proven in [10,11,27] for B → π and B (s) → K form factors that the IR divergences of the quark contributions at twist 2 level cancel exactly with those coming from the evolution kernel V qq . Here, due to the mixing with the twist 2 gluonic contributions, the convolution of V qg of the T H LO will exactly cancel the IR divergence in the first gluonic diagram in Fig.3. At the twist 3 level of O(α s ) the IR divergences of quark diagrams mutually cancel, as shown before in [10,11]. This gives the final proof of the collinear factorization formula at the given order for the heavy-to-light M → η transition form factors. In the calculation of the gluonic contributions to the form factors we have faced the problem of the consistent treatment of the γ 5 in the dimensional regularization. Leading order for the gluonic amplitude is given by one-loop Feynman diagrams in Fig.3. and we have to deal with IR divergence which is a consequence of having massless quarks propagating through the loops. In the calculation of the gluonic contributions to the form factors it appears a Levi-Civita tensor in the projector of the twist-2 two-gluon DA (2.23) and a single γ 5 matrix in the trace which are both quantities with well-defined properties only in D = 4 space-time dimensions. Generalization of these quantities in D dimensions is problematic and different approaches to avoid resulting ambiguities can be found in the literature. Moreover, in our case there is no gluonic contributions which appear at LO on α s that would greatly help in resolving the γ 5 problem at NLO level. The problem was not addressed in the paper where the gluonic amplitude was evaluated for the first time [13] and it is not clear how they resolved the ambiguities.
In the case of the interest it is possible to completely avoid γ 5 problem and all connected complications since the IR divergence is direct consequence of the massless quark lines and putting a small mass m in massless quark propagators regularizes (removes) the divergence. As a consequence, we are not forced to use dimensional regularization and calculation can be performed in four dimensions without any problem. Note that putting mass in quark propagators doesn't spoil any of properties and symmetries of the amplitude contrary to the case when, so called, mass regularization is used on gluon propagators. At the very end of calculation it is necessary to expand final result around zero for the small introduced quark mass m. The IR divergence will now reappear as ln(m 2 ) term and it is straightforward to connect it with 1/(D-4) term in the framework of dimensional regularization.
The obtained expressions are as follows: where the gluonic contribution to f + form factors is − 37m 6 − m 4 56q 2 + 55s + m 2 18q 4 + 76q 2 s + 17s 2 + 3q 6 − 27q 4 s − 11q 2 s 2 − 2s 3 (4.12) and the corresponding contribution to f T form factors has the following form − 59m 6 − m 4 (72q 2 + 85s) + m 2 s(84q 2 + 23s) + 3(6q 6 + 6q 4 s + 6q 2 s 2 − s 3 ) The result for the gluonic O(α s ) contribution to f + form factors was first given in the appendix of [13]. Our result (4.12) does not completely agree with the one presented there. While we agree in the part being proportional to the logarithmic terms, there is a disagreement between the coefficients in the second line of (4.12) and the expression (A.1) from [13]. Since those terms are exactly those which change with the different treatment of γ 5 , and the authors of [13] have not placed any comment how they have resolved the γ 5 ambiguities in the calculation of f gg,+ , we assume that the difference comes from the improper treatment of the γ 5 in [13].
The result for the gluonic O(α s ) contribution to f T , Eq.  Table 4: where the quoted error intervals are coming from the variation of s M 0 and M 2 M only since other uncertainties are canceled in ratios in Eqs.(3.5,3.6,3.7). By comparing our results with the previous LCSR results and the most recent determinations from [28], where in the perturbative part the higher order corrections were included, we see good agreement. The results are also within uncertainties of the lattice QCD calculations of the same decay constants [29].
For the f D and f Ds the experiment gives somewhat lower values [1], f D = 204.6 ± 5.0 MeV, f Ds = 257.5 ± 4.6 MeV , but still consistent within the complete LCSR error results [28].
The renormalization scale is given by the expression µ B ( s) = m 2 B (s) − m 2 b and similarly for D (s) → η ( ) transitions. Therefore, for the renormalization scale we use µ = 3 GeV, for the f 0,+,T Bη ( ) form factors and µ s = 3.4 GeV for f 0,+,T Bsη ( ) and for µ D = 1.4 GeV and µ Ds = 1.5 GeV. As usual, we will check the sensitivity of the results on the variation of above scales and will include it in the error estimation.
The method of extraction of the Borel parameters M and the effective thresholds s 0 for f +,0,T M η ( ) form factors is the same as described in [10]. It relies on the requirement that the derivative over −1/M 2 of the expression of the complete LCSR for a particular form factor, which gives heavy-meson masses m 2 M , does not deviate more than 0.5 − 2.5% from the experimental values for those masses. Additional requirements such as that the subleading twist-4 terms in the LO, are small, less than 10% of the LO twist-2 term, that the NLO corrections of twist-2 and twist-3 parts are not exceeding 30% of their LO counterparts, and that the subtracted continuum remains small, are also satisfied. These demands provide us the central values for the LCSR parameters listed in Tab = ±20, which dependence is explicitly displayed in the errors. The errors are compilation of the variation of parameters added in quadratures. In the errors we explicitly stress SR parameter dependence (s 0 , M ), η − η mixing parameter dependence (mix) and dependences coming from the variation of the rest of parameters (rest={µ, m c,b , a 2 , a 4 }).
The errors of the results are much larger for the transitions B, D → η ( ) where B, D → η q dominates then for B s , D s → η s decays since the error in the parameter h q (2.29) is huge, of O(200%) depending not on η − η mixing parameters but exhibiting a numerical cancellation among terms. If one would use approximation (2.27) applied in [30] instead, the (rest)-errors would be almost an order of magnitude lower and the mean values would be somewhat larger for those decays, which we assume is the main reason, apart from the rest of SU (3) F approximations used there, of the discrepancies with some of the results presented in [30]. We see that the dominant errors in M → η form factors is coming from the variation of b η ( ),g 2 and it amounts to about 15%, while in M → η decays come to 2%. Our findings for calculated B → η ( ) form factors agree very well with those from [13,[31][32][33].
Their q-dependence of the form factors and their ratios is shown in Fig.4-9. From Fig.5 and Fig.8 we see that the gluonic corrections are much larger for B (s) , D (s) → η decays then for M → η, as expected. Also the gluonic corrections are larger in D (s) decays. It is obvious that even in ratios of form factors the gluonic contributions give main error and that it would be difficult to constrain b 2 , unless all M → η ( ) semileptonic transitions are measured, Fig. 6 and 9.
We can now investigate SU(3) F approximations from (2.17). By using the obtained results and the result for f BK from [11] we obtain  We can note that the approximation works quite well although somewhat better for M → η decays than for M → η transitions. There exists LCSR calculations of f + Dsη form factor [34,35]. In these papers the f + Dsη form factor is then obtained by using the relation While their predictions for f + Dsη agree with ours, the use of the above approximative relation which neglects the gluonic contributions gives somewhat larger f + Dsη form factor then the one obtained here, (5.2,5.4).
There exist also recent lattice results on D s → η ( ) form factors [36]. These transitions at the lattice are challenging due to the presence of disconnected quark-line contributions and in [36] only the scalar f 0 Dsη ( ) form factors are calculated, which at q 2 = 0 are equal to the f + . By comparing the results one can see that the lattice predictions give f + Dsη < f + Dsη , which is just opposite in LCSR for all M s → η ( ) transitions. The tendency f + M η < f + M η in LCSR is established for non-strange meson decays, see results in (5.2,5.4).

Phenomenological applications
In this section we comment on some phenomenological results for semileptonic D (s) → η ( ) and B (s) → η ( ) decays which include the calculated from factors. To be able to calculate the branching ratio we need the form factor extracted in whole accessible kinematical regions. For D (s) decays the LCSR are applicable only in the region q 2 m 2 c and for B (s) the region is 0 < q 2 < 12 GeV.
The are many parametrization for calculating the shape of form factors at q 2 = 0. All of them work equally well and therefore we decided to use the most simplest one [37]: where the extrapolation of the form factors is performed just by fitting one parameter α i for each of the decays and using the appropriate vector meson resonances m * H , Tab.6, while the normalization is given by the form factors at q 2 = 0. The fitted parameters α i for D (s) The semileptonic D (s) → η ( ) eν e and B → η ( ) eν e decay rates are calculated by  For the rare B s → η ( ) l + l − (νν) decays we use the effective Standard Model hamiltonian for b → sl + l − (νν) transitions [38] and calculate decay rates as [39] Γ where [39]. For the Wilson coefficients we use the following values C 7 = −0.3031 , C 9 = 4.1696 , C 10 = −4.4641 , C L = 2.74 · 10 −9 . (6.7) Our predicted branching ratios for various M → η ( ) decays are given in Tab.1. By comparing with the existing calculations performed in different models [39,[44][45][46] we agree quite well, expect that we predict somewhat larger branching ratios for B s → η ( ) τ + τ − decays.
Because of the larger errors in B, D → η ( ) decays, M s → η ( ) would be better for extraction of the unknown b η ( ) ,g 2 parameter, but measurements of M s decays still have to achieve sufficient precision, in particular Br(B s → η l + l − ) and Br(B s → ηl + l − ) are challenging with the branching ratio of O(10 −7 − 10 −8 ) but they could be measured at future SuperB and SuperKEK experiments.

Summary
We have investigated B, B s → η ( ) and D, D s → η ( ) form factors (f + , f 0 and f T ) by including m 2 η ( ) corrections in the leading (up to the twist-four) and next-to-leading order (up to the twist-three) in QCD, as well as gluonic contributions to the form factors at the leading twist in the framework of the QCD light-cone sum rules and have also taken SU(3)-flavour breaking corrections and the axial anomaly contributions to the distribution amplitudes consistently into account. The two-gluon twist-2 contributions are calculated for all f + , f 0 and f T form factors.
We have given the values and shapes at q 2 = 0 of all calculated form factors and have shown predicted ratios for some semileptonic B, B s → η ( ) and D, D s → η ( ) decay modes.
With the determined form factors of transitions B, B s → η ( ) it will be possible to analyze consistently nonleptonic decays to charmonia and to test the factorization hypothesis in such transitions which we be a subject of the future investigations.

Acknowledgments
We are grateful to D. Bečirević and K. Passek-Kumerički for useful discussions. The work is supported by the Croatian Science Foundation (HrZZ) project "Physics of Standard Model and Beyond", HrZZ 5169 and by the Munich Institute for Astro-and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe".
A Explicit results for f + , f 0 and f T form factors at the leading order in B, B s → η ( ) and D, D s → η ( ) transitions The leading O(α s ) part of the f + B (s) η ( ) LCSR, (3.4), has the following expression (P = η, η ; r = q for B q → P and r = s for B s → P ; for D, D s the same expressions are valid with the replacement m b → m c ): = cos φf s , and similarly for the two-particle twist-three DAs φ p,σ 3P : φ (q)p,σ 3η In the case of the twist-2 DA, we will express the decay constants F at m 0 = 1 GeV, the energy at which the FKS parameters are determined and Numerically, The short-hand notations introduced for the integrals over three-particle DA's are 2 : where Finally, the leading order LCSR for the penguin form factor, (3.6), reads: The expressions for f +,0,T D (s) from factors follows from above, by replacing m b by m c .

B Parameters used in the calculation
In this appendix we summarize the parameters used in the calculation of f M η ( ) form factors as well as in the calculation of f M decay constants, Tables 2-5. In Table 6. we summarize meson masses, lifetimes and vector resonances used in the calculation of phenomenological predictions for semileptonic M → η ( ) decays.   1.16 ± 0.05 Table 4. Decay constants used in the paper, the values are in MeV. The decay constants of heavy mesons are obtained from the two-point SR at O(α s ) and agree with those from [10,11,28]. The quoted errors are coming only from the variation of s 0 and the Borel parameter M , since other errors will cancel in the calculation of the form factors. For comparison the recent more complete LCSR results from [28] are given.   Table 6. Meson masses and lifetimes. The vector meson resonances m * H are used in the extrapolation formula for q 2 -dependence of the form factors (6.1).