Exclusive Radiative Decays of W and Z Bosons in QCD Factorization

We present a detailed theoretical analysis of very rare, exclusive hadronic decays of the electroweak gauge bosons V=W, Z from first principles of QCD. Our main focus is on the radiative decays V->M+gamma, in which M is a pseudoscalar or vector meson. At leading order in an expansion in powers of Lambda_{QCD}/m_V the decay amplitudes can be factorized into convolutions of calculable hard-scattering coefficients with the leading-twist light-cone distribution amplitude of the meson M. Power corrections to the decay rates arise first at order (Lambda_{QCD}/m_V)^2. They can be estimated in terms of higher-twist distribution amplitudes and are predicted to be tiny. We include one-loop O(alpha_s) radiative corrections to the hard-scattering coefficients and perform the resummation of large logarithms [alpha_s log(m_V^2/mu_0^2)]^n (with mu_0=1 GeV a typical hadronic scale) to all orders in perturbation theory. Evolution effects have an important impact both numerically and conceptually, since they reduce the sensitivity to poorly determined hadronic parameters. We present detailed numerical predictions and error estimates, which can serve as benchmarks for future precision measurements. We also present an exploratory study of the weak radiative decays Z->M+W. Some of the decay modes studied here have branching ratios large enough to be accessible in the high-luminosity run of the LHC. Many of them can be measured with high accuracy at a future lepton collider. This will provide stringent tests of the QCD factorization formalism and enable novel searches for new physics.

1 Introduction a single meson as a laboratory to test and study the QCD factorization approach in a context where power corrections are definitely under control. The enormous rates of W and Z bosons that will become available at future colliders will present us with a new playground for precision electroweak and QCD physics, which will make such studies feasible. With 3000 fb −1 collected during the high-luminosity run at the LHC, one will have produced more than 10 11 Z bosons and 5 · 10 11 W bosons in both ATLAS and CMS. A clean, tagged sample of W bosons is expected to come from top-quark decays [15]. At a future lepton collider such as TLEP, samples of up to 10 12 Z bosons per year can be expected in a dedicated run at the Z pole [16]. Our main focus is on the simplest processes: the hadronic radiative decays Z → Mγ and W → Mγ, where M is a pseudoscalar or vector meson. They offer a perfect way to probe some properties of the leading-order LCDAs of various mesons. The price one needs to pay is that the higher the energy release in the process, the smaller the probabilities for any particular exclusive final state are. The branching fractions we obtain range from few times 10 −8 to few times 10 −11 or even smaller. The big challenge of such a program will be to measure such decays experimentally with some precision. While we do not perform a detailed feasibility study in this work, we speculate that some of these rare modes will be accessible in the highluminosity run of the LHC, at a level that will be useful to probe our theoretical predictions. At future lepton colliders operating on the Z pole, it would be possible to measure several of these decays at or below the 1% level. This may present us with a unique opportunity to extract information about LCDAs in a theoretically clean environment.
Our interest in this subject was raised by recent investigations of the exclusive decays h → V γ [17][18][19][20] and h → ZV [21,22] of the Higgs boson to final states containing a single vector meson. It was proposed to use these decays as a way to probe for possible non-standard Yukawa couplings of the Higgs boson. Such measurements are extremely challenging at the LHC and other future colliders. Observing exclusive hadronic decays of W and Z bosons would provide a proof of principle that this kind of searches can be performed. An encouraging first search for the decays Z → J/ψ γ and Z → Υ(nS) γ has just been reported by ATLAS [23].
From a theoretical perspective, the very rare, exclusive radiative decays of W and Z bosons have received relatively little attention in the literature, and very few accurate predictions for such branching fractions have been obtained. In a pioneering study [24], Arnelos, Marciano and Parsa presented a first detailed analysis of the decays W → P γ and Z → P γ for both light and heavy pseudoscalar mesons in the final state. Strong-interaction effects were parametrized in terms of vector and axial-vector form factors, which were estimated using ideas from perturbative QCD on the asymptotic behavior of form factors at large momentum transfer. Several years later, Manohar studied the decays Z → πW and Z → πγ using a local operator-product expansion [25], which expresses the decay amplitudes as power series in parameters 26 and ω 0 = 2, respectively. If only the leading term is kept, the amplitudes can be related to the pion matrix element of the axial-vector current, which is proportional to the pion decay constant f π . For radiative decays such as Z → πγ this truncation cannot be justified theoretically, and the infinite tower of local operators would need to be resummed. In a very recent work, the method developed by Manohar was used to derive an estimate for the W → πγ branching fraction [15]. While the analyses performed in these papers can provide some order-of-magnitude results, they do not allow to obtain accu-rate predictions with reliable error estimates. In a classic paper [26], Guberina et al. analyzed the radiative decays of the Z boson into heavy quarkonia in the non-relativistic limit. The first relativistic corrections to the Z → J/ψ γ and Z → Υ(1S) γ decay rates were added only recently in [27], where in addition the authors considered for the first time the decay Z → φγ, using an approach closely related to ours. As we will discuss, renormalization effects have a profound impact on the decay amplitudes. When evolved up to the relevant scales of order the Z-boson mass, the LCDAs of heavy quarkonia can no longer be accurately described by the leading term in a non-relativistic expansion.
In the present work, we present a comprehensive analysis of a large class of radiative decays of W and Z bosons using the QCD factorization approach, including for the first time a consistent treatment of O(α s ) corrections and performing the resummation of large logarithms of order α s ln(m 2 Z /µ 2 0 ) n , with µ 0 ≈ 1 GeV, to all orders in perturbation theory. 1 Our approach provides a systematic expansion of the decay amplitudes in powers of the small parameters α s (m Z ) ∼ 0.1 and Λ QCD /m Z ∼ 0.01. We study the structure of the leading power corrections to the Z → Mγ and W → Mγ decay rates and show that they are of second order and hence negligibly small, of order 10 −4 relative to the leading terms. For processes involving heavy quarks, power corrections of order (m Q /m Z ) 2 exist, which are still very small (less than 1%) even for final-state mesons containing b-quarks. Finally, using the most recent experimental data we perform a reanalysis of meson decay constants, which provide crucial input to our phenomenological analysis. Our paper is organized as follows: In Section 2 we derive a factorization theorem for the Z → Mγ and W → Mγ decay amplitudes, in which they are expressed as convolutions of calculable hard-scattering kernels with meson LCDAs. We explain how the kernel functions can be calculated by performing projections of on-shell partonic amplitudes. We then summarize the existing theoretical information on the shapes of the LCDAs for both light and heavy mesons and study their behavior under scale evolution. In Section 3 we apply this approach to derive explicit predictions for the Z → Mγ and W → Mγ decay amplitudes at leading power in Λ QCD /m Z . The relevant convolution integrals of hard-scattering kernels with LCDAs are calculated in analytic form using an expansion in Gegenbauer polynomials. We demonstrate that renormalization-group (RG) evolution from a low hadronic scale up to the electroweak scale of relevance to these processes has the nice effect of significantly reducing the sensitivity to poorly determined hadronic parameters. By studying radiative decays into transversely polarized vector mesons, we present some detailed estimates of power-suppressed effects in the QCD factorization approach. In some old papers, it was suggested that the radiative decay amplitudes into pseudoscalar mesons can be hugely enhanced due to effects of the axial anomaly [28,29]. We explain why such an enhancement does not exist. We then present our numerical predictions for W, Z → Mγ branching fractions, including detailed error estimates. The results span more than three orders of magnitude, and we explain the striking differences seen between the various decay channels in terms of electroweak couplings, differences in decay constants, and enhancement factors occurring for heavy-light mesons. In Section 4 we present a first exploratory study of the weak radiative decays Z → MW , in which a heavy W boson is part of the final state. In this case significantly less energy is released to the final-state meson M, and as a result the QCD factorization approach can be tested at energies of order 10 GeV, about a factor 2 higher than those available in exclusive B-meson decays. We round off our study in Section 5 with some experimental considerations. Our main results are summarized in Section 6. Technical details of our calculations and the extraction of meson decay constants are relegated to three appendices.

Theoretical framework
Our main focus in this work is on the rare, exclusive radiative decays Z → Mγ and W → Mγ, where M denotes a pseudoscalar or vector meson. We assign momentum k to the final-state meson and q to the photon. The leading-order Feynman diagrams for the case of Z → Mγ are shown in Figure 1. The decay plane is spanned by the vectors k and q. We will refer to vectors in this plane as being longitudinal, and to vector orthogonal to it as being transverse. We only consider cases where the mass of the final-state meson satisfies m M ≪ m Z . Up to corrections suppressed as (m M /m Z ) 2 , this mass can then be set to zero. In this limit, we have k µ = En µ and q µ = En µ , where E = m Z /2 is the energy of the final-state particles in the Z-boson rest frame, and n andn are two light-like vectors satisfying n ·n = 2.

