Exclusive Radiative Z-Boson Decays to Mesons with Flavor-Singlet Components

We present a detailed study of the exclusive radiative decays $Z\to\eta^{(\prime)}\gamma$ employing the QCD factorization approach. We derive a factorization formula for the decay amplitudes valid at leading power in an expansion in $(\Lambda_{QCD}/m_Z)^2$, which includes convolutions of calculable hard-scattering kernels with the leading-twist quark and gluon light-cone distribution amplitudes of the mesons. Large logarithms arising in the evolution from the high scale $m_Z$ down to hadronic scales are resummed using the renormalization group, carefully accounting for the effects of the heavy bottom and charm quarks. Our results for the branching ratios are very sensitive to hadronic input parameters, such as the decay constants and mixing angle characterizing the $\eta-\eta'$ system. Using the most recent estimates of these parameters, we obtain the branching ratios $Br(Z\to\eta\gamma)\sim 1.6\cdot 10^{-10}$ and $Br(Z\to\eta'\gamma)\sim 4.7\cdot 10^{-9}$. A measurement of these processes at a future high-luminosity Z factory could provide interesting information on the gluon distribution amplitude.


Introduction
Exclusive decay processes involving individual hadrons in the final state pose a formidable challenge to theoretical physics, because the complicated strong-interaction physics describing hadronic bound states cannot be described using perturbative methods. For the case of hard exclusive processes the QCD factorization approach [1][2][3] provides a systematic framework for factorizing calculable short-distance effects associated with high energy scales from hadronic dynamics, which is described in terms of light-cone distribution amplitudes (LCDAs) of individual hadrons. This approach has been applied successfully for different processes, such as meson form factors at large momentum transfer (see [4,5] for recent discussions) and hadronic weak decays of heavy B mesons [6]. In previous work, we have studied the exclusive radiative decays of Z → M γ and W → M γ into final states containing a single meson M as an ideal testing ground for the QCD factorization approach, arguing that power corrections to the factorized amplitudes are suppressed by (Λ QCD /m Z,W ) 2 and are thus bound to be very small [7]. By including higher-order QCD corrections in the short-distance coefficients and solving their renormalization-group (RG) evolution equations, large logarithms of the form α s ln(m 2 Z /µ 2 0 ) n , where µ 0 ≈ 1 GeV is a typical hadronic scale, can be resummed to all orders of perturbation theory. Applying the same formalism to the exclusive radiative Higgs-boson decays h → V γ, where V is a vector meson, provides access to the Yukawa couplings of the Higgs boson to light quark flavors and thus serves as a powerful probe of physics beyond the Standard Model [7][8][9][10]. When studying the decays Z → M γ in [7] one set of processes was left out, namely those where the final state pseudoscalar meson M = P has a flavor-singlet component in its wave function. An important complication in this case lies in the fact that there exists a new contribution to the decay amplitude at leading order in power counting, where the meson is formed from two collinear gluons instead of a quark-antiquark pair. The existence of this contribution not only gives rise to a more complicated form of the factorization formula but also influences the RG equations satisfied by the short-distance coefficients [11][12][13][14]. In this paper we present a detailed analysis of the decays Z → η ( ) γ in the context of QCD factorization, treating flavor mixing in the Feldmann-Kroll-Stech (FKS) scheme [15] and carefully accounting for the decoupling of the heavy bottom and charm quarks in the evolution from the high-energy scale m Z down to low energies.

