Exclusive Radiative Higgs Decays as Probes of Light-Quark Yukawa Couplings

We present a detailed analysis of the rare exclusive Higgs-boson decays into a single vector meson and a photon and investigate the possibility of using these processes to probe the light-quark Yukawa couplings. We work with an effective Lagrangian with modified Higgs couplings to account for possible new-physics effects in a model-independent way. The h->V\gamma{} decay rate is governed by the destructive interference of two amplitudes, one of which involves the Higgs coupling to the quark anti-quark pair inside the vector meson. We derive this amplitude at next-to-leading order in \alpha_s using QCD factorization, including the resummation of large logarithmic corrections and accounting for the effects of flavor mixing. The high factorization scale \mu~m_h ensures that our results are rather insensitive to the hadronic parameters characterizing the light-cone distribution amplitude of the vector meson. The second amplitude arises from the loop-induced effective h\gamma\gamma* and h\gamma Z* couplings, where the off-shell gauge boson converts into the vector meson. We devise a strategy to eliminate theoretical uncertainties related to this amplitude to almost arbitrary precision. This opens up the possibility to probe for O(1) modifications of the c- and b-quark Yukawa couplings and O(30) modifications of the s-quark Yukawa coupling in the high-luminosity LHC run. In particular, we show that measurements of the ratios Br(h->\Upsilon(nS)\gamma)/Br(h->\gamma\gamma) and Br(h->bb)/Br(h->\gamma\gamma) can provide complementary information on the real and imaginary parts of the b-quark Yukawa coupling. More accurate measurements would be possible at a future 100 TeV proton-proton collider.


Introduction
The discovery of the Higgs boson in 2012 [1,2] has established the existence of a new kind of elementary particle, which couples to the other particles of the Standard Model (SM) in a non-universal way. The SM predictions that the Higgs couplings to heavy gauge bosons and fermions are given by 2m 2 W,Z /v and m f /v, where v ≈ 246 GeV is the Higgs vacuum expectation value, have been confirmed within experimental uncertainties for the W and Z bosons and for the third-generation fermions [3,4]. However, no direct measurements of the Higgs couplings to the light fermions of the first two generations are available at present. It is not difficult to come up with models in which these couplings can deviate significantly from those predicted in the SM. For example, in [5,6] it was proposed that the Yukawa couplings may depend in a non-trivial way on the Higgs field, and that this might explain the hierarchies seen in the spectrum of fermion masses. A more general analysis of different classes of models in which the Higgs couplings to fermions can differ significantly from those of the SM was presented in [7]. Probing the Higgs couplings to light fermions is thus of paramount importance. This includes both flavor-diagonal and flavor-changing interactions. Correlations between the two types of couplings, which to some extent are model dependent, have been studied in [8]. The CMS collaboration has recently reported a slight excess in the search for the flavor-violating decay h → τ ± µ ∓ [9], which, if interpreted as a signal, corresponds to a branching fraction Br(h → τ ± µ ∓ ) = (0.89 + 0. 40 − 0.37 ) %. Not surprisingly, this observation has led to much theoretical speculation.
Our focus in the present work is on the Higgs couplings to light quarks (q = t). In a couple of beautiful papers, it has recently been proposed that one might get access to these couplings by focussing on the rare, exclusive decays h → V γ of the Higgs boson [10,11], where the final state contains a single vector meson V . Such measurements are extremely challenging at the LHC, as the corresponding branching fractions are in the range of few times 10 −6 . Nevertheless, observing these processes is not hopeless in view of the fact that in its highluminosity run the LHC will serve as a Higgs factory. With 3 ab −1 of integrated luminosity, about 1.7×10 8 Higgs bosons per experiment will have been produced [12], and even larger Higgs samples could be obtained at a future facility such as a 100 TeV proton-proton collider. The theoretical description of rare exclusive decays employs the formalism of QCD factorization [13][14][15][16][17], which was originally developed for the analysis of hard exclusive QCD processes and later extended to the more complicated case of exclusive hadronic two-body decays of B mesons [18,19]. In a recent paper, we have systematically developed this approach for the case of the exclusive decays Z → Mγ, W → Mγ and Z → MW [20]. These processes are less susceptible to new physics and can therefore be used to test the QCD factorization approach and extract valuable information about the light-cone distribution amplitudes (LCDAs) of various mesons M. Observing exclusive radiative decays of heavy electroweak gauge bosons in the high-luminosity run of the LHC would provide a proof of principle that this kind of raredecay searches can be performed in a hadron-collider environment. A promising first step in the direction of observing the decays Z → J/ψ γ and Z → Υγ, along with the corresponding Higgs decays h → J/ψ γ and h → Υγ, has recently been reported by the ATLAS collaboration [21].
In this work we extend previous studies of exclusive radiative decays of the Higgs boson in several important ways. We include QCD radiative corrections and resum large logarithms arising from evolution effects between the Higgs mass scale and a low hadronic scale. We comment on the structure of power-suppressed corrections to the factorization formula, study the effects of the off-shellness of the photon in the h → γγ * → γV process and include the power-suppressed h → γZ * → γV contribution. We also take into account the effects of ω−φ mixing. Most importantly, our analysis allows for generic non-standard Higgs couplings to SM fermions and gauge bosons, including CP-odd interactions. We devise a strategy which allows us to eliminate the theoretical uncertainties related to the dominant h → γγ * → γV conversion contribution, including possible new-physics effects, to almost arbitrary precision. This is a crucial prerequisite for achieving the desired sensitivity to the Yukawa couplings of light quarks. Finally, we address all relevant h → V γ decays to both light and heavy mesons in one coherent formalism. For technical details on the QCD factorization approach the reader is referred to [20].
In the following section we start by defining effective couplings of the Higgs boson to SM fermions and gauge bosons. These include both CP-even and CP-odd couplings. In Section 3 we discuss in detail the calculation of the different contributions to the h → V γ decay amplitudes in the QCD factorization approach, distinguishing between "direct" contributions induced by the Yukawa couplings of the Higgs boson to light quarks (q = t) and "indirect contributions" resulting from a h → γγ * /γZ * transition followed by the conversion of the off-shell boson into a vector meson. We show that uncertainties related to the effective hγγ coupling strength, as regards to both theoretical uncertainties and possible new-physics contributions, can be eliminated by studying the ratio of the h → V γ and h → γγ branching fractions. In Section 4 we present a phenomenological analysis both in the SM and in generic new-physics models with modified Higgs interactions. In particular, we point out that measurements of the ratios of the h → Υ(nS) γ, h → bb and h → γγ branching fractions can yield highly complementary information on the real and imaginary parts of the bottom-quark Yukawa coupling. Section 5 contains a summary of our results and the conclusions. Technical details are relegated to five appendices.