Derivation of the factorization formula
For the purposes of this discussion we work in the rest frame of the decaying heavy boson. The decay amplitudes can be calculated from first principles using the QCD factorization approach [1][2][3][4][5], 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 , they can be written as convolutions of calculable hard-scattering coefficients with LCDAs of the meson M. A simple way to derive the corresponding factorization theorem employs the formalism of SCET [10][11][12][13]. It provides a systematic expansion of decay amplitudes in powers of a small expansion parameter λ = Λ QCD /E. The light final-state meson moving along the direction n µ can be described in terms of collinear quark, anti-quark 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 [30,31] where iD µ c = i∂ µ + igA µ c denotes the covariant collinear derivative, and is a collinear Wilson line extending from x to infinity along the directionn. Both fields are of O(λ) in SCET power counting. Adding more component fields to an operator always leads to further power suppression. At leading order in λ, the operators with a non-zero matrix element between the vacuum and a single meson state are thus of the formX c (tn) . . . X c (0) and A µ c⊥ (tn) . . . A c⊥µ (0), where without loss of generality we set x = 0 for one of the fields. 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. The two-gluon operator would only be relevant for decays into mesons containing a flavor-singlet component on their wave functions, such as the pseudoscalar mesons η and η ′ [32]. Such decays will be discussed in a forthcoming publication [33]. It follows that at leading power in the expansion in λ, the Z → Mγ and W → Mγ decay amplitudes into non-singlet final states can be written in the factorized form where µ is the factorization scale, and Γ i ∈ {1, γ 5 , γ µ ⊥ }. The four matrices (/ n/2) Γ i provide a basis of Dirac matrices sandwiched between two collinear quark spinors. The Wilson coefficients C i (t) are process dependent and can be calculated perturbatively. In the last step we have used the definition (1) and combined the two Wilson lines W c (tn) W † c (0) ≡ [tn, 0] into a straight Wilson line extending from 0 to tn. The meson matrix elements of the bi-local operators in the second line define the leading-order LCDAs of pseudoscalar and vector mesons. Specifically, one has where E =n · k/2 denotes the energy of the meson in the rest frame of the decaying boson, f P and f V are the decay constants of pseudoscalar and vector mesons defined in terms of their matrix elements of local (axial-)vector currents, and f ⊥ V (µ) is a scale-dependent vector-meson decay constant defined in terms of a matrix element of the QCD tensor current. The leadingorder LCDAs can be interpreted as the amplitudes for finding a quark with longitudinal momentum fraction x insinde the meson. The factor of γ 5 in the first equation is present for a pseudoscalar meson (M = P ) but absent for a longitudinally polarized vector meson (M = V ). The projection onto a transversely polarized vector meson does not arise at leading power in the radiative decays of W and Z bosons. For a given meson, exactly one of the possible Dirac structures contributes, and we denote the corresponding Wilson coefficient by C M (t, µ).
Defining the Fourier-transformed Wilson coefficient, called the hard function, via we obtain the factorization formula Insertions of additional collinear fields or derivatives yield power-suppressed contributions. In particular, the insertion of an additional collinear gluon field gives rise to three-particle LCDAs.
In order to fully establish the factorization theorem (6) one must show that the convolution integral over the momentum fraction x converges at the endpoints. This question has been addressed in the context of the more complicated processes B → γlν [34] and B → K * γ in [35]. The behavior near the endpoints can be described by means of soft-collinear fields [36,37] with momenta scaling as (n · p sc ,n · p sc , p ⊥ sc ) ∼ E(λ 2 , λ, λ 3/2 ). The contributions of such modes are always power suppressed. In the present case, we find that endpoint singularities are absent at leading and subleading power in the large-energy expansion.
LCDAs play the same role for hard exclusive processes which PDFs play for inclusive ones. While they encode genuinely non-perturbative hadronic physics, they can be rigorously defined in terms of non-local operator matrix elements in QCD [1][2][3][4][5]. These matrix elements can be systematically expanded in terms of structures of different twist. When applied to highenergetic exclusive processes such as the ones considered here, the twist expansion translates into an expansion in powers of Λ QCD /E. There is an extensive amount of literature devoted to the study of distribution amplitudes. For light pseudoscalar mesons, the two-and threeparticle LCDAs up to twist-3 order were studied, e.g., in [38], while the corresponding LCDAs for vector mesons were analyzed, e.g., in [39][40][41]. We stress that, at the scale of the large energies released in decays of W and Z bosons, even charm and bottom quarks can be treated as light quarks, and hence heavy mesons containing these quarks can be described by LCDAs. This will be discussed further below.
In order to apply these results in practical calculations, it is convenient to define momentumspace projection operators, which can be applied directly to the decay amplitudes computed with on-shell external parton states [8,42]. For all two-particle projections onto LCDAs of leading and subleading twist, it is sufficient to assign momenta k 1 = xk + k ⊥ + . . . and k 2 = (1 − x)k − k ⊥ + . . . to the quark and the anti-quark in the meson M, where k is treated as a light-like vector (k 2 = 0). Meson mass effects of order m 2 M enter only at twist-4 level. They have a tiny numerical impact for the decays considered here, and we will consistently set m 2 M → 0 unless noted otherwise. The variables x and (1 − x) denote the longitudinal momentum fractions carried by the quark and the anti-quark in the two-body Fock state of the meson. Each Feynman diagram gives an expression of the form where in the last step we have introduced the light-cone projection operator M M (k, µ) for the meson M, which at higher order contains derivatives with respect to the parton transverse momentum k ⊥ . It is understood that k ⊥ is set to zero after these derivatives have been performed. Up to twist-3 order, the light-cone projector for a pseudoscalar meson can be written in the form [8,42] Here φ P is the leading-twist LCDA of the meson, while φ p and φ σ denote the two-particle LCDAs appearing at twist-3 order. These are scale-dependent functions, which we define in the MS renormalization scheme. The decay constant f P of the meson P is defined in terms of its matrix element of a local axial-vector current The scale-dependent parameter µ P (µ) = m 2 P /[m q 1 (µ) + m q 2 (µ)] governs the normalization of the twist-3 LCDAs. 2 The vectorn in the above expression denotes a longitudinal light-like vector not aligned with k. A convenient choice is to take the photon momentum,n = q. At twist-3 order the projector also contains three-particle LCDAs containing a quark, an anti-quark and a gluon. We will see that the contributions of twist-3 LCDAs are strongly suppressed compared with those of the leading-twist amplitudes. In order to estimate their effects, we will for simplicity neglect the three-particle LCDAs. This is referred to as the Wandzura-Wilczek approximation (WWA) [43]. When this is done, the QCD equations of motion fix the form of the twist-3 LCDAs completely, and one obtains [38] The light-cone projection operators for vector mesons are more complicated. They are given in Appendix A. For our purposes it suffices to quote the projector for a longitudinally polarized vector meson at leading power. It is 2 Note that µ π = m 2 π /(m u + m d ) holds for charged and neutral pions, see e.g. [32].
The function φ V (x, µ) is sometimes called φ V (x, µ) in the literature. We have used that the longitudinal polarization vector is given by ε Vn µ k·n . The vector-meson decay constant f V is defined in terms of the local matrix element Before proceeding, let us comment on the structure of power corrections to the factorization formula (6). Inspecting the explicit form of the projection operator for a pseudoscalar meson in (8), and the corresponding projectors for vector mesons given in (A.1) and (A.4), we observe that consecutive terms in the twist expansion contain even and odd numbers of Dirac matrices in alternating order. Since the gauge interactions in the Standard Model preserve chirality, it follows that for a given helicity amplitude either all terms with an even number of Dirac matrices contribute or all terms containing an odd number, but not both. Consequently, the SCET expansion for the Z → Mγ decay amplitudes with fixed polarizations of all particles is an expansion in powers of (Λ QCD /m Z ) 2 . The power counting changes when quark-mass effects are taken into account. They give rise to chirality-changing vertices, which give corrections suppressed by m Q /m Z to both the amplitudes and the meson projectors. This leads to power corrections of order m Q Λ QCD /m 2 Z and (m Q /m Z ) 2 . For heavy quarks with m Q ≫ Λ QCD , the latter corrections are the dominant ones. However, as long as the relevant quark masses m Q are much smaller than the hard scale m Z of the process, these corrections are still small. The present case is different from the situation encountered in exclusive B-meson decays [6][7][8][9], where the presence of a heavy quark mass, which is of the same order as the energy released in the decay, allows for O(1) chirality-changing interactions. In this case the decay amplitudes receive first-order Λ QCD /m b corrections.