Theoretical Framework
In previous work the QCD factorization formula for exclusive radiative decays Z → M γ was derived for the case of flavor-nonsinglet pseudoscalar or vector mesons M = P or V , which are produced via a quark-antiquark pair [7]. Representative Feynman diagrams contributing at leading and next-to-leading order (NLO) are shown in Figure 1. For a pseudoscalar meson in the final state, the decay amplitude can be written in the general form where k and q denote the meson and photon momenta, e and g are the electromagnetic and weak coupling constants, and θ W is the weak mixing angle. A second form factor, which is allowed by Lorentz invariance, vanishes since the final-state meson P is an eigenstate of the charge-conjugation operator. At leading order in an expansion in powers of (Λ QCD /m Z ) 2 , the form factor for the case of Z → π 0 γ decays reads [7] where f π 0 ≈ 130 MeV is the decay constant of the pion and φ π 0 (x, µ) is its leading-twist LCDA. Q q and v q = 1 2 T q 3 − Q q sin 2 θ W are the electric and weak vector charges of a quark with flavor q. At NLO in QCD perturbation theory, the hard-scattering kernel H q (x, µ) reads where C F = 4/3 is a color factor, and [16] h q (x, µ) = (2 ln x + 3) ln m 2 The convolution integral of the hard-scattering kernel with the LCDA in (2) is independent of the choice of the factorization scale µ. The Z → π 0 γ branching ratio turns out to be strongly suppressed, Br(Z → π 0 γ) = (9.80 ± 1.03) · 10 −12 [7], because the relevant combination of quark charges is numerically very small. We use sin 2 θ W = 0.23126(5) as extracted from the neutral-current couplings of the Z boson [17]. It is thus more promising to search for the decays Z → η ( ) γ. Two complications arise in this case. First, the physical η and η mesons are complicated mixtures of quark-antiquark states with different flavor. Second, and more profoundly, the flavor-singlet quark-antiquark state mixes with a pure gluon state under renormalization, and indeed the η and η mesons can also be produced via a two-gluon LCDA.