Effective Higgs couplings
In our analysis we assume SM couplings for all particles other than the Higgs boson. For the Higgs couplings to SM quarks and gauge bosons we adopt the effective Lagrangian where s W ≡ sin θ W and c W ≡ cos θ W are the sine and cosine of the weak mixing angle. We use s 2 W = 0.23126 ± 0.00005 as determined from the neutral-current couplings of the Z boson evaluated at µ = m Z [22].F µν = 1 2 ǫ µναβ F αβ is the dual field-strength tensor, and we use a sign convention where ǫ 0123 = 1. Our choice of factoring out a loop factor in the second line is made for later convenience. For SM extensions in which the new particles are heavy, the coefficients Figure 1: Two-loo Barr-Zee diagram (left) and effective one-loop contribution (right) to the EDM of the electron arising from the CP-odd couplingsκ f ,κ γγ andκ γZ in the effective Lagrangian (1). of these higher-dimensional operators are suppressed by two powers of the new-physics scale. We emphasize that the above is not a complete list of operators. For instance, we have not included higher-dimensional operators of the form hW µν W µν and hW µνW µν , whose coefficients are already strongly constrained by data. These operators would enter our analysis only via the h → γγ * /γZ * one-loop amplitudes, and without loss of generality their effects can be absorbed into the coefficients κ γγ , κ γZ andκ γγ ,κ γZ .
Both the CP-even couplings κ i and the CP-odd coefficientsκ i are real. In the SM κ W = κ Z = κ f = 1, whileκ f and all the remaining couplings in the second line vanish. Our κ q and κ q parameters for quarks are related to the corresponding Yukawa couplings by In the last step we have introduced rescaled parameters normalized to the mass of the b quark. This will turn out to be a useful definition for the quarks of the first two generations. For our analysis we need the quark masses at the scale µ = m h , where m h = (125.09 ± 0.24) GeV is the Higgs-boson mass [23]. We define the running quark masses at next-tonext-to-leading order (NNLO) in the MS scheme, starting from the low-energy values given in [22]. This yields m b (m h ) = (2.79 ± 0.02) GeV, m c (m h ) = (622 ± 12) MeV, m s (m h ) = (52.8 ± 1.4) MeV, m d (m h ) = (2.66 ± 0.11) MeV, and m u (m h ) = (1.21 ± 0.08) MeV. For the top quark we obtain m t (m h ) = (166.8 ± 0.7) GeV starting from m t (m t ) = (163.4 ± 0.7) GeV, which we have derived from the present world average obtained in [24] using the conversion tables provided in [25]. Present LHC data are largely insensitive to the Yukawa couplings to light quark flavors. From a global χ 2 fit to the measured Higgs rates, the authors of [11,26] have derived the bounds |κ u | 2 + |κ u | 2 < 1.3 and |κ d,s,c | 2 + |κ d,s,c | 2 < 1.4 at 95% confidence level (CL). The corresponding bounds for the original parameters are |κ u | 2 + |κ u | 2 < 3000, Bounds on the CP-violating Higgs couplings to third-generation fermions have been studied in [27]. Under the assumption that the Higgs couples to the electron in the standard way (κ e = 1,κ e = 0), the strongest constraints come from the electric dipole moment (EDM) of the electron. They arise from two-loop Barr-Zee diagrams such as the one shown on the left in Figure 1 and yield |κ t | < 0.01, |κ b | < 1.9 and |κ τ | < 2.4 at 90% CL. Bounds from the neutron EDM are approximately one order of magnitude weaker, but they do not reply on assumptions about the Higgs couplings to the electron. The h → bb and h → τ + τ − rate measurements at the LHC can be used to place upper limits on the combinations |κ b,τ | 2 + |κ b,τ | 2 , which imply the stronger bounds (at 95% CL) |κ b | < 1.44 and |κ τ | < 1.24 from CMS [3], and |κ b | < 1.3 and |κ τ | < 1.5 from ATLAS [4]. Bounds from the electron EDM can also be derived for the local operators multiplied byκ γγ andκ γZ in (1). Evaluating the one-loop contributions to the electron EDM shown on the right in Figure 1, we obtain in the MS scheme where x Z = m 2 Z /m 2 h . This contribution is logarithmically UV divergent, because the inner structure of the effective hγγ and hγZ vertices is not resolved. The term proportional toκ γγ κ e agrees with a calculation performed in [28]. The subtraction scale µ should be identified with the scale of the new physics, which is responsible for these non-standard interactions. Setting µ = Λ NP = 1 TeV for an estimate, and assuming SM-like Higgs couplings to the electron, we obtain from the present experimental bound |d e | < 8.7 · 10 −29 e cm (at 90% CL) [29] the constraint κ γγ + 0.09κ γZ < 0.006 (90% CL) .
Barring a fine tuning of the two contributions, this implies that |κ γγ | < 6 · 10 −3 and |κ γZ | < 0.07. If the new-physics scale lies above 1 TeV then these bounds become stronger. With Λ NP = 10 TeV, for example, they improve by a factor of 2.

Radiative hadronic decays of Higgs bosons
Our focus in this work is on the rare, exclusive radiative decays h → V γ, where V denotes a vector meson with momentum k and the photon carries momentum q. We will refer to vectors orthogonal to the plane spanned by k and q as being transverse. Up to tiny corrections of order (m V /m h ) 2 , the mass of the vector meson can be set to zero. The leading-order Feynman diagrams are shown in Figure 2. In the first two graphs, the Higgs boson couples to the quark and anti-quark pair from which the meson is formed. We refer to this as the "direct" contribution to the decay amplitude [30,31]. It can be calculated from first principles using the QCD factorization approach [13][14][15][16][17], because the energy released to the final-state meson is much larger than the scale of long-distance hadronic physics [11,32,33]. At leading power in an expansion in Λ QCD /m h , the direct contribution can be expressed as a convolution of a calculable hard-scattering coefficient with the leading-twist LCDA of the vector meson V . The corresponding factorization formula was derived in [20] using the formalism of soft-collinear effective theory [34][35][36][37]. It was shown that, for a given helicity amplitude, the power corrections to the leading term are suppressed by (Λ QCD /m h ) 2 for light mesons and (m Q /m h ) 2 for mesons containing heavy quarks of flavor Q. Even for the b-quark these power corrections are negligible. The third diagram in Figure 2 shows a different production mechanism, in which the vector meson is produced via the conversion of an off-shell photon or Z boson produced in a h → γγ * /γZ * transition [10]. We refer to this as the "indirect" contribution. It involves the hadronic matrix element of a local current and thus can be expressed in terms of the decay constant f V of the vector meson. The direct contribution is sensitive to the Yukawa coupling of the Higgs boson to the quarks which make up the vector meson. We shall find that in the SM the direct and indirect contributions to the h → V γ decay amplitude interfere destructively. They are of similar size for V = Υ, while the direct contributions are smaller than the indirect ones by factors of about 0.06 for V = J/ψ, 0.002 for V = φ, and few times 10 −5 for V = ρ 0 and ω. The sensitivity to the Yukawa couplings thus crucially relies on the precision with which the indirect contributions can be calculated. We will come back to this point below.
The most general parametrization of the h → V γ decay amplitude is where both the final-state meson and the photon are transversely polarized. From (5), the decay rate is obtained as Here α = 1/137.036 is the fine-structure constant evaluated at q 2 = 0 [22], as appropriate for a real photon. We choose to normalize the decay amplitude in (5) to the vector-meson decay constant f V , which is defined in terms of a matrix element of a local vector current. Since we consider neutral, flavor-diagonal mesons, the definition of the decay constants (and of other hadronic matrix elements) is complicated by the effects of flavor mixing. In complete generality, such a neutral meson V can be regarded as a superposition of flavor states |qq . We can thus define flavor-dependent decay constants f q V via A certain combination of these flavor-specific decay constants can be measured in the leptonic decay V → e + e − . The corresponding decay amplitude involves the matrix element of the electromagnetic current J µ em = q Q qq γ µ q. We thus define Here c V q denote the flavor coefficients in the naive constituent-quark model, where and so on. With these definitions, the electromagnetic decay rate is given by In [20], the so-defined decay constants f V have been extracted using the most recent experimental data. Our motivation for factoring out these decay constants in (5) is that they can be determined experimentally without adopting any specific model for flavor mixing. As we will discuss in Section 3.2, the dominant terms in the indirect contributions to the h → V γ decay amplitudes are proportional to precisely these quantities.