Systematics of the Gegenbauer expansion
The leading-twist LCDAs obey an expansion in Gegenbauer polynomials of the form [1,5] which can be inverted to give a M n (µ) = 2(2n + 3) 3(n + 1)(n + 2) The Gegenbauer moments have a diagonal scale evolution at leading order in perturbation theory. They are non-perturbative hadronic parameters, which can only be accessed using data or a non-perturbative approach such as light-cone QCD sum rules (see e.g. [39][40][41]) or lattice QCD [44]. In Table 1 we collect the values for the decay constants and the first two Gegenbauer moments a M 1,2 for light pseudoscalar and vector mesons. Our notation is such that K ( * ) ∼ (qs) with q = u, d, and x is the momentum fraction of the light quark q.
pion were studied in [52,53] using a QCD sum-rule approach employing non-local vacuum condensates. These authors find a π 2 = 0.20, a π 4 = −0.14, a π 6 = 5 · 10 −3 , and a π 8 = a π 10 = 4 · 10 −3 at the scale µ 0 = 1 GeV. Their value of a π 2 is consistent with the result given in Table 1, while higher moments a π n with n ≥ 6 are estimated to be negligibly small. On the other hand, in more recent work [54] the authors have performed fits to the first eight Gegenbauer moments of the pion LCDA using data on the π 0 γ * γ form factor obtained by the BaBar and Belle collaborations [55,56]. They find a π 2 = 0.10 (0.14), a π 4 = 0.10 (0.23), a π 6 = 0.10 (0.18) and a π 8 = 0.034 (0.050) at µ 0 = 1 GeV for Belle (BaBar), which suggests that a π 6 and a π 8 may not be insignificant. In our phenomenological analysis we will vary a M 4 (µ 0 ) between −0.15 and +0.15 for all light mesons and use this to estimate the effect of unknown higher Gegenbauer moments. With this treatment, the relevant combination of Gegenbauer coefficients given in relation (44) below agrees with all of the above models within our quoted uncertainties.
It is an important question to ask what can be said on general grounds about the behavior of the Gegenbauer expansion. It is commonly assumed, and is supported by power-counting analyses in SCET, that the leading-twist LCDAs vanish at the endpoints x = 0 and x = 1, such that the integrals converge. This statement implies that the infinite sums n a M n and n (−1) n a M n converge. Barring accidental cancellations, this requires that for large n the coefficients a M n fall off faster than 1/n, and this condition should hold for all values of µ 0 . From a physical point of view, high-rank Gegenbauer polynomials C (3/2) n (2x − 1) with n ≫ 1 resolve structures on scales ∆x ∼ 1/n. For a light meson M, it is reasonable to assume that the LCDA φ M (x) does not exhibit pronounced structures at scales much smaller than O(1), in which case the coefficients a M n must decrease rapidly at large n. The LCDAs of heavy mesons are an exception to this rule, since the presence of the heavyquark mass introduces a distinct scale. For a quarkonium state M ∼ (QQ) composed of two identical heavy quarks, the LCDA peaks at x = 1/2 and has a width that tends to zero in the limit of infinite heavy-quark mass. The second moment of the LCDA around x = 1/2 can be related to a local matrix element in non-relativistic quantum chromo-dynamics (NRQCD), the effective field theory describing heavy quarkonia states [57,58]. This framework provides a systematic expansion of hadronic matrix elements in powers of the small velocity v ∼ α s (m Q v) of the heavy quark in the quarkonium rest frame. One obtains [59] To derive this result one uses that in the heavy-quark limit x = p Q ·n 2m Q V ·n = 1+vz 2 , wheren µ is a light-like vector, and p µ Q = m Q V µ + k µ denotes the momentum of the heavy quark inside the quarkonium state with velocity V µ . The various vectors are defined such that V ·n = 1 and V · k = 0. In the rest frame of the quarkonium state we can choose V µ = (1, 0), n µ = (1, −e z ), and k µ = (0, m Q v), where the 3-vector v is the residual velocity of the heavy quark inside the (QQ) bound state. The factor 1/3 on the right-hand side of (15) is due to rotational invariance in the rest frame. Numerical values for the NRQCD matrix element v 2 for the J/ψ and Υ(1S) states have been obtained from an analysis of the leptonic decay rates Γ(J/ψ → e + e − ) and Γ(Υ(1S) → e + e − ) including first-order α s corrections and nonperturbative contributions proportional to v 2 . In this way, the values v 2 J/ψ = 0.225 +0.106 −0.088 [60] and v 2 Υ(1S) = −0.009 ± 0.003 [61] have been extracted, the latter one being inconsistent with the fact that the second moment in (15) must be positive. Both estimates suffer from the fact that the two-loop [62,63] and three-loop [64] perturbative corrections to the NRQCD predictions for these decay rates are known to be huge, precluding a reliable extraction of nonperturbative parameters. Based on the power-counting rules of NRQCD one would naively expect that v 2 J/ψ ∼ 0.3 and v 2 Υ(1S) ∼ 0.1, and we will use these estimates, along with a 50% relative error assigned to them, in our phenomenological analysis. For our calculations we need the first inverse moments of the LCDA with respect to x or (1 − x). Expanding the inverse moments about x = 1/2, it is immediate to derive the model-independent relation [20] As a reasonable model at the low scale µ 0 = 1 GeV we adopt the Gaussian ansatz where the polynomial in front of the Gaussian factor ensures that the LCDA vanishes at the endpoints x = 0, 1. The normalization constant N σ ≈ 1 can be expressed in closed form in terms of an error function. For a heavy-light meson state M ∼ (qQ) composed of a light quark and a heavy anti-quark, the LCDA peaks at a small value x ∼ Λ QCD /m M , where x refers to the momentum fraction of the light spectator quark. The appropriate effective field theory for heavy-light bound states is called heavy-quark effective theory (HQET), see [65] for a review. In the context of this theory, it is possible to show that the first moment of the LCDA is determined by the ratioΛ M /m M , 326 ± 17 -0.10 ± 0.05 0.091 ± 0.023 Table 2: Hadronic input parameters for pseudoscalar and vector mesons containing heavy quarks. Scale-dependent quantities are defined at µ 0 = 1 GeV. The values for f D and f Ds are taken from [45]. The values for f B and f B are taken from two recent, unquenched lattice calculations [70,71], which obtain identical central values but quote very different error estimates. We quote the averages of the uncertainties given by the two groups. The values of the J/ψ and Υ(nS) decay constants can be derived from data, as explained in Appendix B.
where m M denotes the heavy-meson mass andΛ M = m M − m Q (with m Q being the pole mass of the heavy quark) is a hadronic parameter. One obtains , where the one-loop radiative corrections have been calculated in [67] and are numerically significant. In our analysis below we need the first inverse moment of the LCDA with respect to x, which is of order m M /Λ QCD and cannot be related to a local HQET matrix element. One defines [6] where the hadronic parameter λ M (µ 0 ) ∼ Λ QCD is independent of the heavy-quark mass, and the dots denote corrections that are power-suppressed relative to the leading term. The parameter λ M is poorly known at present. A QCD sum-rule estimate for the B meson yields λ B (1 GeV) = (460±110) MeV [68], and we will use this value in our phenomenological analysis for both B and D mesons. Concerning B s and D s mesons, we shall use the estimate λ Ms −λ M ≈ 90 MeV from [69] and increase the error to ±150 MeV. As a plausible model at a low scale µ 0 = 1 GeV we take [66] where the normalization constant N σ ≈ 1 can be determined in closed form. For heavylight mesons M ∼ (Qq) containing a heavy quark and a light anti-quark, one simply replaces x ↔ (1 − x) in the above relations.
In Table 2 we collect the values for the decay constants and the width parameters for heavy pseudoscalar and vector mesons, which will be used in our phenomenological analysis. In the cases of (qQ) and (QQ) bound states, Gegenbauer moments of roughly n 1/σ give important contributions to the LCDAs, because they are required to resolve the narrow structures of the LCDAs near the peak region. For example, at µ 0 = 1 GeV the first 5 (6) Gegenbauer coefficients of the B-meson (Υ-meson) LCDA are larger in magnitude than 0.1, and the first 7 (12) Gegenbauer coefficients are larger than 0.01. We will discuss in the next section that the effects of QCD evolution from a low scale µ 0 up to a high scale reduces the high-rank Gegenbauer moments much stronger than Gegenbauer moments of low rank. For example, at µ = m Z only the first 3 (2) Gegenbauer coefficients of the B-meson (Υ-meson) LCDA are larger than 0.1, and the first 6 (8) Gegenbauer coefficients are larger than 0.01. As a result, the shapes of the LCDAs for mesons containing heavy quarks are significantly affected by RG evolution. For the case of the B-meson LCDA, this effect was studied in [72]. Consequently, the low-scale predictions for the inverse moments considered here are strongly modified at µ = O(m Z ).