Factorization in the presence of flavor-singlet contributions
The Z → P γ decay amplitudes, where P = π 0 , η, η , . . . denotes a neutral pseudoscalar meson, can be calculated from first principles using the QCD factorization approach [1][2][3], because the energy E released to the final-state meson is much larger than the scale of long-distance hadronic physics. At leading power in an expansion in λ = Λ QCD /m Z , the form factors F P can be written as convolutions of calculable hard-scattering coefficients with leading-twist LCDAs of the meson P . To derive the corresponding factorization formula we employ the formalism of soft-collinear effective theory (SCET) [18,19], which provides a systematic expansion of decay amplitudes in powers of λ. For the purposes of this discussion we work in the rest frame of the decaying Z boson and assign momenta k µ = En µ and q µ = En µ to the meson and photon, respectively, where E = m Z /2, while n µ = (1, 0, 0, 1) andn µ = (1, 0, 0, −1) are two light-like vectors. Up to power corrections of order (m P /m Z ) 2 the meson mass can be set to zero. The light final-state meson moving along the direction n µ can be described in terms of collinear quark, antiquark and gluon fields. These particles carry collinear momenta p c that are approximately aligned with the direction n. Their components scale like (n·p c ,n·p c , p ⊥ c ) ∼ E(λ 2 , 1, λ). Note that p 2 c ∼ Λ 2 QCD , as appropriate for an exclusive hadronic state. The collinear quark and gluon fields are introduced as gauge-invariant objects dressed with Wilson lines. Explicitly, one defines [20,21] where iD µ c = i∂ µ +g s A µ c denotes the covariant collinear derivative, and W c is a collinear Wilson line extending from the location of the field to infinity along the directionn. Both fields are of O(λ) in SCET power counting, while other components of the gluon field are of higher order. Adding more component fields to an operator thus always leads to further power suppression in λ. At leading order, the operators with a non-zero matrix element between the vacuum and a single meson state are thus bilinears of the formX c . . . X c and A µ c⊥ . . . A ν c⊥ . Since the effective collinear fields are gauge invariant by themselves, composite operators built out of these fields can be non-local along the light-like directionn without leading to any power suppression [18][19][20][21]. The matrix elements of the bilocal quark-antiquark operators between a meson state and the vacuum can be parameterized in terms of the leading-twist quark LCDA. Specifically, one defines the flavor-specific quark-antiquark LCDAs where [tn, 0] = W c (tn) W † c (0) is a Wilson line extending from 0 to the point tn, and q can be any quark flavor. The quark LCDAs are normalized such that 1 0 dx φ q P (x, µ) = 1. The flavor-specific decay constants f q P entering in (7) are defined in terms of the local matrix elements Note that, due to the axial anomaly, the flavor-diagonal axial currents are not conserved. As a consequence the decay constants f q P (µ) are scale-dependent quantities with an anomalous dimension that starts at two-loop order [22].
For the matrix element of the bilocal two-gluon current we define where T F = 1/2, ⊥ µν (n ·n) = µναβn α n β , and we use the convention that 0123 = +1. The trace in this expression acts in color space. The normalization to the flavor-singlet sum of the light-flavor decay constants, f uds and hence its normalization integral vanishes. Using the integral representation [21] it is straightforward to show that (9) is equivalent to the more conventional definition [23] 1 whereG µν = 1 2 µναβ G αβ is the dual field-strength tensor, and [tn, 0] AB denotes a Wilson line in the adjoint representation of the gauge group. The gauge-invariant bilocal matrix elements in (7) and (9) can be multiplied by functions of the coordinate t. After Fourier transformation to momentum space, these function become the hard-scattering kernels.
The diagrams in Figure 1 produce all n f active quark flavors with an amplitude proportional to Q q v q . This feature remains true when QCD corrections are included. We can decompose the result into a flavor-singlet and a flavor-nonsinglet contribution. In addition, at NLO in α s there exist diagrams of the form shown in Figure 2, in which the meson P is produced via a two-gluon state. In these graphs all possible quark flavors including the top quark contribute in the loop. At a high matching scale µ ∼ m Z the heavy-particle scales m t and m Z are integrated out and absorbed into short-distance coefficient functions. Adding up the various contributions, we obtain the factorization formula 1 Our distribution amplitude φ g P is equal to C F /6 times the function φ (g) M defined in (18) of [24], − C F /12 times the function φ P g defined in (A.10) of [25], and −C F /6 times the LCDA assumed in Sections 3 and 4 of [25] and in [26]. where n f = 5 is the number of active quark flavors in the effective theory below the scale µ ∼ m Z , and in general we define The terms in the first line in (12) correspond to the flavor-singlet contributions to the form factors, involving both quark and gluon LCDAs. The terms in the second line are the combined flavor-nonsinglet contributions, with q c (n f ) q = 0. In Table 1 we collect the relevant coefficients for different numbers of active flavors. We will discuss in Section 2.2 how (12) is evolved down to a low value µ 0 = 1 GeV of the factorization scale. In this process the heavy bottom and charm quarks are integrated out, and the factorization formula must be matched onto an analogous formula in a low-energy effective theories with n f = 3 active flavors.
In order to compute the hard-scattering kernels in (12) at NLO in perturbation theory we evaluate the diagrams shown in Figure 1 (plus six other one-loop graphs) for H (S) q (x, µ), and those shown in Figure 2 for H g (x, µ), using dimensional regularization with d = 4 − 2 spacetime dimensions and working in the MS scheme. In practice, one evaluates these diagrams with on-shell external quark and gluon states carrying momenta xk µ and (1−x)k µ , and then applies projections onto the meson LCDAs. The relevant projections for the quark LCDAs have been discussed in [7]. For the gluon case one computes the partonic amplitude in the form where A, B are color indices. The corresponding hadronic amplitude is then obtained by stripping off the gluon polarization vectors and contracting the indices with the projector [23] The individual loop graphs contain divergences. The ultraviolet divergences cancel in the sum of all diagrams, while infrared divergences are subtracted when we renormalize the meson LCDAs (including a finite renormalization of the axial current). For the flavor-nonsinglet case this has been discussed in detail in [7]. In the flavor-singlet case one needs to account for the mixing between quark and gluon LCDAs. Organizing the products q f q P φ q P (x) and f uds P φ g P (x) into a two-component vector , we can express the bare functions in terms of renormalized functions via Up-type quarks Down-type quarks Table 1: Flavor-number dependent coefficients entering the factorization formula (12).
where the matrix Z f φ (x, y, µ) of renormalization factors is given by Explicit expressions for the kernel functions V ij (x, y) were obtained in [11][12][13][14] and are collected in Appendix A. From these equations the counterterms required to cancel the IR divergences in the bare hard-scattering coefficients are derived. At NLO in α s (but not beyond), the flavor-singlet hard-scattering kernel H S q (x, µ) is given by the same expression as H q (x, µ) in (3). For the gluon kernel we find where the first term accounts for the contributions of the n f = 5 light quark flavors, while the second term describes the contribution of the heavy top quark. We obtain where r t = m 2 Z /m 2 t . The function h g agrees with the corresponding kernel for the γγ * → gg scattering amplitude calculated in [25]. The function h t g can be derived from the expression for the charm-quark contribution to the γγ * → η ( ) form factors presented in [24]. We find that the top-quark contribution is extremely small and can be neglected for all practical purposes.
It will be convenient to transform the factorization formula (12) from momentum space to Gegenbauer moment space, as this turns the convolution integrals into simple sums. To this end one expands the meson LCDAs into Gegenbauer polynomials, such that 2 2 Our Gegenbauer moments b P n are equal to 2 9 times the Gegenbauer moments c (g) n,M used in [24], and − 1 135 times the Gegenbauer moments B g P n and a g P n used in [25] and [26].
C-parity implies that the quark LCDA is an even function under x ↔ (1 − x), while the gluon LCDA is odd. When the Gegenbauer expansions are used, the factorization formula (12) in the 5-flavor effective theory can be recast in the form where the sums run over even integers n = 0, 2, 4, . . . , and we define a P,q is the moment-space expression for the gluon kernel H g (x, µ). Note that we have extracted an overall factor 6 for convenience. These coefficients can be calculated using a technique described in [7]. At one-loop order we obtain Here H n+1 = n+1 k=1 1 k are the harmonic numbers. The coefficients C n (µ) and C S n (µ) start to differ from two-loop order on.