Direct contributions to the form factors
We first consider the direct contribution to the decay amplitude shown in the first two diagrams in Figure 2, in which the Higgs boson couples to the quarks contained inside the vector meson. In order to calculate the corresponding contributions to the form factors F V i defined in (5), one calculates the corresponding partonic amplitudes with on-shell quark and anti-quark states and then projects these amplitudes onto the leading-twist LCDA of a transversely polarized vector meson V . The relevant light-cone projector reads [38] where the dots stand for higher-twist contributions. In analogy with (7) and (8), we define a set of flavor-specific transverse decay constants via The quantity f ⊥ V entering (11) is then defined by the combination The transverse decay constants are scale-dependent quantities, since the QCD tensor current has a non-zero anomalous dimension. The leading-twist LCDA φ ⊥ V (x, µ) can be interpreted as the amplitude for finding a quark with longitudinal momentum fraction x inside the meson. 1 It depends on the choice of the factorization scale µ employed in the factorization formula.
1 Strictly speaking, one should introduce flavor-specific LCDAs φ q⊥ V (x, µ) and define the product f q⊥ V (µ) φ q⊥ V (x, µ) in terms of a matrix element of a non-local quark current with flavor q, in analogy with relation (4) in [20]. Because the LCDAs are normalized to 1, and given the present large uncertainties in the shapes of these functions, it is a safe approximation to employ SU(3) symmetry and replace φ q⊥ The two leading-order graphs shown in Figure 2 are supplemented by the diagrams in Figure 3, which arise at O(α s ). Including these effects is crucial in order to control the scale dependence of the transverse decay constant f ⊥ V , the Yukawa coupling y q and the LCDA φ ⊥ V . It will also allow us to resum large logarithms of the form α s ln(m 2 h /µ 2 0 ) n to all orders in perturbation theory. Here µ 0 ≈ 1 GeV is a typical hadronic scale, at which model predictions for f ⊥ V and φ ⊥ V are obtained. We work in dimensional regularization and subtract UV and IR divergences in the MS scheme. The product of the bare decay constant times the LCDA of a transversely polarized vector meson is related to the product of the corresponding renormalized quantities via where at one-loop order with C F = 4/3. The relevant one-loop Brodsky-Lepage kernel reads [13,17] For the decays h → V γ, which are mediated by (pseudo-)scalar currents, an overall UV divergence remains, which is cancelled by the counterterm for the Yukawa coupling, derived from y q,bare = µ ǫ Z y (µ) y q (µ) with When dealing with pseudo-scalar currents we employ the 't Hooft-Veltman (HV) scheme [39], 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 γ µ ⊥ . 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 [40]. This is accomplished by performing the finite renormalization P = Z P HV P HV of the pseudo-scalar current P =qγ 5 q, where [41] Z P By evaluating the relevant Feynman graphs in Figures 2 and 3, we find that the direct contributions to the form factors in the amplitude decomposition (5) are given by where we have defined The reduced form factors F V are given by with This result agrees with a previous calculation performed in [32] apart from a typo. 2 We have expressed the Yukawa coupling in terms of the running b-quark mass using the second relation in (2). We focus primarily on the cases V = J/ψ and Υ(nS), where to an excellent approximation the vector meson contains a single quark flavor q, and henceκ J/ψ ≈κ c and κ Υ(nS) ≈ κ b , and similarly for the CP-odd parametersκ V . For the light mesons V = ρ 0 , ω and φ, on the other hand, flavor-mixing effects can be important. This concerns, in particular, the possibility of a small admixture of an |ss flavor component in the wave functions of ρ 0 and ω, which can be important due to the smallness of the Yukawa couplings to the up and down quarks in the SM. Since the ρ 0 meson is a member of an isospin triplet, its flavor mixing into |ss can only be caused by electromagnetic interactions or isospin-violating effects in QCD. Both types of effects are estimated to be very small, and hence we expect that |f s⊥ ρ 0 /f ⊥ ρ 0 | ≪ 1. To good approximation we can thus use the naive relation The situation is different for the case of the ω meson, whose mixing into an |ss flavor state can be non-negligible. In Appendix A we derive explicit expressions for the parametersκ ω andκ φ in a simple flavor-mixing scheme for the ω−φ system. Assuming that |κ s | ≫ |κ u,d | like in the SM, and working in the SU(3) limit and to first order in the small mixing angle θ ωφ , we obtainκ 2 These authors use the pole mass instead of the running quark mass in the prefactor, which adds We are grateful to the authors for confirming this mistake.
In the SM the contributions from the up and down quarks almost precisely cancel inκ ω , and hence the contribution induced by ω − φ mixing is likely to be the dominant one. Existing estimates for the mixing angle θ ωφ derived from mass-independent analyses include θ ωφ ≈ 0.05 [42] and θ ωφ ≈ 0.06 [43,44]. On the other hand, in a more recent mass-dependent analysis the values θ ωφ (m 2 ω ) ≈ 0.008 and θ ωφ (m 2 φ ) ≈ 0.081 were obtained [45]. We conclude from this discussion thatκ φ ≈κ s to good approximation, while a more accurate description of flavormixing effects would be required before the quantityκ ω can be interpreted reliably in terms of quark Yukawa couplings.
In the factorization formula (20) all non-perturbative hadronic physics is contained in the decay constants and the LCDA. The quantity multiplying the LCDA under the integral is the hard-scattering coefficient, which can be calculated in perturbation theory. It depends on the momentum distribution of the quark inside the hadronic bound state. QCD-based model calculations of LCDAs are typically performed at a low hadronic scale µ 0 ∼ 1 GeV. If such a low value is chosen for the factorization scale µ in (20), the hard-scattering coefficient contains large logarithms of the form α s ln(m 2 h /µ 2 0 ) n with ln(m 2 h /µ 2 0 ) ≈ 9.7, which should be resummed to all orders in perturbation theory. To perform this resummation, it is convenient to use the expansion of the LCDA in the basis of Gegenbauer polynomials, which reads [13,17] We can then rewrite (20) as a sum over Gegenbauer moments, finding with The factor in brackets in (25) precisely compensates the scale dependence of the product m b (µ) f ⊥ V (µ), while the quantity I V is formally scale invariant. Using a technique developed in [20], we obtain for the hard-scattering coefficients C n in moment space the closed expression where H n+1 = n+1 and a V ⊥ n (µ) up to that scale. The solution of the corresponding renormalization-group (RG) equations at NLO of RG-improved perturbation theory is discussed in Appendix B. At leading order in QCD the Gegenbauer moments a V ⊥ n (µ) in (24) are Table 1: Hadronic input parameters for light vector mesons. The decay constants f V are extracted from data on the electromagnetic decay widths V → l + l − [20], while the ratios f ⊥ V /f V are derived from a compilation of theoretical predictions. The values of the Gegenbauer moments at the scale µ 0 = 1 GeV are taken from [48,49]. The last two columns show the effective charges Q V and v V defined in (8) and below (32).
renormalized multiplicatively, such that [13,30] Here µ 0 denotes a low reference scale, while µ = O(m h ) is the hard factorization scale, to which the LCDAs are evolved. The NLO corrections to these relations are discussed in Appendix B.
They have a negligible impact on our numerical results. All of the anomalous dimensions are strictly positive (for n = 0), which implies that a V ⊥ n (µ) → 0 in the formal limit µ → ∞. In this limit the leading-twist LCDAs approach the asymptotic form 6x(1 − x). Similarly, it follows from relation (B.2) in Appendix B that the transverse decay constants of vector mesons vanish in the asymptotic limit, i.e. f ⊥ V (µ) → 0 for µ → ∞. It has been emphasized in [20] that RG evolution effects render our predictions rather insensitive to the precise values of the Gegenbauer moments. From (26), we obtain While all Gegenbauer moments have O(1) coefficients at the high-energy scale, the coefficients of the higher moments are strongly reduced when one expresses the answer in terms of moments normalized at the low scale µ 0 = 1 GeV. In order to obtain numerical predictions for the reduced form factors we need as hadronic input parameters the decay constants f V and f ⊥ V and the Gegenbauer moments a V ⊥ 2n of the various vector mesons. As mentioned earlier, the decay constants f V can be extracted from experimental data, and up-to-date values have been derived in [20]. The ratios f ⊥ V /f V needed in (25) must be obtained using some non-perturbative approach, such as lattice QCD, lightcone QCD sum rules or the non-relativistic effective theory NRQCD for heavy quarkonia [46,47], which 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. Details of such determinations are reviewed in Appendix C. In Table 1 we compile the relevant input Table 2: Hadronic input parameters for heavy quarkonium states. The decay constants f V are extracted from data on the electromagnetic decay widths V → l + l − [20], while the ratios f ⊥ V /f V are derived from NRQCD scaling relations. The width parameters σ V are obtained from relation (30), where the first error is of parametric origin and the second one parameterizes the uncertainty due to higher-order effects. The last two columns show the effective charges Q V and v V . parameters for light vector mesons. Because of the lack of information about higher Gegenbauer moments we can only keep few terms in the infinite sum (24). The systematics of the Gegenbauer expansion has been discussed in [20], where it was pointed out that the higher moments a V ⊥ n with n ≫ 1 fall off faster than 1/n. Indeed, high-rank Gegenbauer polynomials C (3/2) n (2x − 1) with n ≫ 1 would resolve structures on scales ∆x ∼ 1/n. For a light vector meson V , it is reasonable to assume that the LCDA φ ⊥ V (x) does not exhibit pronounced structures at scales much smaller than O(1). To estimate the impact of higher moments we use a V ⊥ 4 (µ 0 ) = ±0.15 for our error estimates. Relation (29) suggests that the effect of yet higher-order terms is small.
The LCDAs of heavy mesons exhibit a different behavior, since the presence of the heavyquark mass introduces a new scale. For a quarkonium state V ∼ (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 NRQCD matrix element called v 2 V [50]. Including the one-loop QCD corrections calculated in [32], we obtain A critical discussion of the extraction of the parameters v 2 V for different quarkonium states is presented in Appendix C. Using the values compiled there, but with increased error estimates, we obtain the ratios of decay constants and the width parameters σ V (µ 0 ) at the low scale µ 0 = 1 GeV shown in Table 2. As a reasonable model at the scale µ 0 we adopt the form [20] where the polynomial in front of the Gaussian factor ensures that the LCDA vanishes at the endpoints x = 0, 1. In order to estimate the uncertainties related to the functional form and to capture the effects of unknown higher-order corrections to relation (30)  error of ±25% on the σ V parameters. Given this form, we compute the first 20 Gegenbauer moments at the low scale µ 0 , evolve them up to the factorization scale µ ≈ m h using (28), and use these results in evaluating relation (26).
We are now ready to present our numerical results for the direct contributions to the reduced form factors F V in (25) for various vector mesons, including detailed error estimates. They are collected in Table 3. The different sources of theoretical errors contain a perturbative uncertainty (subscript "µ"), which we determine by varying the factorization scale µ between m h /2 and 2m h . Once the NLO corrections are included our results are very stable under scale variations. The scale uncertainties are larger for the imaginary parts than for the real parts of the form factors, since these start at O(α s ) and there is thus no compensation of the scale dependence. We emphasize, however, that the imaginary parts only have a small impact on our numerical predictions for the decay rates. We also include the uncertainty in the value of the b-quark mass, which has a very small impact. The uncertainties related to hadronic parameters include the ratio f ⊥ V /f V (subscript "f ") and uncertainties in the shapes of the LCDAs, as modeled by the values of the Gegenbauer moments a V ⊥ 2 and a V ⊥ 4 (light mesons) and the width parameter σ V (heavy mesons). These hadronic uncertainties are the dominant sources of errors. The last column in the table shows the results obtained when all errors are added in quadrature. These numbers will be used for our phenomenological analysis in Section 4. We observe that the spread of the results for the form factors F V for different vector mesons is rather small. The theoretical uncertainties on the real part of F V are typically between 4% and 9%. The only exception is F J/ψ , for which the uncertainty in the ratio of decay constants is about 16%. It would probably be possible to reduce this uncertainty by performing a more detailed NRQCD analysis. Power-suppressed corrections to our results (20) and (25) can be organized in an expansion in (Λ QCD /m h ) 2 for light mesons and (m V /m h ) 2 for heavy mesons [20]. They are at most of order 10 −2 for V = Υ(nS), 10 −3 for V = J/ψ, and 10 −4 for all light mesons. It is thus safe to work with the leading-order terms. 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 the top row in Figure 4. For flavor-diagonal final-state mesons, the first diagram involves a factor m W v α π ∼ 0.7 · 10 −3 instead of the Yukawa coupling y q in the diagrams in Figure 2, while the second diagram involves a factor y t α π ∼ 2 · 10 −3 . This is smaller (by roughly a factor 10) than the charm-quark Yukawa coupling and of the same order as the strange-quark Yukawa coupling. If the goal is to reach sensitivity to the strange-quark Yukawa coupling in the SM, then these electroweak corrections should be calculated. However, such a level of sensitivity will be out of reach at the LHC. In the bottom row in Figure 4 we show other examples of neglected diagrams, which involve a h → gg(γ) transition followed by the conversion of the two gluons into the final-state meson. The diagram on the left corresponds to a two-gluon LCDA of the vector meson, which does not exist at leading twist due to the Landau-Yang theorem. The graph on the right is analogous to the second diagram shown in the first row, but with the internal W bosons replaced by gluons. A naive estimate indicates that the two types of effects should be of similar magnitude.

Indirect contributions to the form factors
We now proceed to study the photon-and Z-pole contributions to the decay amplitude shown in the third diagram in Figure 2. The crossed circle in this diagram represents the off-shell h → γγ * and h → γZ * amplitudes, which in the SM are induced by loop graphs involving a virtual charged fermion or a W boson (in unitary gauge), as shown in Figure 5. Since the = Figure 5: One-loop SM contributions to the effective hγγ and hγZ vertices.
indirect contributions to the decay amplitudes are numerically dominant over the direct ones and our goal is to compute them with the highest possible accuracy, we include the effect that the intermediate gauge boson is slightly off shell (k 2 = m 2 V ), and we keep the full dependence on the meson mass even though this is a very small effect.
The exact one-loop expressions for the off-shell h → γγ * and h → γZ * amplitudes have been derived in [51]. Using these results and extending them to the case of the pseudo-scalar Higgs couplings in (1), we find where x V = m 2 V /m 2 h accounts for the effects of the off-shell boson, are the vector couplings of the Z boson to fermions. It is a safe approximation to neglect flavor-mixing effects for the subleading contribution from h → γZ * → γV . At one-loop order, the loop functions are given by The first two terms in each coefficient are the contributions from the quarks and leptons, the third term in C γγ and C γZ arises from gauge-boson loops, and the last term accounts for possible new-physics contributions parameterized by the operators shown in the second line of (1). We have introduced the dimensionless variables τ f = 4m 2 f /m 2 h (for f = q, l) and τ W = 4m 2 W /m 2 h . We use the running quark masses m q (m h ) when evaluating the variables τ q , which is appropriate in view of the large momentum transfer in the loop. Explicit expressions for the loop functions A f , A γV W and B f are given in Appendix D. In the SM we have κ q = κ l = κ W = 1 and κ γγ = κ γZ = 0. The effective Higgs couplingsκ i entering in (34) all vanish in the SM.
Since for small values of τ f the fermion loop functions A f (τ f , x V ) and B f (τ f , x V ) are proportional to τ f , it suffices for all practical purposes to keep the contributions from the thirdgeneration fermions. The effects of the off-shellness of the photon that converts into the finalstate vector meson gives rise to very small corrections. The relevant variable x V = m 2 V /m 2 h varies between 3.8 × 10 −5 for V = ρ and 6.9 × 10 −3 for Υ(3S). Note also that the contribution of the h → γZ * amplitude in (32) is by itself strongly power suppressed. Numerically, we obtain and similar expressions hold for the CP-odd coefficientsC γγ andC γZ . Notice that the contributions from second-generation fermions are very small even if one assumes that their Yukawa couplings are as large as those of the b-quark (i.e., forκ c,s = 1). In the SM, we have C γγ (0) = −3.266 + 0.021i and C γZ (0) = −2.046 + 0.005i, while the CP-odd coefficientsC γγ andC γZ vanish. In the expressions for the form factors in (32) these coefficients are weighted by different factors. Note also that to good approximation the indirect contributions to the form factors are proportional to 1/m V and hence they are larger for lighter mesons.
Let us now discuss the impact of QCD and electroweak radiative corrections on the above results. Gluon exchanges between the two quark lines that make up the final-state meson in the third diagram in Figure 2 are of non-perturbative nature and are accounted for by the meson decay constant f V . Radiative corrections connecting the effective hγγ and hγZ vertices in the last diagram in Figure 2 to the final-state quarks give rise to graphs which no longer receive the enhancement proportional to m 2 h /m 2 V from the photon propagator. Also, at least two gluons would need to be exchanged by color conservation. The effect of such diagrams is negligible. We thus only need to worry about QCD corrections to the h → γγ * and h → γZ * amplitudes, which arise from two-loop graphs in which a gluon is exchanged inside the quark loop. These corrections have been calculated in numerical form in [52] and analytically in [53]. In practice the corrections are only relevant for the top-quark contribution. Their effect is to enhance the decay amplitude by a few percent. Two-loop electroweak corrections to the h → γγ amplitude in the SM were calculated in [54,55] (see also [56]). One finds that they are small and negative. For m h = 125.09 GeV, the effects of QCD and electroweak corrections nearly cancel each other, leaving a total correction of about −0.2% [55]. In our phenomenological analysis we will ignore radiative corrections on the indirect contributions computed here, as their combined effect is well below the 1% level. We will device a strategy where the dominant contribution proportional to C γγ (0), including radiative corrections, drops out. Radiative corrections then only have a tiny impact of the coefficients of the subleading terms.

Reduction of theoretical uncertainties
If we assume that the Higgs couplings to the electron are close to those predicted by the SM, the CP-odd form factors F V 2 are more than two orders of magnitude smaller than F V 1 . For a first discussion we can thus focus on the form factors F V 1 . Keeping only the numerically significant terms, we find where the last term in each line represents the direct contribution. We have dropped the small imaginary parts of the latter, whose impact is tiny. Replacingκ φ ≈κ s ≈ 0.019κ s and κ c ≈ 0.223κ c one obtains equivalent expressions in which the modified Higgs couplings are expressed as corrections to the SM Yukawa couplings. The challenge is to detect the small impact of the direct contributions in the last two cases.
To this end, it is essential to have absolute confidence in the precision with which the indirect contributions can be calculated in the SM, and to be able to subtract these contributions in a reliable way without assuming that the SM is correct. The latter task can be accomplished because the off-shellness of the photon in the h → γγ * contribution as well as the h → γZ * contribution in the third graph in Figure 2 are both very small effects. It is therefore possible to eliminate the main dependence of the indirect contributions on the new-physics parameters by considering the following ratio of decay rates: (37) Taking such a ratio has the additional advantage that one becomes insensitive to the unknown total width of the Higgs boson, and hence one obtains directly the ratio of branching fractions. One can even go one step further and eliminate the sensitivity to the decay constant f V by using (10) and considering the ratio The only remaining hadronic uncertainties are now contained in the calculation of the reduced form factors F V , which we have collected in Table 3. The explicit expressions for the various quantities entering the right-hand side of (37) are r CP =C γγ (0)/C γγ (0) and where we work to leading order in the small ratios m 2 V /m 2 Z and x V = m 2 V /m 2 h . Since the individual terms in these expressions are all normalized to C γγ (0), it is convenient to define an effective parameter κ eff γγ by normalizing this coefficient to its SM value. Specifically, we write Here and below the omitted contributions have coefficients equal to zero up to the indicated number of digits. From the empirical fact that the Higgs couplings to W bosons and to the third-generation fermions agree with their SM values within errors, it follows that κ eff γγ cannot be too different from its SM value, except perhaps for a possible new-physics contribution parameterized by κ γγ . Note that a tiny imaginary part of κ eff γγ , which can be present for non-standard values of the κ i parameters, would have no noticeable impact on our analysis. The ratio r CP entering in (37) vanishes in the SM and is entirely due to the various CP-odd new-physics parametersκ i . We find Under the assumption that the Higgs couplings to the electron are approximately SM like, the upper bounds |κ t | < 0.01 and |κ γγ | < 0.006 mentioned above relation (3) imply that the first and last term in the numerator, which have the largest coefficients by far, are at most 0.004 and 0.002 in magnitude, respectively. It follows that the magnitude of r CP can at most be of order 1%, and hence the impact of this parameter in (39) is likely to be negligible. We emphasize that the parameter r CP can in principle be probed experimentally by studying h → γγ decays, in which both photons undergo nuclear conversion to electron-positron pairs [57]. In practice, such a measurement appears to be very challenging. We now present our numerical results for the CP-even coefficients ∆ V for the various mesons. The complete expressions are collected in Appendix E. They contain direct contributions proportional to the relevantκ q parameters and indirect contributions, which are due to the power-suppressed h → γZ * → γV contribution and the effect of the off-shellness of the photon in the h → γγ * → γV contribution. These latter terms are significantly smaller than the theoretical uncertainties in the direct terms. In our phenomenological analysis we will keep these small effects but evaluate them in the SM. This gives rise to the expressions ∆ ρ 0 = (0.068 ± 0.006) + i(0.011 ± 0.002) κ ρ 0 κ eff γγ + 0.00002 , and Approximate expressions forκ ρ 0 ,κ ω andκ φ have been given in (22) and (23). The constant terms in the above results show the tiny power-suppressed corrections. Only for the Υ(nS) states they reach the level of percent. Our complete expressions for the CP-odd coefficients∆ V are also given in Appendix E. It is a good approximation to only keep the direct contributions in these terms, which are likely to give rise to the dominant effects. Their coefficients are the same as in the expressions above, but withκ q replaced byκ q and κ b replaced byκ b . It is interesting to compare our result for the quantities ∆ V with corresponding expressions obtained by other authors. From [11] one can extract ∆ ρ 0 = (0.095 ± 0.020) (2κ u +κ d )/3, ∆ ω = (0.092 ± 0.021) (2κ u +κ d ) and ∆ φ = (0.130 ± 0.027)κ s , while from [33] one can obtain ∆ J/ψ = (0.392±0.053)κ c , ∆ Υ(1S) = (1.048±0.046)κ b , ∆ Υ(2S) = (1.138±0.053)κ b and ∆ Υ(3S) = (1.175 ± 0.056)κ b . These values are systematically higher than ours due to the fact that these authors have not (or not fully) included QCD radiative corrections and RG evolution effects in the direct contributions. For the Υ(nS) states it is important to keep the small imaginary parts of the direct contributions, since in the SM the real parts almost perfectly cancel in the combinations 1 −∆ V in (37). The result for ∆ ω obtained in [11] misses the contribution from ω−φ mixing and contains a sign mistake in front ofκ d . Note also that our predictions for the ∆ V parameters of light mesons are significantly more accurate than those obtained in [11].

Phenomenological results
We begin by quoting our benchmark results for the h → V γ branching fractions in the SM. For a Higgs mass of m h = (125.09 ± 0.024) GeV, the SM value of the h → γγ branching ratio is (2.28 ± 0.11) · 10 −3 [58]. Using this result, we obtain for the decays into light vector mesons Br(h → ρ 0 γ) = (1.68 ± 0.02 fρ ± 0.08 h→γγ ) · 10 −5 , Br(h → ωγ) = (1.48 ± 0.03 fω ± 0.07 h→γγ ) · 10 −6 , where we quote separately the uncertainties due to the vector-meson decay constant f V and the h → γγ branching ratio, the latter being the dominant source of uncertainty. Our predictions are systematically lower and more accurate than those obtained in [11], where the values Br(h → ρ 0 γ) = (1.9 ± 0.15) · 10 −5 , Br(h → ωγ) = (1.6 ± 0.17) · 10 −6 and Br(h → φγ) = (3.0 ± 0.13) · 10 −6 are quoted. While the first two results are compatible with ours within errors, there is a significant difference for the important mode h → φγ. For decays into heavy vector mesons, we find In these cases there is an extra source of theoretical uncertainty related to the calculation of the direct contribution to the decay amplitude. Note that there is an almost perfect cancellation between the direct and indirect contributions to the h → Υ(nS) γ decay amplitudes, and as a consequence the resulting branching ratios are roughly three orders of magnitude smaller than the h → J/ψ γ branching fraction. For comparison, we note that the branching ratios found in [33] read (2.79 + 0. 16 − 0.15 ) · 10 −6 for J/ψ, (0.61 + 1.74 − 0.61 ) · 10 −9 for Υ(1S), (2.02 + 1.86 − 1.28 ) · 10 −9 for Υ(2S) and (2.44 + 1.75 − 1.30 ) · 10 −9 for Υ(3S). We find good agreement with the results reported by these authors except for the decay h → Υ(1S) γ, where their value is about a factor 7 smaller than ours. The reason is that we do not neglect the imaginary part of the direct contribution to ∆ Υ(1S) in (42), which prevents 1 − ∆ Υ(1S) 2 from becoming arbitrarily small.
Our predictions may also be compared with the upper limits obtained from a recent first analysis of these rare decays reported by the ATLAS collaboration. They are [21]. It will require an improvement by a factor 500 to become sensitive to the h → J/ψ γ mode in the SM, while the SM branching fractions for the decays h → Υ(nS) γ are out of reach at the LHC. Nevertheless, as we will discuss below, these decay modes allow for very interesting new-physics searches. With 3 ab −1 of integrated luminosity, about 1.7 × 10 8 Higgs bosons per experiment will have been produced by the end of the highluminosity LHC run [12]. If the J/ψ is reconstructed via its leptonic decays into muon pairs, the effective branching ratio in the SM is Br(h → J/ψ γ → µ + µ − γ) = 1.8 · 10 −7 , meaning that about 30 events can be expected per experiment. If also the decays into e + e − can be used, then ATLAS and CMS can hope to collect a combined sample of about 120 events. A detailed discussion of the experimental prospects for reconstructing these events over the background can be found in [10]. Concerning the h → φγ decay mode, a reconstruction efficiency ǫ φγ = 0.75 was assumed for the φγ final state in [11], which appears to us as an optimistic assumption. In the SM one expects about 400ǫ φγ events per experiment in this mode, meaning that the two experiments can hope to look at a combined sample of several hundred events. Likewise, in the SM one expects about 2900ǫ ρ 0 γ events per experiment in the decay mode h → ρ 0 γ.
In Figure 6 we show our predictions for the ratio of branching fractions (times 1000) defined in (37) in the plane of the parametersκ V /κ eff γγ andκ V /κ eff γγ . We focus on the most interesting cases V = φ, J/ψ and Υ(1S). The corresponding plots for V = ρ 0 , ω would look very similar to that for V = φ (apart from the overall scale of the branching fractions), while the plots for higher Υ(nS) resonances would look very similar to that for the Υ(1S) meson. For orientation, we mention that a value of 0.4 in these plots corresponds to a h → V γ branching fraction of about 10 −6 , assuming that the h → γγ branching fraction is SM like. This assumption will be implicit whenever we quote absolute branching ratios below; otherwise the quoted numbers must be rescaled by Br(h → γγ)/Br(h → γγ) SM . The structure of our results (42) implies that the rescaled Yukawa couplings always enter normalized to κ eff γγ . Hence, if a deviation from the SM is observed in any of these modes, then this could be caused either by a new-physics effect on the h → γγ branching ratio (parameterized by κ eff γγ ) or by a non-standard Yukawa coupling. The former effect would however be correlated among all decay channels. We observe that the h → φγ branching ratio is rather insensitive to the CP-odd parameterκ φ , the reason being that this parameter enters quadratically and only via the direct contribution, which by itself is a small correction. An analogous statement holds (with less accuracy) for the case h → J/ψ γ. It is thus a reasonable approximation to study these decay modes under the To an excellent approximation, λ 2 bγ measures the deviation of the ratio Br(h → bb)/Br(h → γγ) from its SM value. The Higgs bosons must be produced via the same production mechanism in both cases, so that possible new-physics effects in Higgs production cancel out. Since the h → bb mode is measured at the LHC in the rare V h and tth associated-production channels, at present no accurate direct measurements of λ bγ are available. However, from the model-independent global analyses of Higgs couplings performed by ATLAS and CMS, in which all couplings to SM particles (including the effective couplings to photons and gluons) are rescaled by corresponding κ i parameters and also invisible Higgs decays are allowed, one obtains λ bγ = 0.63 ± 0.27 for CMS (see Figure 17 in [3]) and λ bγ = 0.67 ± 0.32 for ATLAS (see Table 9 in [4]). At 95% CL this (roughly) implies λ bγ < 1.3. Within this allowed region, the h → Υ(1S) branching ratio varies by more than two orders of magnitude and can take values as large as 1.3 · 10 −6 . This might be accessible in the high-luminosity run of the LHC. If the Υ(1S) meson is reconstructed via its decays into muon or electron pairs, one could then hope for a sample of about 20 events with 3 ab −1 combining the ATLAS and CMS data sets.
In order to better assess the theoretical uncertainties in our predictions, we now study the projections of the results onto the axis where the CP-odd couplings vanish. For the light mesons (V = ρ 0 , ω, φ), settingκ V = 0 has basically no impact on the branching ratios. In Figure 7 we show the ratio of branching fractions defined in (37) as a function of the CP-even couplingsκ V for h → ρ 0 γ and h → φγ. The width of the bands indicates the theoretical uncertainty. We have not included the small uncertainty in the values of the decay constants f V , because they can be eliminated using relation (38), and we assume that by the time the h → V γ modes will be discovered the decay constants will have been measured more precisely than today. The corresponding plot for h → ωγ would look identical to the left plot, but with a different vertical scale. While the theoretical uncertainties are small in all cases, we observe that the sensitivity of the branching ratios to the modified Yukawa couplings is unfortunately rather week. For example, a hypothetical 10% measurement of the h → ρ 0 γ branching ratio at the SM value would imply that |κ ρ 0 /κ eff γγ | < 0.8, which is to say that a certain combination of the up-and down-quark Yukawa couplings is bounded not to exceed 80% of the b-quark Yukawa coupling. A 1% measurement would be required to obtain the more interesting bound |κ ρ 0 /κ eff γγ | < 0.08, which is still more than 100 times the SM value forκ ρ 0 given in (22). The situation is not much better for the case h → φγ. With a 10% measurement of the branching fraction at the SM rate, one would be able to conclude that −0.55 <κ φ /κ eff γγ < 0.62. With a 1% measurement one would obtain the bounds −0.04 <κ φ /κ eff γγ < 0.08, which would come close to the SM valueκ φ ≈ 0.02. Such a measurement is however out of the reach of the LHC.
We now turn to the more interesting cases of radiative Higgs decays into heavy quarkonium states. In Figure 8 we show our predictions as a function of the physical parameters κ c (not κ c ) and κ b , again assuming that the CP-odd couplingsκ c andκ b vanish. In the latter case the impact of a possible CP-odd coupling on the branching fraction can be significant, and in the case of a measurement of a non-standard rate one should keep this possibility in mind. From the left plot in the figure we conclude that a 20% measurement of the h → J/ψ γ branching ratio at the SM value would allow one to constrain −0.51 < κ c /κ eff γγ < 3.07, which would provide quite interesting information on the CP-even charm-quark Yukawa coupling. With a 10% measurement this range could be shrunk to 0.32 < κ c /κ eff γγ < 1.53, and with a 5% measurement one could reach 0.75 < κ c /κ eff γγ < 1.19. Such accurate measurements serve as an interesting physics target for a future 100 TeV proton-proton collider. The situation with the h → Υ(nS) γ decay modes is different and quite interesting. In the SM the corresponding branching fractions shown in (45) are so small that these decays would be unobservable. The strong suppression arises from an almost perfect cancellation between the direct and indirect contributions to the decay amplitudes, which results from the fact that in the SM Re ∆ Υ(nS) ≈ 1 within a few percent, see (43). Thanks to this fortuitous fact, these decays offer a much enhanced sensitivity to the effects of new physics. For instance, the SM value of the h → Υ(1S) γ branching ratio of 4 × 10 −9 can be enhanced by a factor of more than 200 for κ b /κ eff γγ ≈ −1 or κ b /κ eff γγ ≈ 3. The first of these possibilities would yield a h → bb rate consistent with current LHC measurements. For example, with a hypothetical 25% measurement Br(h → Υ(1S) γ)/Br(h → γγ) = (0.4 ± 0.1) · 10 −3 one would conclude from the figure that −1.21 < κ b /κ eff γγ < −0.64, which would be a very significant piece of information and a spectacular sign of new physics.
One may ask whether the current bounds obtained by the ATLAS collaboration already have a significant impact on the Higgs couplings. Unfortunately this is not the case. We find that the upper values reported in [21] imply approximately both at 95% CL. For comparison, we note that m t /m c ≈ 268 and m t /m b ≈ 60. In other words, the current bounds derived from exclusive h → V γ decays imply that the couplings of the charm and bottom quarks to the Higgs boson are not much stronger than the top-quark Yukawa coupling (the more optimistic value |κ c | < 220 was quoted in [26]). We emphasize again that any experimental information on the rare radiative decays h → Υ(nS) γ should be interpreted in terms of an allowed region in the two-dimensional plane of the couplings κ b /κ eff γγ andκ b /κ eff γγ . The one-dimensional projection shown in Figure 8 is meant for illustrative purposes only. It is interesting to speculate about some possible scenarios that may arise at the end of the high-luminosity LHC run with an integrated luminosity of 3 ab −1 . Existing estimates of the precision achievable on the Yukawa coupling to bottom quarks and the effective Higgs coupling to photons (our parameter |κ eff γγ |) suggest that, at 95% CL, the quantity λ bγ defined in (46) can be measured with a precision at least as good as 15% [12]. In Figure 9 we consider two possible future scenarios: In the first scenario the ratio Br(h → bb)/Br(h → γγ) is measured at its SM value and an upper limit of about 0.5 · 10 −6 is placed on the h → Υ(1S) γ branching ratio. Even though no observation of this rare decay is accomplished, the assumed upper bound still has a non-trivial impact, as it limits the allowed values of κ b /κ eff γγ to the right half-plane. In particular, this would exclude the possibility that κ b /κ eff γγ ≈ −1. In the second scenario the parameter λ bγ is measured close to its current value but with higher accuracy, while a rough 50% measurement of the h → Υ(1S) γ branching ratio at about (0.91 ± 0.46) · 10 −6 is obtained. In this case one could exclude the SM value κ b /κ eff γγ = 1 and limit the allowed values of this ratio to the left half-plane. These speculative results nicely indicate the power of future searches for the rare exclusive decay h → Υ(1S) γ in the high-luminosity phase of the LHC. The statistics of such a search can be approximately doubled by including also the Υ(2S) γ and Υ(3S) γ final states.
Several ideas for constraining the absolute value of the charm-quark Yukawa coupling in the high-luminosity run at the LHC or at a future 100 TeV collider have been put forward in [26]. If such a measurement can indeed be made, it implies a circular allowed region centered at (0, 0) in the plane of the parameters κ c /κ eff γγ andκ c /κ eff γγ . A measurement of the h → J/ψ γ branching ratio would be mainly sensitive to κ c /κ eff γγ and hence confine the couplings to a curved band, which intersects this region (see the center plot in Figure 6). In this way, it may be possible to perform an analysis similar to that shown in Figure 9 for the Higgs couplings to the charm quark.
If at a future 100 TeV proton-proton collider one succeeds to collect large data samples of the rare decays h → V γ → l + l − γ with V = J/ψ or Υ(nS), one might even speculate about the possibility to measure the polarization of both the vector meson and the photon by considering events in which the photon undergoes a nuclear conversion to an electron-positron pair, in analogy with what was proposed for h → γγ decay in [57] and for the B → K * γ process in [59]. With such a measurement it would be possible to differentiate between the two structures in (5), which in 3-vector notation correspond to the products ε * V ·ε * γ and ε * V ×ε * γ of polarization vectors. In this way one would become sensitive to the sign of the ratioκ q /κ q of the CP-odd and CP-even Yukawa couplings, thus breaking the degeneracy beetween the upper and lower half-planes in Figure 9.

Conclusions
We have performed a state-of-the-art analysis of the rare exclusive decays h → V γ in the context of generic extensions of the SM with modified Higgs couplings. These decays are characterized by a destructive interference between two decay topologies. The direct contribution is governed by diagrams where the Higgs boson decays into a quark anti-quark pair, from which the vector meson is formed. It is proportional to the (modified) Yukawa coupling of the Higgs to that quark flavor. Using QCD factorization and techniques developed in [20], we have derived a closed analytic expression for the direct amplitudes at next-to-leading order in α s as an infinite sum over Gegenbauer moments renormalized at the scale µ ∼ m h . In this way large logarithmic corrections are resummed, and the sensitivity of our predictions to hadronic input parameters is reduced. Power corrections to the factorization formula are suppressed by (Λ QCD /m h ) 2 or (m V /m h ) 2 and can be safely neglected. The second, indirect contribution to the decay amplitude arises from diagrams where the Higgs boson decays into a photon and an off-shell neutral gauge boson, which converts into the vector meson. Due to the smallness of the Yukawa couplings, this indirect contribution is the dominant one in all cases but h → Υ(nS) γ. In order to reduce the sensitivity to possible new-physics effects in the effective (loop-induced) hγγ vertex and thus get access to the quark Yukawa couplings, we consider the ratio of the h → V γ and h → γγ branching fractions, in which the indirect contribution drops out up to very small power corrections. In our analysis we account for the effects of flavor mixing, which can be important for the light mesons ρ 0 , ω and φ.
We have derived numerical predictions for the h → V γ branching fractions and studied the possibility of using such measurements as probes of new-physics effects on the light-quark Yukawa couplings. In the SM, the branching ratios we find typically lie in the range of few times 10 −6 , several orders of magnitude below the ATLAS bounds for the h → J/ψ γ and h → Υ(nS) γ branching ratios reported in [21]. Nevertheless, the very high yield of Higgs bosons expected in the high-luminosity phase of LHC operation, combined with a dedicated experimental effort, could make measurements of these rare processes possible. We estimate that with 3 ab −1 of integrated luminosity it will be possible to probe for O(1) modifications of the real part of the charm-quark Yukawa coupling and O(30) modifications of the real part of the strange-quark Yukawa coupling.
We have emphasized that the decays h → Υ(nS) γ provide a golden opportunity to probe for new-physics effects on the bottom-quark Yukawa couplings. Due to a fortuitous cancellation between the direct and indirect contributions to the corresponding decay amplitudes, the SM branching fractions for these modes are more than two orders of magnitude smaller than the h → J/ψ γ branching ratio. They are unobservable at the LHC and at any conceivable future collider. However, we show that with O(1) modifications of the b-quark Yukawa couplings these branching fractions can be enhanced to an observable level. Any measurement of such a decay would be a clear signal of new physics. Moreover, as we have shown, a combined measurement of the two ratios Br(h → Υ(nS) γ)/Br(h → γγ) and Br(h → bb)/Br(h → γγ) can provide complementary information on the real and imaginary parts of the b-quark Yukawa coupling. This will allow one to probe the interesting option that κ b ≈ −1 has the opposite sign as in the SM.
Several extensions of our work are possible and worth pursuing. Decays into flavor off-diagonal neutral mesons can be used to probe for flavor-violating couplings of the Higgs boson, which at tree level are forbidden in the SM. Another straightforward generalization is the application of our approach to the weak radiative decays h → M 0 Z and h → M ± W ∓ , where M can be a pseudoscalar or vector meson, as well as decays into final states containing two mesons. This is left for future work. Exclusive radiative decays of the Higgs boson are an experimentally challenging endeavor, because the expected branching fractions are very small and the final states not easy to reconstruct. Nevertheless, we have argued that studies of these processes may not be impossible in the high-luminosity phase at the LHC and, even more so, at a future 100 TeV proton-proton collider. This would present us with a unique laboratory to study the real and imaginary parts of the Yukawa couplings of bottom and charm quarks, and probe for enhanced Yukawa couplings of the lighter quarks. We cannot think of any other way in which such direct studies of Yukawa couplings could be performed.
By definition, the quantities f ⊥ ω and f ⊥ φ defined below (12) are then given by From (19) it then follows that In the limit of SU (3) flavor symmetry the ratio f ⊥ φ I /f ⊥ ω I can be replaced by 1. Note that in general the mixing angle θ ωφ and the matrix elements in (A.2) may depend on the momentum transfer k 2 . If this is the case, the values of f ⊥ ω I and f ⊥ φ I entering in (A.3) and (A.4) are different. All parameters entering the quantities δ V in (A.7) must then be evaluated at k 2 = m 2 V .

B RG evolution equations
The running quark masses and transverse decay constants obey the RG equations 3 where γ m and γ T are the anomalous dimensions of the quark mass and the QCD tensor current in (12). At two-loop order these objects were first obtained in [60] and [61], respectively. At NLO in RG-improved perturbation theory, the evolution equations have the solutions The relevant one-and two-loop coefficients of the anomalous dimensions read where C F = 4/3, C A = 3 and T F = 1/2 for SU c (3). Moreover, are the first two coefficients of the QCD β-function. Above, µ 0 is a low reference scale, at which the quark masses and decay constants are calculated using some non-perturbative approach. The evolution of the leading-twist LCDAs at NLO order has been studied in [62][63][64][65][66]. Starting at two-loop order the scale dependence of a V ⊥ n (µ) receives contributions proportional to a V ⊥ k (µ) with k = 0, . . . , n. Defining the vector a = (1, a V ⊥ 2 , a V ⊥ 4 , . . . ) containing the even Gegenbauer moments, one finds a(µ) = U (µ, µ 0 ) a(µ 0 ), where the evolution matrix U has a triangular structure with entries and U nk (µ, µ 0 ) = 0 for k > n (with even k, n ≥ 0). Explicit expressions for E n (µ, µ 0 ) and d k n (µ, µ 0 ) can be found, e.g., in Appendix A of [67]. The relevant two-loop anomalous dimensions for the Gegenbauer moments of a transversely polarized vector meson were calculated in [68,69]. Numerically, we obtain for the evolution from the low reference scale µ 0 = 1 GeV up to the high scale µ = m h where we have adjusted the number of light quark flavors when crossing the thresholds at µ = m c and m b . Due to the high value of the electroweak scale the mixing effects are strongly suppressed; for example, we obtain a V ⊥ 4 (m h ) ≈ 0.304a V ⊥ 4 (µ 0 ) − 0.007a V ⊥ 2 (µ 0 ) + 0.0008. When NLO evolution effects are included, the coefficients of the various terms shown in the second row of (29) get rescaled by factors of 1.001, 0.994, 0.984, 0.975, 0.969, respectively. Given that all present estimates of the hadronic parameters a V ⊥ n are afflicted with large theoretical uncertainties, it is sufficient for all practical purposes to use the leading-order solution (28).  Table 4: Compilation of theoretical predictions for the ratio f ⊥ V (µ)/f V at µ = 2 GeV for light vector mesons. Values marked with a dagger are obtained by taking ratios of individual results for the two decay constants. In our combinations we adopt more conservative error estimates than in some of the original references.

C Transverse vector-meson decay constants
The direct contributions to the h → V γ form factors in (20) involve the ratio f ⊥ V (µ)/f V of the vector-meson decay constants of tensor and vector currents, as defined in (8) and (12). It is reasonable to assume that non-perturbative evaluations of this ratio are subject to smaller theoretical uncertainties than calculations of the individual decay constants. For light mesons, predictions for the ratio of decay constants have been obtained using light-cone QCD sum rules and lattice QCD. Table 4 shows a compilation of available results, normalized at the reference scale µ = 2 GeV. The numbers in boldface show our own combinations of these results.
For heavy quarkonia, the ratio of decay constants can be calculated using NRQCD. Including the leading relativistic corrections [33] and one-loop QCD effects [32], we obtain where m Q is the pole mass of the heavy quark, and the dots denote higher-order terms. Numerical values for the NRQCD matrix element v 2 V for the J/ψ and Υ(nS) states with n = 1, 2, 3 have been obtained from an analysis of the electromagnetic decays V → e + e − including O(α s ) corrections and the leading relativistic effects [74,75]. The results obtained in this way, along with the adopted values of the pole masses of the heavy quarks, are compiled in Table 5. There are significant uncertainties related to the values of the heavy-quark pole masses and the NRQCD matrix elements v 2 V , which in our opinion have been underestimated 1.4 ± 0.2 0.225 + 0.106 − 0.088 [74] 0.91 ± 0.14 Υ(1S) 4.6 ± 0.1 −0.009 ± 0.003 [75] 1.09 ± 0.02 Υ(2S) 4.6 ± 0.1 0.090 ± 0.011 [75] 1.08 ± 0.02 Υ(3S) 4.6 ± 0.1 0.155 ± 0.018 [75] 1.07 ± 0.03 Table 5: NRQCD parameters for heavy quarkonia extracted from their electronic decay widths Γ(V → e + e − ), and resulting values for the ratio of decay constants.
in these analyses. It is well known that the concept of a pole mass is ill defined beyond perturbation theory and affected by renormalon ambiguities [76,77]. The values used in [74,75] are "one-loop pole masses", but the intrinsic uncertainties in these values are parametrically of order Λ QCD . Also, as emphasized in [78], the NRQCD expressions for the electromagnetic decay rates Γ(V → e + e − ) receive very large perturbative corrections at two-and three-loop order, and this prevents an accurate extraction of the non-perturbative parameters v 2 V . In view of these remarks, the errors assigned on the b-quark mass and on the NRQCD matrix elements for the Υ(nS) states seem overly optimistic. In order to be conservative, we increase these errors by a factor of 2. This yields the results shown in Table 2.

E Coefficient functions ∆ V and∆ V
Here we list the complete expressions for the CP-even coefficients ∆ V and the CP-odd coeffi-cients∆ V defined in (39). For the former ones, we obtain  (19). The corresponding expressions for the Υ(2S) and Υ(3S) mesons will not be shown explicitly, since they are very similar to that for the Υ(1S) state, see (42). In the above results the first term shows the direct contribution. The remaining terms, which originate from the power-suppressed h → γZ * → γV contribution and the effect of the off-shellness of the photon in the h → γγ * → γV contribution, are extremely small. Even for ∆ φ and assuming a SM-like Higgs couplings to strange quarks, the theoretical uncertainty in the direct contribution is an order of magnitude larger than the power-suppressed terms. Only for Υ(1S) the power-suppressed terms reach the level of 10 −2 , but still this contribution is smaller than the theoretical uncertainty in the direct contribution. Our complete expressions for the CP-odd coefficients∆ V arẽ ∆ ρ 0 − r CP = (0.068 ± 0.006) + i(0.011 ± 0.002) κ ρ 0 + (0.429 + 0.003i)κ t + (0.306 + 0.002i)κ γγ κ eff It is a good approximation to only keep the direct contributions, which are likely to give rise to the dominant effects.