Radiative corrections and RG evolution
In order to improve the accuracy of our predictions and be in a position to meaningfully discuss the setting of the factorization scale µ, we include the O(α s ) radiative corrections to the leading-twist contributions in our analysis, finding that RG evolution effects are very important. The reason is that logarithms of the form α s ln(m 2 Z /µ 2 0 ) n , where µ 0 ∼ 1 GeV denotes the scale at which non-perturbative calculations of the LCDAs are performed, are numerically large and must be resummed to all orders of perturbation theory. We perform the calculation of the loop diagrams shown in Figure 2 using dimensional regularization with d = 4 − 2ǫ space-time dimensions. The individual on-shell graphs contain both UV and IR divergences. For the decays Z → Mγ and W → Mγ, which are mediated by vector and axial-vector currents, the UV divergences cancel in the sum of all diagrams. The remaining IR poles are cancelled when we renormalize the LCDAs. To this end, we express the bare LCDAs in terms of the renormalized ones, At one-loop order, one obtains [1,5] where is the one-loop Brodsky-Lepage kernel. For symmetric functions g(x, y), the plus distribution is defined to act on test functions f (x) as Besides the subtraction of 1/ǫ poles using dimensional regularization in the MS scheme, one must carefully address the question of how to define γ 5 in d = 4 dimensions. Some of the amplitudes considered in this work involve traces of Dirac matrices containing a single insertion of γ 5 . It is well known that for such traces the naive dimensional regularization scheme with anti-commuting γ 5 is algebraically inconsistent. Here we employ the 't Hooft-Veltman (HV) scheme [73], in which γ 5 = iγ 0 γ 1 γ 2 γ 3 anti-commutes with the four matrices γ µ with µ ∈ {0, 1, 2, 3}, while it commutes with the remaining (d − 4) Dirac matrices γ µ ⊥ . 3 While this definition is mathematically consistent, it violates the Ward identities of chiral gauge theories by finite terms, which must be restored order by order in perturbation theory [75]. In the present case, this is accomplished by performing the finite renormalization A µ = Z HV A µ HV of the axial-vector current, where [76] In addition, the leading-twist LCDA of a pseudoscalar meson, which is defined in terms of a matrix element of a non-local axial-vector current on the light-cone, receives a finite renormalization of the form where [77] This redefinition is important to restore the proper normalization of the LCDA φ P (x, µ).
Integrating relation (25) over x, we find that with Z HV given in (24). The integral turns the matrix element of the non-local axial-vector current into the corresponding local matrix element.
Our final expressions for the decay amplitudes will contain the scale-dependent, leadingtwist LCDAs φ M (x, µ) with M = P, V . These functions satisfy the integro-differential evolution equation where , and hence the Gegenbauer moments a n (µ) defined in (13) are multiplicatively renormalized at this order. They obey the RG equation [ where The evolution of the leading-twist LCDAs at two-loop order has been studied in [78][79][80][81]. The RG equation for the Gegenbauer moments becomes more complicated at this order, since the scale dependence of a M n (µ) receives contributions proportional to a M k (µ) with k = 0, . . . , n [80][81][82]. The evolution equation can still be solved analytically using an iterative scheme. Explicit results for the moments up to n = 12 can be found in [14]. However, given that all present estimates of the hadronic parameters a M n are afflicted with large theoretical uncertainties, it is sufficient for all practical purposes to use the leading-order solution (29). It reads where β 0 = 11 3 N c − 2 3 n f is the first coefficient of the QCD β function. Here µ 0 ∼ 1 GeV denotes a low scale, at which the Gegenbauer moments are derived from a non-perturbative approach, while µ is a high scale to which the LCDAs are evolved. In our analysis this scale is set by the mass of the decaying electroweak boson. Note that one must adjust the values of β 0 whenever µ crosses a flavor threshold. All of the anomalous dimensions are strictly positive, which implies that a M n (µ) → 0 in the formal limit µ → ∞. Indeed, for large n the evolution supplies an additional suppression factor (1/n) K with K = C F αs π ln µ 2 µ 2 0 . In this limit, the leading-twist LCDAs approach the asymptotic form 6x(1 − x). Figure 3 shows the RG evolution of the LCDAs of the kaon, J/ψ meson and B meson from a low scale µ 0 = 1 GeV up to a high scale m Z . We use the Gegenbauer moments and width parameters collected in Tables 1 and 2. For light mesons we truncate the Gegenbauer expansion (13) at n = 2. For heavy mesons we use the model LCDAs given in (17) and (19), compute their first 20 Gegenbauer moments, evolve the corresponding coefficients a M n from µ 0 to m Z , and reconstruct the LCDAs at the high scale from (13). The dotted line in the plots shows the asymptotic form 6x(1 − x). Evolution effects alter the shapes of the various distributions in a significant way. At the electroweak scale, the LCDAs are significantly closer to the asymptotic form 6x(1−x) than at a low hadronic scale. Consequently, RG effects render our predictions more insensitive to poorly determined hadronic input parameters. Notice, in particular, that the LCDA of the J/ψ meson at µ = m Z is as close to the asymptotic form as the kaon LCDA. In practice, the LCDAs of heavy mesons at a scale much larger than the heavy-quark mass can be well described in terms of a Gegenbauer expansion truncated after a few Gegenbauer moments.

Flavor wave functions of neutral mesons
The couplings of photons and of the electroweak gauge bosons W and Z to fermions are flavor dependent. While the flavor content of charged mesons is unambiguous, for neutral mesons complications arise from the fact that a given meson can be a superposition of different flavor components. We write the flavor wave function of the neutral final-state meson M 0 in the form For heavy mesons containing charm or bottom quarks such effects can safely be neglected. The heavy mesons η c and J/ψ have c c = 1, while η b and Υ have c b = 1. Mixing effects can however be important for light mesons. Following [32], we assume isospin symmetry of all hadronic matrix elements, but we differentiate between the matrix elements of mesons containing up or down quarks and those containing strange quarks. The π 0 and ρ 0 mesons are members of an isospin triplet and have flavor content (|uū − |dd )/ √ 2. Things get more complicated when we consider the mesons η, η ′ and ω, φ, however. In the SU(3) flavor-symmetry limit, the pseudoscalar meson η is 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 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 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 the two flavor states [83,84]. This assumption implies, in particular, that the vector mesons ω and φ are pure (|uū + |dd )/ √ 2 and |ss states, respectively, as is indeed the case to very good approximation. The anomaly introduces an effective mass term for the system of η and η ′ states, which is not diagonal in the flavor basis {|η q , |η s }. Since this is by assumption the only mixing effect, one obtains a mixing scheme with a single mixing angle in the flavor basis.
As explained in [32], the η and η ′ mesons have a leading-twist two-gluon LCDA besides the LCDAs corresponding to the quark-anti-quark Fock states η q and η s . The two-gluon LCDA contributes to the Z → η (′) γ decay amplitudes at order α s , through fermion box graphs with Zγgg as external particles. A detailed analysis of these decays will be presented elsewhere [33].

Radiative decays of electroweak gauge bosons
We now apply our general approach to study the rare, exclusive radiative decays Z → Mγ and W → Mγ, where M denotes a pseudoscalar (P ) or vector meson (V ). The leading-order Feynman diagrams contributing to the first process were already shown in Figure 1. We only consider cases where the mass of the final-state meson is much smaller than the mass of the decaying boson. Up to corrections of order (m M /m Z,W ) 2 this mass can then be set to zero.

Radiative hadronic decays of Z bosons
We begin our analysis with the decays Z 0 → M 0 γ. We find that, at leading order in the expansion in Λ QCD /m Z , only pseudoscalar or longitudinally polarized vector mesons can be produced. The corresponding decay amplitudes can be written in the general form where the upper (lower) sign refers to the case where M = P (V ). Here θ W is the electroweak mixing angle. Both the photon and the Z boson are transversely polarized with respect to the decay axis. The second term inside the brackets can be written more compactly as ε ⊥ Z ·ε ⊥ * γ , and below we use this as a short-hand notation. We use a convention where ǫ 0123 = −1. For neutral mesons that are eigenstates of the charge-conjugation operation, C invariance implies [24] The decay amplitudes are then proportional to the vector product ε Z × ε * γ of the transversely polarized photon and Z boson. However, in new-physics models in which the Z boson has flavor-changing neutral-current (FCNC) couplings, Z → Mγ decays into mesons that are not flavor diagonal (and hence not eigenstates of C) can occur. In this case relation (34) no longer holds. In complete generality, the decay rates, summed (averaged) over the polarization states of the photon (Z boson), are obtained as Here α = 1/137.036 is the fine-structure constant evaluated at q 2 = 0 [45], as appropriate for a real photon, and v denotes the Higgs vacuum expectation value, which enters through the relation (g/ cos where we have used α(m Z ) = 1/127.940 ± 0.014 and sin 2 θ W = 0.23126 ± 0.00005, with the weak mixing angle determined from the neutral-current couplings of the Z boson evaluated at µ = m Z [45]. The form factors F M i are given in terms of overlap integrals of calculable hard-scattering coefficients with LCDAs.
Evaluating the diagrams shown in Figures 1 and 2, we find that the relevant hard-scattering coefficients for the decays V → Mγ (with V = Z, W ) are given by where Our result for h + agrees with a corresponding expression derived in the context of a study of meson-photon transition form factors at high Q 2 performed in [85]. The expression for h − is new. The relevant convolutions of the hard-scattering coefficients H ± (x, m V , µ) with LCDAs give rise to the master integrals (we define a M 0 (µ) ≡ 1) with The integrals I M ± arise from the diagrams shown in Figure 2, in which the photon is attached to the quark inside the meson. Diagrams in which the photon is attached to the anti-quark give rise to the integralsĪ M ± . In evaluating the integrals we have used the Gegenbauer expansion (13). The two types of integrals are related to each other by the fact that the Gegenbauer polynomials C (3/2) n (2x − 1) transform into themselves times a factor (−1) n under the exchange of x ↔ (1 − x). Notice that at tree level the master integrals involve the infinite sums over Gegenbauer moments with equal coefficients. Employing a technique explained in Appendix C, we have succeeded to derive a closed expression for the one-loop coefficients c (n + 1)(n + 2) + 2 (n + 1) 2 (n + 2) 2 − 9 .  (30), it is straightforward to check that the master integrals in (39) are independent of the factorization scale µ. Indeed, the coefficient of the logarithm in (41) is equal to −γ n /(2C F ). Note also that the imaginary parts associated with the logarithm do not contribute to the decay rates at O(α s ).
In Figure 4, we study the scale dependence of individual terms in the sums over Gegenbauer moments in (39) at leading (dashed red lines) and next-to-leading order (solid lines) in perturbation theory. At leading order the µ dependence of the Gegenbauer moments, shown explicitly in (31), is left uncompensated, and hence a significant scale dependence arises. At next-to-leading order this dependence is compensated by the logarithmic terms contained in the one-loop corrections (41) to the hard-scattering coefficients C (+) n (m Z , µ) (blue lines) and C (−) n (m Z , µ) (orange lines). The resulting next-to-leading order curves exhibit excellent stability under variations of the factorization scale in the interval m Z /2 < µ < 2m Z .
In terms of the master integrals defined in (39), the form factors F M i are given by where and the coefficients Q ′ M and related to Q M by exchanging v q ↔ a q . Corrections to the results (42) arise only at twist-4 level and are suppressed by (Λ QCD /m Z ) 2 or (m M /m Z ) 2 . They are phenomenologically irrelevant. In the above expressions Q q denotes the electric charge of a quarks in units of e, while v q = 1 2 T q 3 − sin 2 θ W Q q and a q = 1 2 T q 3 (not to be confused with the Gegenbauer moments) are its vector and axial-vector couplings to the Z boson. Our finding that the form factors for pseudoscalar and vector mesons in (42) x), and hence for these mesons the odd Gegenbauer moments a M 2n+1 vanish. This leads to relation (34). The non-zero form factor F M 1 involves the infinite sum over the even Gegenbauer moments times some flavor-dependent coefficients Q M , which we collect in Table 3.
Explicit predictions for the leading-twist LCDAs derived by means of non-perturbative methods are typically obtained at a low hadronic scale µ 0 ∼ 1 GeV. When these predictions are used in (42), the expressions for the radiative corrections involve large logarithms ln(m 2 Z /µ 2 0 ) ≈ 9, which must be resummed to all orders in perturbation theory to obtain reliable predictions. This resummation is most readily performed by evaluating the result (42) We use the three-loop expression for the running coupling as provided by the RunDec program [86], normalized to α s (m Z ) = 0.1185 ± 0.0006 [45] and with heavy-quark thresholds at m b (m b ) = 4.163 GeV and m c (m c ) = 1.279 GeV [87]. The Gegenbauer moments at the high scale µ = m Z in the first line can be related to hadronic input parameters calculated at the low scale µ 0 = 1 GeV using the relations (31). In this process the coefficients of the higher moments get successively smaller. Decays into a transversely polarized vector meson are only allowed at twist-3 order. This presents us with an opportunity to study the structure of power corrections with a specific test case. We adopt the approximation where three-particle LCDAs are neglected. We then evaluate the diagrams in Figure 1 using the projector for a transversely polarized vector meson given in Appendix A. The decay amplitude can be decomposed in a form analogous to (33), The Z boson must be longitudinally polarized, and its polarization vector can be written as ε µ Z = (q − k) µ /m Z . The extra factor of m V /m Z compared with (33) makes the power suppression of these amplitudes explicit. The corresponding decay rate, summed (averaged) over the polarizations of the final-state (initial-state) particles, is given by The general expressions for the form factors F ⊥ i in terms of overlap integrals over the various twist-3 LCDAs appearing in the projector (A.4) for a transversely polarized vector meson are given in Appendix A. They can be simplified a lot by using relations implied by the equations of motion in the limit where three-particle LCDAs are neglected. Assuming for simplicity that quark-mass effects can be neglected, we obtain An analogous expression with Q V replaced by Q ′ V and a relative minus sign between the two terms inside the parenthesis holds for F ⊥ 2 . Since the leading-twist LCDAs of neutral mesons are symmetric in x ↔ (1 − x), it follows that F ⊥ 2 = 0, and hence once again only the term involving the Levi-Civita tensor in (45) contributes to the decay amplitude. Using the Gegenbauer expansion of the LCDA φ V (x, µ) in (13), we obtain a V 2n (µ) (n + 1)(2n + 1) .
Since we have not evaluated radiative corrections to the form factors, we do not control the scale dependence of the Gegenbauer moments a V n . However, it is clear that in order to avoid large logarithms we should again set µ ≈ m Z in the final result. Note that the result for F ⊥ 1 is similar to that for F V 1 in (42), but the coefficients of higher Gegenbauer moments are more strongly suppressed. For a rough estimate, we may assume that the Gegenbauer moments only have a minor impact on the final results (i.e. a V n (m Z ) ≫ 1), in which case it follows that the rates for Z → V γ decays with transversely and longitudinally polarized vector mesons are related by This ratio is of order 10 −4 . The outcome of this discussion is that power corrections in the expansion in Λ QCD /m Z are completely negligible for phenomenological applications. Figure 5: Non-local (left and center) and local (right) contributions to the W + → M + γ decay amplitudes.

Radiative hadronic decays of W bosons
The exclusive radiative decays W + → M + γ are, at first sight, very similar to the decays Z 0 → M 0 γ. Indeed, the contributions from the first two diagrams shown in Figure 5 can be obtained from the corresponding contributions to the Z-boson decay amplitudes by means of simple substitutions. The charged currents are flavor non-diagonal, and hence the final-state meson M + has a definite flavor structure described by a wave function |u idj . Note that now different electric-charge factors arise, depending on whether the photon is attached to the uptype quark or down-type anti-quark. The charged currents are purely left-handed, and hence we must replace in the equations of the previous section. However, a careful analysis shows that the first two diagrams in Figure 5 give rise to an extra contribution with a different tensor structure. It reads where the upper (lower) sign refers to the case of a pseudoscalar (longitudinally polarized vector) meson in the final state. Note that this contribution is independent of the LCDA of the final-state meson. It vanishes for an on-shell (transverse) photon, but is not compatible with U(1) em gauge invariance. Since the W boson has a direct coupling to the photon, an extra contribution to the W + → M + γ decay amplitudes exists, which arises from the third diagram in Figure 5, in which the final-state meson is produced by the conversion of an off-shell W boson. This graph has no analog in the Z-boson case. The corresponding contribution to the decay amplitude involves the meson matrix element of a local current, which to all orders in QCD is given in terms of a meson decay constant. We find where we keep the exact dependence on the vector-meson mass m M for the time being. The second relation can be simplified by considering the cases of longitudinal and transverse polarization separately. The polarization vector for a longitudinally polarized vector meson can be decomposed as which satisfies the conditions k · ε V = 0 and (ε V ) 2 = −1. The polarization vector for a transversely polarized vector meson is defined such that k · ε ⊥ V = q · ε ⊥ V = 0. We then obtain The second amplitude is non-zero only if the W boson is longitudinally polarized, and we have used a decomposition analogous to (53) in the final result. The local amplitudes for M = P, V are such that they combine with the extra term in (51) to give a gauge-invariant result proportional to ε ⊥ W · ε ⊥ * γ [24,25]. It follows from this discussion that, in analogy with (33), the leading-power amplitudes for the decays W + → M + γ can be written in the general form Summing (averaging) over the polarization states of the photon (W boson), we obtain the corresponding decay rates In close analogy with (42), we find that the form factors are given by The contribution −2 to F M 2 arises from the local contribution in (52). Corrections to these results are suppressed by (Λ QCD /m W ) 2 or (m M /m W ) 2 . The corresponding amplitudes for the decays W − → M − γ are obtained by replacing V ij → V * ij in (55) and by replacing the charge factors Q u ↔ Q d in (57). In addition, one must take into account that the odd Gegenbauer moments of the meson M − have the opposite sign as those of M + . This can be accounted for by replacing I M ± ↔Ī M ± . As a result, the form factor F M 1 remains invariant, while the form factor F M 2 changes sign. The decay rate in (56) stays invariant under these replacements. At low values of the factoriszation scale µ, the Wilson coefficients C (±) n (m W , µ) in (57) contain large logarithms of the form ln(m 2 W /µ 2 ), which can be resummed to all orders in perturbation theory by evaluating the scale-invariant quantities in (57) at the scale µ = m W . We obtain Decays into a transversely polarized vector meson are once again only allowed at twist-3 order. We decompose the decay amplitude in a form analogous to (45), such that The corresponding decay rate, summed (averaged) over the polarizations of the final-state (initial-state) particles, is given by The form factors can be calculated in analogy with the discussion in the previous section. The final results are where for simplicity we neglect contributions proportional to the quark masses. The local contribution in (54) adds −2(Q u − Q d ) to F ⊥ 2 . When expressed in terms of Gegenbauer moments, these results take the form In the limit where the Gegenbauer moments are neglected, we find in analogy with (49) that the ratio is strongly suppressed.

Absence of enhanced contributions from the axial anomaly
In has been suggested in the literature that the decay amplitudes for Z → P γ and W → P γ receive a very large enhancement due to an analog of the axial anomaly, which gives a contribution √ 2α/(πf π ) to the π 0 → γγ decay amplitude that does not vanish in the chiral limit [28,29]. We now explain why such a contribution does not exist in our case.
We consider the axial current A µ q =qγ µ γ 5 q and evaluate the triangle diagrams shown in Figure 6, where instead of two photons we take the external particles to be a photon and a Z boson. When taking the divergence of the current and considering the chiral limit m q → 0, we find that the loop diagrams vanish if one naively assumes that γ 5 anti-commutes with all Dirac matrices γ µ . However, a more careful regularization prescription adopting the HV scheme shows, like in the case of two external photons, that a finite remainder exists. It corresponds to a local operator for two gauge bosons. In operator language, we find that where for completeness we have included the terms proportional to the quark mass on the right-hand side. We do not show contributions involving two electroweak gauge bosons (ZZ or W + W − ) or two gluons, since they are irrelevant to our discussion. After electroweak symmetry breaking, the heavy Z boson acts like an external source, which is invariant under U(1) em gauge transformations. We see that anomalous contributions to the divergence of the axial-vector current not only involve the photon field, but also the Z boson. Indeed, such anomalous terms even arise for the charged currents A µ ij =d j γ µ γ 5 u i . Evaluating the corresponding triangle graphs with an external γW state, we obtain where we omit a contribution involving two electroweak gauge bosons (W Z). Does the existence of the anomalous di-boson terms in the above relations imply that there exist enhanced contributions to the Z → P γ and W → P γ decay amplitudes, in analogy to the famous case of the π 0 → γγ amplitude? This possibility was suggested in [28,29], where the authors speculated about a huge enhancement of the rates for the radiative decays Z → π 0 γ and W + → D s γ (see also the recent paper [15], which considers the decay W + → π − γ). We now demonstrate that such an enhanced contribution does not exist, focussing for concreteness on the case of a neutral pseudoscalar meson meson P . Let us parameterize a hypothetical anomalous contribution to the Z → P γ decay amplitude in the form where the general structure of the amplitude is consistent with (33). Let us furthermore parameterize the amplitude coupling an initial-state Z boson to the axial current A µ q and a photon as At lowest order, this amplitude is obtained from the diagram shown on the left in Figure 7. Inserting a complete set of hadron states that can be interpolated by the axial current A µ q , and summing over quark flavors, we find that the amplitude iM µαβ (k, q) contains the following contribution from the the single-hadron state P : This can be read off from the graph shown on the right in the figure. In the last step we have taken the chiral limit, in which the meson becomes massless. The key feature of the anomalous contribution would be that it exhibits a 1/k 2 pole in this limit. We have calculated the oneloop contributions to the amplitude M µαβ by evaluating the loop diagrams in Figure 7. We find that the contribution associated with the tensor structure shown in (68) is proportional to the expression (with p = k + q) where p 2 ≡ p 2 + i0 and k 2 ≡ k 2 + i0. For the case of the π 0 → γγ amplitude p 2 = 0 for the external photon, and in the limit ǫ → 0 one obtains −1/k 2 , which indeed exhibits a pole. Form the residue of this pole one can derive the anomaly-mediated π 0 → γγ decay amplitude. For
our case, on the other hand, p 2 = m 2 Z is equal to the mass of the decaying heavy gauge boson, in which case the above expression does not exhibit a 1/k 2 pole, but is instead proportional to 1/m 2 Z . Hence we conclude that A = 0 in (68). Note that in the limit k 2 → 0 one obtains from (69) 1 which is precisely of the form of our (bare) hard-scattering coefficients.

Phenomenological results
We are now ready to present detailed numerical predictions for the various radiative decay modes. We start with the decays of the Z boson, using relation (35). Besides the input parameters already mentioned, we need the Z-boson mass m Z = (91.1876 ± 0.0021) GeV and total width Γ Z = (2.4955±0.0009) GeV [45]. When squaring the decay amplitudes, we expand the resulting expressions consistently to first order in α s . The imaginary parts of the form factors in (42) do not enter at this order. Our results are presented in Table 4. Significant uncertainties in our predictions arise from the hadronic input parameters, in particular the meson decay constants (see Appendix B) and the various Gegenbauer moments. Their impact is explicitly shown in the table. Our error budget also includes a perturbative uncertainty, which we estimate by varying the factorization scale by a factor of 2 about the default value µ = m Z . All other uncertainties, such as those in the values of Standard Model parameters, are negligible. Note also that power corrections from higher-twist LCDAs are bound to be negligibly small, since they scale like (Λ QCD /m Z ) 2 for light mesons and at most like (m M /m Z ) 2 for heavy ones. The predicted branching fractions range from about 10 −11 for Z 0 → π 0 γ to about 10 −7 for Z 0 → J/ψ γ. In the last row, the symbol Υ(nS) means that we sum over the first three Υ states (n = 1, 2, 3). Strong, mode-specific differences arise foremost from the relevant flavor-dependent coefficients in Table 3, as well as from differences in the values of the decay constants. The combined uncertainties in the predictions for the branching fractions are typically of order 10% and are dominated by the uncertainties in the shapes of the LCDAs. The only exception are the decays Z 0 → Υγ, for which the relevant hadronic overlap integral is constrained by the model-independent relation (16). In our analysis we neglect two-loop QCD corrections, whose effects should be covered by the error we estimate from scale variations, and one-loop QED or electroweak radiative corrections, a few examples of which are shown in Figure 8. Their impact should be much smaller than the theoretical uncertainties inherent in our predictions. Consider, as a concrete example, the contribution of the first diagram, which only contributes to the Z → V γ amplitudes. Notice that the photon propagator 1/k 2 with k 2 = m 2 V is cancelled, because the Z → γγ * amplitude vanishes if both photons are on-shell [88]. As a result, there is no enhancement factor and the diagram is suppressed, compared with the leading contributions shown in Figure 1, by a factor α/π ∼ 2 · 10 −3 . This naive estimate is confirmed by the result of a detailed calculation of this contribution to the Z → J/ψ γ decay amplitude, which found that its effect leads to a reduction of the leading contribution by 0.2%, corresponding to a 0.4% correction of the branching ratio [27]. In the same paper, the authors have presented predictions for three of the Z → V γ decay modes along with theoretical error estimates. They are Br(Z → φγ) = (11.7 ± 0.8) · 10 −9 , Br(Z → J/ψ γ) = (9.96 ± 1.86) · 10 −8 , and Br(Z → Υ(1S) γ) = (4.93 ± 0.51) · 10 −8 . The last two branching ratios are consistent with our findings within errors. Note that in the NRQCD approach adopted by these authors the decay constants of the heavy quarkonia are themselves derived from an expansion about the non-relativistic limit. This introduces additional uncertainties, which can be avoided if the decay constants are extracted from data, as discussed in Appendix B. The analysis of the decay Z 0 → φγ presented in [27] uses an approach similar to ours but only includes the leading logarithmic evolution effects from the hadronic scale µ 0 = 1 GeV to the high scale µ ∼ m Z . Also, the value f φ = (235±5) MeV based on 2006 data is used for the φ-meson decay constant, which is inconsistent with our updated value f φ = (210 ± 5) MeV. Rescaling their result with the squared ratio of decay constants leads to Br(Z 0 → φγ) = (9.3 ± 0.6) · 10 −9 , which is consistent with our result but has a smaller uncertainty. The non-logarithmic O(α s ) corrections included here for the first time reduce the branching ratio by a significant amount. We also find that the present ignorance about the precise shape of the φ-meson LCDA gives rise to a larger theoretical uncertainty.
We now proceed to present our predictions for exclusive radiative decays of W bosons. In this case we need the input parameters m W = (80.385 ± 0.015) GeV and Γ W = (2.0897 ± 0.0008) GeV, as well as the relevant entries of the quark mixing matrix, which are |V ud | = 0.97425 ± 0.00022, |V us | = 0.2253 ± 0.0008, |V cs | = 0.986 ± 0.016, |V cd | = 0.225 ± 0.008, |V cb | = (41.1±1.3)·10 −3 , and |V ub | = (4.13±0.49)·10 −3 [45]. Starting from relation (56), we obtain the results shown in Table 5. In this case the pattern of the different decay modes reflects mainly the pattern of the relevant CKM matrix elements, and to a lesser extent the differences in the decay constants. The Cabibbo-allowed decays W → πγ, ργ, and D s γ have branching fractions of order few times 10 −9 to few times 10 −8 , where decays into heavy mesons are enhanced due to the structure of the relevant overlap integral in (18). The Cabibbo-suppressed modes W → K ( * ) γ and the strongly CKM-suppressed decay W → Bγ have correspondingly smaller branching ratios. The uncertainties inherited from CKM elements are shown where they are significant. In a recent paper, the W ± → π ± γ branching ratio was estimated to be 0.64 · 10 −9 [15], which is about 6.3 times smaller than the value we obtain (see below).
In the last two columns in Tables 4 and 5 we show different approximations to our results. The first one (labelled "asymptotic") gives the central values of the branching ratios (in the appropriate units) obtained if the asymptotic form 6x(1 − x) of the meson LCDA is employed. As we have explained, RG evolution effects from the low hadronic scale µ 0 = 1 GeV up to the electroweak scale have the effect of strongly suppressing the contributions from higher Gegenbauer moments. Indeed, we observe that using the asymptotic form provides reasonable approximations in most cases (especially for the Z → Mγ modes). The corresponding expressions for the decay rates read Decay mode Branching ratio SM background Table 6: Branching fractions for FCNC transitions Z → Mγ, which could arise from physics beyond the Standard Model. The different theoretical uncertainties have been added in quadrature. The last column shows our estimates for the irreducible Standard Model background up to which one can probe the flavor-changing couplings v ij and a ij . Here λ ≈ 0.2 is the Wolfenstein parameter.
where V ij is the relevant CKM matrix element for the production of the charged meson M + . The dominant corrections to the Z → Mγ branching fractions arise from the second Gegenbauer moment a M 2 , which is positive for light mesons and negative for heavy quarkonia. The dominant corrections to the W → Mγ branching fractions with kaons, D mesons or B mesons in the final state arise from the first Gegenbauer moment a M 1 . It gives a large positive contribution in all cases. In the case of heavy mesons this effect is particularly pronounced. The approximate results in (71) are fully consistent with corresponding (tree-level) expressions derived in [24]. The result for the Z 0 → π 0 γ decay rate derived in [25] is lower than ours by a factor 4/9, and the formula for the W ± → π ± γ decay rate derived in [15] differs from (71) by a factor 2/9. The origin of the discrepancy is related to the fact that the theoretical approach used in these papers is based on an expansion in a parameter ω 0 = 2, and the numerical estimates are obtained by keeping only the leading term in the expansion -a fact that was admitted in these papers. In Appendix D we trace the source of the discrepancy in more detail. In the last column in Tables 4 and 5 (labelled "LO") we present the branching ratios one would obtain at tree level using the model predictions for the LCDAs at the low scale µ 0 . In this approximation the one-loop QCD corrections, which contain large logarithms of the form α s ln(m 2 Z,W /µ 2 0 ), are omitted. For most decays the corresponding results overshoot the values obtained at next-to-leading order by significant amounts; only for decays into heavy quarkonia they underestimate the branching fractions.
Future precision measurements of the exclusive radiative decays Z → Mγ would serve as powerful tests of the Standard Model and of the framework of QCD factorization. The branching ratios for decays into vector mesons shown in Table 4 are proportional to |a q | 2 , where a q denote the axial-vector couplings of the Z boson to the various quarks. The couplings |a b | and |a c | have been measured at LEP with 1% accuracy [89], but no similarly accurate direct measurements of the couplings to light quarks are available. Given the theoretical precision of our predictions, it would be possible to determine these couplings with about 6% accuracy. The decays Z → Mγ can in principle also be used to search for non-standard FCNC couplings of the Z boson. If such couplings exist, then the diagrams shown in Figure 1 can lead to finalstate mesons of mixed flavor, such as K 0 , D 0 , B 0 and B s . It is straightforward to calculate (v bs ± a bs ) 2 < 5.5 · 10 −7 (v bs ) 2 − (a bs ) 2 < 1.4 · 10 −7 Table 7: Indirect constraints on the flavor-changing Z-boson couplings v ij and a ij (at 95% confidence level) derived from neutral-meson mixing [90][91][92].
the corresponding decay rates in our approach, starting from the general relations (35) and (42). We parameterize the non-standard vector and axial-vector couplings of the Z boson by v ij and a ij , respectively, where i, j are the quark flavors of the final-state meson. Our predictions for the corresponding branching fractions are given in Table 6. At higher order some Standard Model background to these searches exists, since electroweak loop graphs such as those shown in the last two diagrams in Figure 8 can give rise to flavor-changing transitions. Naive dimensional analysis shows that the contributions of these diagrams, relative to the contributions from the graphs in Figure 9 (in units of the new-physics couplings v ij and a ij ), scale like (α/π) |V ik V * kj |/ sin 2 θ W , where k can be any one of the three possible generation indices. The relevant loop functions depend on the dimensionless ratios m 2 k /m 2 Z and m 2 W /m 2 Z , which are either of O(1) or can be set to zero. Consequently one can only probe the newphysics couplings v ij and a ij up to some irreducible Standard Model background, which is estimated in the last column in Table 6.
Possible FCNC couplings of the Z boson are heavily constrained by precision flavor physics, in particular by bounds on the ∆F = 2 mixing amplitudes. It is a straightforward exercise to match our parameters v ij and a ij onto the Wilson coefficients in the general effective ∆F = 2 Hamiltonian as defined, e.g., in [90]. We obtain All other coefficients are zero at tree level. Using the bounds compiled in [90] as well as updated results reported in [91,92], we find the upper bounds on various combinations of v ij and a ij parameters shown in Table 7. The strongest bounds exist for the coefficients C 5 of mixed-chirality operators (right column). They can be avoided by assuming that v ij = ±a ij , such that the flavor-changing couplings are either purely left-handed or purely right-handed. Under this assumption one finds from the table that |v sd | < 8.5 · 10 −5 , |v cu | < 7.4 · 10 −5 , |v bd | < 1.0 · 10 −4 and |v bs | < 3.7 · 10 −4 , and the same bounds apply to |a ij |. If these indirect bounds are used, then the branching fraction shown in Table 6 are predicted to be at most a few times 10 −15 (a few times 10 −14 for the case of Z 0 → B s γ), meaning that they will be unobservable at the LHC and all currently discussed future facilities. We find it nevertheless worthwhile to illustrate the general idea of such new-physics searches. First of all, it should be emphasized that the indirect bounds derived from K −K, D−D and B d,s −B d,s mixing are to some extent model dependent, since one cannot tell whether the flavor violation originates from the couplings of the Z boson or from some other new particle. It is conceivable that in some (admittedly fine-tuned) models flavor-violating couplings of the Z boson can be compensated by the effects of some other, heavy boson. Also, in deriving the bounds on a particular Wilson coefficient C i one assumes that a single new-physics operator is present at a time and sets the coefficients of all other operators to zero. The method presented here, on the other hand, is unique in that it allows one (in principle) to probe for flavor-changing couplings of the Z boson directly and in a model-independent way, based on tree-level couplings of an on-shell particle. It should thus be seen as a complementary way to search for such effects. This method can also be generalized to the interesting case of flavor-changing exclusive Higgs-boson decays [19], for which the corresponding indirect bounds have been studied in [93]. While similar at first sight to the radiative Z-boson decays studied in Section 3, these decays are nevertheless interesting for several reasons. Unlike the photon, the final-state W boson can be longitudinally polarized, and hence several different helicity amplitudes contribute to the decay. Also, the trilinear ZW W coupling in the Standard Model gives rise to an additional contribution to the decay amplitude, in which the final-state meson is produced via the conversion of a W boson. This term is analogous to the "local" contribution we encountered in our study of the radiative decays of W bosons. Indeed, the leading-order Feynman graphs contributing to the Z 0 → M + W − decay amplitudes, shown in Figure 9, are analogous to those in Figure 5. Finally, and most interestingly, the decays Z → MW offer an opportunity to test the QCD factorization approach at a scale significantly lower than the Z-boson mass. A factorization theorem of the form (6) can only be derived if the momentum of the final-state meson in the rest frame of the decaying particle is much larger than its mass, since only then the constituents of this meson can be described in terms of collinear quark and gluon fields in SCET. The relevant condition is This condition is satisfied as long as m M ≪ We can thus use the factorization approach to calculate the branching fractions for Z → MW decays with a light, strange or charm meson in the final state, but not with a B meson. It also follows that the non-perturbative corrections to the factorization formula are organized in powers of The local contribution to the decay amplitudes involves the meson matrix elements of local currents. In close analogy with (52), we find Figure 9: Non-local (left and center) and local (right) contributions to the Z → M + W − decay amplitudes.
These results are exact even as far as the dependence on the vector-meson mass is concerned. The second relation can be simplified by considering the cases of longitudinal and transverse polarization of the vector meson separately. The longitudinal polarization vector can be written as Using this result, the second amplitude can be simplified to read while the decay amplitudes for the processes Z → V + ⊥ W − ⊥ and Z ⊥ → V + ⊥ W − , in which the final-state meson is transversely polarized, are suppressed relative to (76) by factors of m V /m Z and m V /m W , respectively. These power-suppressed amplitudes contribute to the decay rate at O(m 2 V /m 2 W ) and are thus negligible. They will be neglected below. Up to very small corrections, the two amplitudes in (74) thus have an identical structure.
Similar to (55), we can write the most general form of the decay amplitudes in the form The last structure inside the parenthesis contributes only if the W and Z bosons are longitudinally polarized. It has no analog in the case of radiative Z or W decays.
Summing (averaging) over the polarizations of the final-state (initial-state) particles, and setting the meson mass to zero, we obtain from (77) the decay rate Notice the close similarity of this result with (56). The differences are that we must now evaluate the coupling α at the electroweak scale, where for simplicity we do not differentiate between m Z and m W . Also, a phase-space suppression factor arises, and the third form factor F M 3 yields a contribution that is absent in (56). In deriving the above relation we have used that g 2 = e 2 / sin 2 θ W and sin 2 θ W = 1 − m 2 W /m 2 Z . Evaluating the first two Feynman graphs in Figure 9, we obtain Here r = m 2 W /m 2 Z = cos 2 θ W , and Z q = v q + a q are the left-handed couplings of quarks to the Z boson. Explicitly, we have Z u = 1 2 − 2 3 sin 2 θ W and Z d = − 1 2 + 1 3 sin 2 θ W . From the results shown in (74) and (76), we see that the local contribution from the third diagram in Figure 9 adds the term 2 to the form factors F M 2 . Note that the relevant combination entering the total decay rate is independent of the form of the LCDA. It follows that the third (longitudinal) term in the expression (78) for the decay rate yields 1/(2r). The integrals over the LCDA in our expressions for the form factors F M 1,2 can be evaluated analytically to any given order in the Gegenbauer expansion, but in practice it is easier to evaluate them numerically. Because the functions in the denominators are slowly varying with x, we find that higher-order Gegenbauer moments in the expansion (13) of the LCDA φ M (x, µ) give very small contributions. This is the essence of the approach by Manohar [25], which we illustrate in Appendix D. We obtain where the leading term 1.911 can be written in closed form as Decay mode  (79) is an interesting project, which we leave for future work. We note, however, that in this case the value of the factorization scale should be taken lower than the mass of the decaying Z boson. Two natural choices would be the typical momentum transfer, µ 2 ∼ 2k · q = (m 2 Z − m 2 W ) ≈ (43 GeV) 2 , and (twice) the energy of the final-state meson in the Z-boson rest frame, µ ∼ (m 2 Z − m 2 W )/m Z ≈ 20 GeV. In our analysis we shall use latter value, but since in the tree-level calculation presented here the scale dependence enters only through the tiny corrections involving a M 1,2 (µ) in (81) this choice has no noticeable effect.
Our predictions for the various Z → MW branching ratios are collected in Table 8. They are smaller than the corresponding branching fractions for W → Mγ decays by more than a factor 20 for light mesons and by more than a factor 60 for heavy mesons. Note, however, the curious fact that the Z 0 → π ± W ∓ branching ratio is about 15 times larger than the Z 0 → π 0 γ branching ratio. The latter quantity is tiny, because it is proportional to Q 2 π ≈ 7 · 10 −4 . We do not quote uncertainties related to the shapes of the meson LCDAs, which have a negligible impact on our result, see (81). We stress that to the quoted errors one should add an uncertainty accounting for our neglect of QCD radiative corrections. We expect that they can change the branching ratios by 10-20%.
It is interesting to compare our results with those obtained by Manohar in [25]. To this end, we used (36) to eliminate the parameter v and expand the function given in (82) in powers of the weak mixing angle. This yields (with s 2 W ≡ sin 2 θ W ) where we do not show terms of O(s 6 W ) and contributions from higher Gegenbauer moments. The leading contributions in this expression are in agreement with those found in [25], where an anlogous formula with the parenthesis replaced by 6 (1+c 2 W ) 2 = 3 2 + 3 2 s 2 W + 9 8 s 4 W + . . . is given.

Experimental considerations
Having obtained detailed and accurate theoretical predictions for a large class of very rare, exclusive radiative decays of Z and W bosons, we now address the question of how to search for such decays experimentally. While we do not perform an exhaustive feasibility study here, we discuss several ideas related to possible experimental analysis strategies. The goal is to get a feeling for how difficult it will be to observe some of the decay modes discussed in this work at present and future particle accelerators, and what accuracy one will be able to reach. We first consider the LHC, where about 10 11 Z bosons and 5 · 10 11 W bosons will have been produced in both ATLAS and CMS by the time the high-luminosity run with an anticipated integrated luminosity of 3000 fb −1 has been completed. We also consider a future lepton collider such as TLEP [16], where one can hope to produce about 10 12 Z bosons per year with a dedicated run on the Z pole, while 10 7 W -boson pairs per year would be produced in a run at the W W threshold. A run just above the tt threshold would also produce very large samples of W -boson pairs. At the LHC one needs to worry about triggering and reconstruction. The trigger for the decays Z → Mγ can be based on photons and muons. The energy of the final-state photon is comparable to the energy of the photons produced in the Higgs-boson decay h → γγ, where it has been demonstrated that such events can be triggered on. Muons, on the other hand, are produced only in some cases, in particular in the decay modes containing a vector meson.
Reconstructing these event appears to be challenging but not impossible. Probably the most promising modes are Z → J/ψ γ and Z → Υ(nS) γ followed by a fully leptonic decay of the heavy quarkonium state. The corresponding rates, however, are very small. If only the muon channel can be used, then the combined Z → V γ → µ + µ − γ branching fractions are of order 5 · 10 −9 for J/ψ and 1.5 · 10 −9 for Υ(1S). Thus, we can expect to trigger on several hundred Z → J/ψ γ → µ + µ − γ events and up to one hundred Z → Υ(1S) γ → µ + µ − γ events at the LHC. While there is no significant physics background we can think of, the combinatorial background may be substantial. Given these challenges, it is encouraging that ATLAS has recently reported first upper bounds (at 95% CL) on the branching fractions Br(Z → J/ψ γ) < 2.6 · 10 −6 , Br(Z → Υ(1S) γ) < 3.4 · 10 −6 , Br(Z → Υ(2S) γ) < 6.5 · 10 −6 , and Br(Z → Υ(3S) γ) < 5.4·10 −6 [23]. Further dedicated experimental studies of these decays would be worthwhile. The yield of Z → Υγ events can be enhanced by about a factor 2 if one combines the Υ(nS) channels with n = 1, 2, 3, which have similar leptonic branching fractions (ranging between 1.9% and 2.5%). The case Z → Υ(4S) γ may be particularly interesting, since the Υ(4S) resonance can decay to a pair of B mesons, which gives rise to displaced vertices. It might be possible to achieve a larger effective rate for this decay mode by using highly efficient b-tagging methods. Observing Z → Mγ decays into other final states seems to be difficult at the LHC. Ideas for reconstructing highly energetic φ, ρ and ω mesons have been presented in [19] in the context of a study of exclusive h → V γ decays. For light mesons decaying into two photons, such as π 0 and η, it might be possible to tell that there are more than one photon in the final state provided one of the photons is converted into an e + e − pair.
The situation with W → Mγ decays at the LHC seems less promising. In this case the final state contains a charged hadron, and it is not clear to us how one could reconstruct it in the high-multiplicity environment of a hadron collider. Perhaps the most promising case is the decay W → D s γ, of which over 10,000 events should be produced. The problem is how to tag the D s mesons. A very interesting, dedicated study of several exclusive radiative W → Mγ decays has been performed in [15], to which we refer the reader for more details.
It looks much more promising to search for rare exclusive Z → Mγ decays at a future lepton collider such as TLEP, in particular if one envisions a dedicated high-luminosity run on the Z pole. The advantage of such a "Z factory" is that it would produce yields of up to 10 12 Z bosons per year in a very clean environment. The detectors of such a future experiment are expected to have excellent particle-identification systems. This would make it possible to perform precision studies of many of the decay modes we have discussed. In particular, all Z → Mγ decays except for Z → π 0 γ should be accessible, and in several cases it should be possible to measure the branching ratios at the percent level. There is even hope for observing some of the weak radiative decays Z → MW , whose branching fractions can be of order several times 10 −10 . With dedicated runs above the W W or tt thresholds one would produce samples of W bosons large enough to search for the Cabibbo-allowed W → Mγ decay modes. It would be most rewarding to perform detailed feasibility studies for measurements of rare exclusive Z and W decays at such a facility.

Summary and conclusions
Based on the formalism of QCD factorization, we have performed a comprehensive and systematic analysis of the very rare, exclusive radiative decays of Z and W bosons into final states containing a single pseudoscalar or vector meson. The basis of our study is a factorization theorem derived in soft-collinear effective theory, which expresses the decay amplitudes as convolutions of calculable hard-scattering kernels with light-cone distribution amplitudes (LCDAs), in a systematic expansion in powers of (Λ QCD /m Z,W ) 2 and (m M /m Z,W ) 2 . For the first time, we have included the complete set of one-loop QCD radiative corrections. Large logarithms involving the ratio of the electroweak scale to the typical hadronic scale µ 0 , at which model predictions for the LCDAs are obtained, are resummed to all orders of perturbation theory. We have also estimated the leading power-suppressed effects, finding that their impact on the branching ratios is typically of order 10 −4 compared with the leading terms. Larger power corrections up of to 1% can only arise for mesons containing heavy b quarks. The exclusive decays Z → Mγ and W → Mγ therefore offer an ideal laboratory for testing the QCD factorization approach in a controlled and theoretically clean way. Our main phenomenological results for the relevant branching ratios are collected in Table 9. We have not considered in this work the interesting decays Z → η (′) γ. Their analysis is complicated by the fact that the η and η ′ mesons contain a flavor-singlet Fock component, and to predict the decay rates at O(α s ) one needs to take into account the two-gluon LCDA of these mesons. This will be discussed in a separate publication [33].
Our results form the basis of a rich, novel program of electroweak precision studies of the Standard Model, which offers powerful tests of the QCD factorization approach in a situation where it should deliver precise predictions. With the statistics obtainable in the high-luminosity run at the LHC and, in a much cleaner environment, at future lepton colliders, it will be possible to observe several of the very rare decays discussed here and measure their Decay mode
Several generalizations and extensions of our work are possible and worth exploring. Our formalism can be applied in a straightforward way to obtain high-precision predictions for exclusive radiative (and weak radiative) decays of the Higgs boson, extending previous treelevel analyses presented in [17][18][19][20][21][22]. One goal of such studies is to search for enhanced Yukawa couplings and flavor-changing interactions of the Higgs boson. In this context, new-physics studies analogous to those presented in Section 3.4 are particularly interesting. Without further conceptual developments, our formalism can also be extended to calculate the rates for purely hadronic decays, such as Z, W, h → M 1 M 2 or even decays with more than two particles in the final state. These extensions are left for future work.
The physics case for studying some of the very rare, exclusive decays of heavy electroweak bosons is compelling to us. There is some beautiful physics to be explored here, both from the theoretical and the experimental points of view. We hope that our detailed exploratory survey will raise sufficient interest that some dedicated feasibility studies for discovering such decays at the high-luminosity LHC and future lepton colliders will be performed.
In this approximation, the twist-3 two-particle amplitudes can be expressed in terms of the twist-2 LCDA φ ⊥ V .
The light-cone projector for a transversely polarized vector meson is yet more complicated. Up to twist-3 order, one obtains ⊥ (y, µ) ∂ ∂k ⊥µ + 3-particle LCDAs . (A.4) Note that, compared with [42], the terms multiplying the Levi-Civita tensor in the second line have the opposite sign, because we are using a different sign convention for this object. In the approximation where three-particle LCDAs are neglected, the equations of motion yield the relations [40][41][42] g (v) ⊥ (x, µ) WWA = 1 2 x 0 dy φ V (y, µ) x 0 dy φ V (y, µ) − g In this approximation, the twist-3 two-particle amplitudes can be expressed in terms of the twist-2 LCDA φ V . As an application of these results, we present the general expressions for the form factors F ⊥ i entering the Z → V ⊥ γ and W + → V + ⊥ γ decay amplitudes in (45) and (59). In the first case, we obtain (withx ≡ 1 − x, and ignoring quark-mass effects for simplicity) ⊥ (x) 4 ⊥ (x) 4 ⊥ (y) .

(A.6)
We omit the scale dependence of the various quantities for simplicity. In the case of (59) we must omit the prefactors Q (′) V /6 and instead assign charge factors Q u and Q d to the two terms under each integral over x. These results can be simplified significantly by using the relations (A.5). The final expressions have been given in (47) and (61).

B Determinations of meson decay constants
We follow [46] and determine the relevant meson decay constants from experimental data. The decay constants of charged pseudoscalar mesons can be determined from their semileptonic decays P − → l −ν l (γ). This analysis is performed by the Particle Data Group [45], and it leads to the values for f π , f K , f D and f Ds shown in Tables 1 and 2.
The decay constants of light charged mesons can also be obtained from the one-prong hadronic decays of the τ lepton. The corresponding decay rates are given by For our analysis the value ( 3 n=1 f 2 Υ(nS) ) 1/2 = (930 ± 4) MeV is also needed. As discussed in [69], the combined effects of ρ -ω and ω -φ mixing lower f ω by about 7.5 MeV and f φ by about 13 MeV, while they have a negligible impact on f ρ . Note also that there are some uncertainties related to the sizable width of the ρ resonance, which we have ignored here. They might explain why the value of f ρ 0 extracted here is larger than the value of f ρ − extracted above from τ decay. Combining the various extractions and using conservative error estimates, we obtain the values shown in Table 1.

C Gegenbauer expansion of the convolution integrals
In order to derive closed analytic expressions for the basic convolution integrals I M ± andĪ M ± defined in (39), we employ the definition of the Gegenbauer polynomials in terms of the generating function 1 n (x) t n . (C.1) From (13) and the second relation in (39), it follows that the coefficients C (±) n (m V , µ) defined in (40) are the coefficients of t n of the integrals where h(t) = 2 t 2 (C. 3) The coefficients c (±) n (m V /µ) in (40) are the coefficients of t n in the series expansion of h(t) around t = 0. It is then a straightforward exercise to derive expression (41).

D Connection with the approach by Manohar
The approach put forward by Manohar in [25] was originally developed as a tool to study the exclusive decay Z → W π by means of a local operator-product expansion. It corresponds to the following series expansion of the propagator in the first diagram in Figure 9: Using this expansion, the convolution integrals over the pion LCDA in (79) can be expressed in terms of local operator matrix elements, and the leading terms are determined by the normalization of the LCDA. For the phenomenological estimates in [25] and [15] only the leading term was kept. When the same approach is used for radiative decays, the W -boson mass in the above relations must be set to zero, in which case one obtains Keeping the leading term is now unjustified, but if it is done this corresponds to replacing 1/x → 2 under the convolution integrals in (39). In the asymptotic limit, where all Gegenbauer moments are set to zero, these integrals are therefore too small by a factor 2/3. This explains the discrepancies in the predictions for the Z → Mπ and W → Mπ branching fractions mentioned in Section 3.4.