Renormalization-group evolution and resummation
Due to the large energy released to the final-state particles, the Z → P γ decay amplitudes receive contributions from quantum fluctuations with virtualities ranging from the high scale m Z down to the hadronic scale µ 0 ∼ 1 GeV characteristic for light mesons. Phenomenological input for the decay constants f q P (µ) and the meson LCDAs φ q P (x, µ) and φ g P (x, µ) is usually provided at such a low scale. The enormous scale hierarchy gives rise to large logarithms of the form α s ln(m 2 Z /µ 2 0 ) n in the expressions for the hard-scattering kernels -see e.g. (22) -which need to be resummed to all orders of perturbation theory. A subtlety is that the appropriate effective theory at different scales µ between m Z and µ 0 contains different number of active quark flavors. At the high scale µ Z ∼ m Z the appropriate effective theory contains n f = 5 active flavors, which up to corrections of order (m 2 q /m 2 Z ) can be treated as massless. When the factorization scale is lowered below the threshold of the b-quark mass one must match relation (20) onto a corresponding expression in an effective theory containing n f = 4 light flavors. The flavor-singlet and flavor-nonsinglet contributions need to be rearranged in this step in order to account for the change in the coefficients collected in Table 1. At the same time, the contribution from the b-quark stops running below µ b ∼ m b and hence it needs to be evaluated at that scale. A similar procedure takes place when crossing the charm threshold at a scale µ c ∼ m c . The final factorization theorem at the low scale µ 0 ∼ 1 GeV can thus be written in the form Numerical estimates of the hadronic input parameters -the various decay constants and Gegenbauer moments -will be provided in Section 2.3. The scale evolution of the hard-scattering coefficients C (S) n (µ) and D n (µ), and of the decay constants and LCDAs, is controlled by RG evolution equations. The terms shown in the first and second line in the factorization formula (12) are separately scale independent. However, in the first line there is a non-trivial mixing between the quark-and gluon-initiated contributions. Relation (15) implies the RG evolution equation where the anomalous-dimension matrix is given by The quantity Z [1] f φ is the coefficient of the 1/ pole in the matrix Z f φ . The eigenfunctions of the one-loop evolution kernels collected in Appendix A are Gegenbauer polynomials of rank 3/2 and 5/2, as used in (19). At one-loop order the RG equation (24) thus implies an evolution equation in the space of Gegenbauer moments, which is diagonal in n, namely where [11][12][13][14] γ qq n = 2C F 4H n+1 − 3 −

(n + 1)(n + 2)
, γ qg n = −T F n f 40n(n + 3) 3(n + 1)(n + 2) , γ gq n = −C F 12 5(n + 1)(n + 2) , The corresponding evolution equation for flavor-nonsinglet combinations of the Gegenbauer moments of quark LCDAs is multiplicative and governed by the anomalous-dimension coefficients γ qq n . The solution of equation (26) is discussed in Appendix A. One finds that both eigenvalues of the moment-space anomalous-dimension matrix are positive and grow logarithmically at large n, and the same holds for the coefficients γ qq n . It follows that all Gegenbauer moments vanish for asymptotically large values of the renormalization scale, a P,q n (µ) → 0 and b P n (µ) → 0 for µ → ∞, irrespective of their values at a low hadronic scale. In momentum space, this implies the asymptotic forms The high value of the hard matching scale µ ∼ m Z in Z → P γ decays ensures that one is rather close to the asymptotic limit [7], and this fact reduces the sensitivity of our predictions to the Gegenbauer moments of the LCDAs, which are currently not known with good accuracy. RG invariance of the form factors in (20) implies that the flavor-singlet hard-scattering coefficients obey the evolution equation The flavor-nonsinglet coefficients C n (µ) obeys an analogous equation without mixing, in which the anomalous-dimension matrix is replaced by the coefficients γ qq n . The solution to these evolution equations obtained at leading order can be written in the form Explicit expressions for the evolution functions U S n (µ 1 , µ 2 ) and U n (µ 1 , µ 2 ) are collected in Appendix A. We use (30) to perform the evolution of the Wilson coefficients in the intervals between flavor thresholds. Concretely, we calculate the initial conditions for the coefficient functions at a high scale µ Z ∼ m Z from (21). Notice that these relations are free of large logarithms at this scale. We then evolve the coefficients to a matching scale µ b ∼ m b using (30) and evaluating the evolution functions with n f = 5 active quark flavors. At the scale µ b the b quark is integrated out from the effective theory, and in this process we need to perform a careful matching of the coefficients in the 5-flavor and 4-flavor theories. The relevant matching conditions are We then continue to run the coefficients down to a matching scale µ c using (30) with evolution functions corresponding to n f = 4 active quark flavors. At the scale µ c the charm quark is integrated out. The relevant matching conditions are Finally, we run the coefficients down to the low-energy scale µ 0 using (30) with evolution functions corresponding to n f = 3 active quark flavors. The values we obtain for the Wilson coefficients of the first few Gegenbauer moments at the scale µ 0 = 1 GeV are collected in Table 2. Two points are worth mentioning. First, for n ≥ 2 the coefficients D n of the gluon LCDA are of the same magnitude as those of the quark LCDAs, which is remarkable given that the former ones start at O(α s ), whereas the latter ones start at tree level. Second, we observe that the differences between the flavor-singlet coefficients C S n and the flavor-nonsinglet coefficients C n are numerically very small. This shows that the mixing of quark and gluon LCDAs under RG evolution is a small effect.
Beyond one-loop order the solution of the RG equation (26) takes a much more complicated form, since Gegenbauer moments of different rank mix under renormalization. The corresponding expressions for the singlet and nonsinglet cases can be derived from [27,28] and are collected in the appendix of [24]. A dedicated study of NLO evolution effects for the case of exclusive radiative Higgs-boson decays has been performed in [10]. It turns out these effects have only a modest impact on our results. For example, the values of the nonsinglet coefficients C n (µ 0 ) shown in Table 2 are reduced by factors of 1, 0.93, 0.92, 0.91 for n = 0, 2, 4, 6 when NLO evolution effects are taken into account. The most important impact on the singlet coefficients is due to the anomalous dimension of the flavor-singlet axial current, which starts at two-loop order [22]. As a result, when NLO evolution effects are included the value of the leading singlet coefficient C S 0 (µ 0 ) differs from the value given in Table 2 by the factor Given the present large uncertainties in the values of the Gegenbauer moments with n = 0, it is a reasonable approximation to account for the effects of NLO evolution by treating κ NLO as a global prefactor in the expressions for the form factors and otherwise use the coefficients compiled in Table 2. The deviation of κ NLO from 1 provides a measure of the remaining perturbative uncertainty in our predictions.

Hadronic input parameters
It remains to obtain the required hadronic input for the decay constants and Gegenbauer moments of the η and η mesons. Following [23], we assume isospin symmetry of all hadronic matrix elements, but we differentiate between the matrix elements of operators containing up and down quarks and those containing strange quarks. In the SU (3) flavor-symmetry limit, the pseudoscalar meson η would be a flavor octet and η a flavor singlet. However, it is known empirically that SU (3)-breaking corrections to these assignments are large. In the following we shall thus not rely on SU (3) flavor symmetry but instead introduce another assumption, expected to be accurate at the 10% level. In the absence of the axial U (1) anomaly, the flavor states |η q = (|uū + |dd )/ √ 2 and |η s = |ss mix only through OZI-violating effects, which are known phenomenologically to be small. It is therefore reasonable to assume that the axial anomaly is the only effect that mixes these two flavor states. This is the basis of the FKS mixing scheme [15]. Since this is by assumption the only mixing effect, the FKS scheme amounts to a scheme with a single mixing angle, which relates the physical mass eigenstates to the flavor states via In the FKS mixing scheme, one introduces decay constants f q and f s and light-cone distribution amplitudes φ q (x, µ) and φ s (x, µ) for the flavor states |η q,s [23]. The physical η and η mesons are then described as coherent superpositions of these states. It follows that the flavor-specific decay constants defined in (8) are given by Analogous relations hold for the flavor-specific LCDAs defined in (7), i.e.
Moreover, in the FKS scheme one sets φ g η (x, µ) = φ g η (x, µ) ≡ φ g (x, µ) for the gluon LCDAs [23,24,26]. As a consequence, the hadronic parameters characterizing the quark and gluon LCDAs are the Gegenbauer moments a q n and a s n of φ q and φ s and the Gegenbauer moments b n of φ g , defined as in (19).
The decay constants f q , f s and the mixing angle ϕ in the FKS scheme have been determined in [15] from a fit to experimental data, finding Here f π = (130.4±0.2) MeV is the pion decay constant [17]. A more recent analysis exploiting additional data but only a subset of the processes investigated in the original paper finds [29] f q = (1.09 ± 0.03) f π , f s = (1.66 ± 0.06) f π , ϕ = 40.7 • ± 1.4 • .
The central values of the flavor-singlet decay constants f uds η ( ) , which provide the normalization for the leading contributions to the form factors in (23) to the choice of FKS parameters will be reflected in the spread of our predictions for the Z → ηγ branching ratio.
The Gegenbauer moments of the LCDAs can be extracted using fits to data for the γ * γ → η ( ) transition form factors at different Q 2 reported by the CLEO [30] and BaBar [31] collaborations. The authors of [24] have assumed SU (3) flavor symmetry and have chosen the first few Gegenbauer moments of the quark LCDAs φ q and φ s according to some popular QCD sum-rule calculations for the pion LCDA. The first Gegenbauer moment of the gluon LCDA φ g has then been extracted from the fit to the data. In this context three benchmark models were identified, to which we refer below as models (i)-(iii). In all three cases the FKS parameters in (37) have been used. The authors of [26], on the contrary, have extracted the first Gegenbauer moments of both the quark and the gluon LCDAs from fits to the data. Below we will consider three of their fit scenarios. Model (iv) refers to their default fit, while model (v) corresponds to a fit exclusively to BaBar data. In both cases the FKS parameters  Table 3: Gegenbauer moments of quark and gluon LCDAs at the scale µ 0 = 1 GeV in different benchmark models obtained from analyses of the γ * γ → η ( ) transition form factors. Models (i)-(iii) correspond to the models in Table 2 of [24], while models (iv)-(vi) refer to the first, third and sixth model in Table 2 of [26].
in (37) are assumed. Model (vi) corresponds to a combined fit to CLEO and BaBar data using the FKS parameters in (38). Table 3 collects the values of the Gegenbauer moments in the six benchmark models, translated to our notations.
To evaluate the form factors in (23) we also need the decay constants f c P and f b P describing the intrinsic charm and bottom contents of the η and η mesons. Following [23], we estimate these parameters using relations among the FKS parameters implied by the axial anomaly. This yields Numerically, we obtain f c η ≈ −1.
2 MeV, f c η ≈ −2.9 MeV and f b η ≈ −0.1 MeV, f b η ≈ −0.3 MeV. These values are of the same order as those obtained using similar methods in [15,[32][33][34][35]. Due to the very small effect of the intrinsic charm and bottom contributions on the form factors we only keep the leading terms with n = 0 in these contributions and set all Gegenbauer moments a P,c n and a P,b n with n > 0 to zero.

Phenomenological predictions
We are now ready to present our numerical results. Before we quote our predictions for the Z → η ( ) γ branching ratios, we demonstrate the sensitivity of the form factors in (23) to the choice of the FKS parameters, the decay constants f c,b P , the Gegenbauer moments of the LCDAs and the factorization scale. Using the mixing parameters in (37) as our default values, we obtain for the real parts of the form factors The small imaginary parts of O(α s ) can be neglected, since they do not contribute to the decay rates at NLO in RG-improved perturbation theory. The parameters f q , f s and ϕ used to compute these values come with uncertainties, which we did not take into account in the expressions above for readability. We demonstrate their effect by looking only at the leading terms, for which we find Re F η = κ NLO 16.3 ± 1.5 ϕ ± 1.0 fq ± 1.6 fs MeV + . . . , We now recompute these results using the FKS parameters collected in (38). In this case we find and the parametric uncertainties in the leading terms read The contributions from the gluon LCDA to the form factors are important. For the η meson the coefficients of the gluon Gegenbauer moments b 2 and b 4 are significantly larger than those of the quark Gegenbauer moments of the same order. In the case of the η meson the individual quark Gegenbauer moments come with large coefficients, but there are large cancellations between the terms involving a q n and a s n with the same n. If one assumes SU (3) flavor symmetry as in [24], then the contributions of the quark and gluon Gegenbauer moments at the same order come again with similar coefficients, but obviously SU (3)-breaking corrections could have a large impact. If future theoretical efforts based on lattice gauge theory or refined QCD sum rules can provide a better control of the Gegenbauer moments of the quark LCDAs, then both decay channels could potentially provide valuable information about the moments of the gluon LCDA.  Table 4: Central values of the Z → η ( ) γ branching ratios in units of 10 −9 , obtained using six different models of hadronic input parameters, see Table 3. Models (i)-(v) use the mixing parameters in (37) Given our results for the form factors F P , the Z → P γ branching ratios are obtained as where α = 1/137.036 is the electromagnetic fine structure constant evaluated at q 2 = 0, Γ Z = 2.4955 GeV is the Z-boson width, and v = 245.36 GeV is the Higgs vacuum expectation value [17]. Table 4 shows our predictions for the Z → η ( ) γ branching fractions obtained using the hadronic parameters corresponding to the six models collected in Table 3 along with κ NLO = 0.89. The sensitivity of the branching ratios to the hadronic input parameters opens up the possibility of probing these parameters using the decays Z → η ( ) γ. This is particularly interesting since the different LCDA parameters from Table 3 are all obtained from fits to the same low-energy experimental data, but still differ a lot from each other depending on the assumptions and methods used to perform these fits. Owing to the very large value µ Z ∼ m Z of the factorization scale inherent to Z → P γ decays, the analysis of these processes is much cleaner theoretically and less affected by uncertainties due to power corrections or higher-order perturbative corrections. The branching ratios we obtain are of similar magnitude as those for the corresponding decays to light vector mesons found in [7], e.g. Br(Z → ργ) ≈ 4.2 · 10 −9 and Br(Z → φγ) ≈ 8.6 · 10 −9 . A future high-luminosity e + e − collider operating at the Z pole could produce samples of about 10 12 Z bosons per year [36]. This would yield at best O(150) events in the Z → ηγ channel and O(4500) events in the Z → η γ channel without considering backgrounds and reconstruction efficiencies.

Conclusions
We have presented a detailed analysis of the rare radiative decays Z → ηγ and Z → η γ within the QCD factorization approach, working at NLO in RG-improved perturbation theory. In particular, we have included the additional contributions that arise since the final-state mesons η and η contain a flavor-singlet component in their wave function, which can be formed from a two-gluon state. We have derived the corresponding QCD factorization formula in the framework of SCET, generalizing the derivation presented in [7] to the flavor-singlet case. In the presence of the gluon contribution the scaling behavior of the hadronic matrix elements is changed, because quark and gluon LCDAs mix under renormalization. This complicates the resummation of large logarithms. Since the flavor singlet is composed of a different number of dynamical quark flavors at the high scale µ ∼ m Z compared to the low hadronic scale µ ∼ 1 GeV, the singlet and nonsinglet matrix elements need to be rearranged when crossing quark flavor thresholds, giving rise to non-trivial matching conditions. Also, the intrinsic charm and bottom content of the η ( ) mesons gives rise to non-trivial effects.
Many of the hadronic input parameters needed for a phenomenological analysis of our calculations are not yet well determined. In the literature, these parameters have been obtained from different approaches leading to sometimes quite different results. Considering a set of six benchmark models obtained in [24,26], the values we obtain for the Z → ηγ branching ratio range from 0.1 · 10 −10 to 1.7 · 10 −10 , while those for the Z → η γ branching fraction vary between 3.1 · 10 −9 and 4.8 · 10 −9 . We find that especially the Z → ηγ channel is very sensitive to the η −η mixing parameters, so that a measurement of this decay mode could be used to favor one particular set of mixing parameters over another. A similar statement can be made with regard to the different sets of LCDA shape parameters. While all models in Table 3 seem to give consistent descriptions of the γ * γ → η ( ) transition form factors, we find that some of them give rather distinct predictions for the Z → η ( ) γ branching ratios. Measurements of these rare radiative Z-boson decays could thus be used as complementary information to attain a better control over the hadronic matrix elements describing the η−η system. If in the future the quark LCDAs and meson mixing parameters can be determined more accurately, using more precise low-energy data and advanced theoretical tools such as lattice gauge theory, a study of the Z → η ( ) γ decay modes at a future lepton collider would provide a very good opportunity to directly access the Gegenbauer moments of the gluon LCDA in a theoretically clean environment.