Order-$v^4$ corrections to Higgs decay into $J/\psi + \gamma$

The process $H \to J/\psi + \gamma$, where $H$ is the Higgs particle, provides a way to probe the size and the sign of the Higgs-charm coupling. In order to improve the theoretical control of the decay rate, we compute order-$v^4$ corrections to the decay rate based on the nonrelativistic QCD factorization formalism. The perturbative calculation is carried out by using automated computer codes. We also resum logarithms of the ratio of the masses of the Higgs and the $J/\psi$ to all orders in the strong coupling constant $\alpha_s$ to next-to-leading logarithmic accuracy. In our numerical result for the decay rate, we improve the theoretical uncertainty, while our central value is in agreement with previous studies within errors. We also present numerical results for $H \to \Upsilon(nS) + \gamma$ for $n=1,2$, and 3, which turn out to be extremely sensitive to the Higgs-bottom coupling.


I. INTRODUCTION
The investigation of the Higgs sector of the Standard Model is one of the most important areas of particle physics today.While measuring the Higgs boson self couplings will reveal important information about electroweak symmetry breaking, the determination of the Yukawa couplings between the Higgs H and the Standard Model fermions is a direct probe of the origin of fermion masses.While current measurements of Higgs production at the LHC provide some constraint on the Higgs-top and Higgs-bottom Yukawa couplings [1], a determination of the Higgs-charm coupling is still out of reach.
The possibility of measuring the Higgs-charm coupling at the high-luminosity LHC (HL-LHC) has been studied in two different processes [2,3].One way is to measure the Higgs decay into cc, by identifying charm jets in the final state.Another way is to measure the decay of the Higgs to a charmonium and a photon [4].Compared to H → cc, the process H → charmonium+γ has an advantage that the charmonium provides a clean final state through its electromagnetic decays.Higgs decay into charmonium+γ also allows a simultaneous measurement of the size and the sign of the Higgs-charm coupling.The current upper limits for the branching ratio Br(H → J/ψ + γ) and the cross section σ(pp → ZH) × Br(H → cc) at 95% confidence level are both about 2 orders of magnitude larger than the Standard-Model predictions [2,3].
It is crucial that the decay rate Γ(H → J/ψ + γ) is in good theoretical control in order that the measurement of the rate leads to a determination of the Higgs-charm coupling.
Recently there have been many efforts to improve the theoretical prediction of the decay rate within the Standard Model [4][5][6][7][8].Especially, approaches based on nonrelativistic effective field theories allow a systematic improvement of theoretical accuracy [4,5,7,8].In the nonrelativistic QCD (NRQCD) effective field theory [9], decay and production processes involving a heavy quarkonium are given by a double series in α s and v, where v is the typical velocity of a heavy quark Q in a heavy quarkonium; for charmonium, v 2 ≈ 0.3, and for bottomonium, v 2 ≈ 0.1.Currently, the decay rates Γ(H → V +γ) for V = J/ψ or Υ(nS) for n = 1, 2, and 3 have been computed to relative order α s v 0 and v 2 accuracy [5,7,8,10].In Refs.[5,7,8], the large logarithms of m 2 H /m 2 V that appear in higher order corrections in α s , where m H is the Higgs mass and m V is the mass of the quarkonium V , have been resummed to all orders in α s by combining the NRQCD and the light-cone formalisms [11][12][13].
In this paper, we improve the accuracy of the Standard-Model prediction of the decay rates Γ(H → V + γ) for V = J/ψ or Υ(nS) for n = 1, 2, and 3 by computing the orderv 4 correction to the decay rate in the NRQCD factorization formalism.In our numerical analysis, we do not consider the ψ(2S) meson, because the nonrelativistic expansion may not converge well, and open flavor threshold effects may be important for this state.We work in the limit m 2 V /m 2 H → 0, where the calculation simplifies dramatically.In this limit, H → V + γ occurs through two distinct processes that we refer to as direct and indirect processes, see Fig. 1.In the direct process, the Higgs boson decays into a heavy quark Q and a heavy antiquark Q through the Yukawa interaction, and the Q Q pair forms a quarkonium after emitting a photon.We compute this amplitude to order-v 4 accuracy.We also resum the logarithms of m 2 H /m 2 V to all orders in α s , using the light-cone formalism, to next-to-leading logarithmic (NLL) accuracy.That is, we resum the leading and next-to-leading logarithmic corrections of the forms α n s log n (m 2 H /m 2 V ) and α n s log n−1 (m 2 H /m 2 V ), respectively, for all orders n ≥ 1 in α s .We note that the light-cone formalism applies only to the leading-order piece in the expansion in powers of m 2 V /m 2 H [11,12].In the indirect process, the Higgs boson first decays into a γ and a γ * , and the γ * evolves into a quarkonium [4].We compute the indirect amplitude to the same accuracy as the direct amplitude.We compute the direct and indirect amplitudes separately because, in the limit m 2 V /m 2 H → 0, we find simplifications in the indirect process that let us compute the indirect amplitude accurately from the known calculation of the H → γγ decay amplitude and the leptonic decay rate of the meson V .Also, the logarithms of m 2 H /m 2 V do not appear in the indirect amplitude.We note that, although the indirect amplitude involves one photon coupling more than the direct amplitude, this is compensated by the heavy-quark-Higgs Yukawa coupling in the direct amplitude.In the case of the charm quark the Yukawa coupling is y c ≈ 0.005, which indeed makes the direct amplitude numerically smaller than the indirect one.In the case of the bottom quark y b ≈ 0.018 and the two amplitudes are numerically close.See Sec.IV.
The remainder of this paper is organized as follows.In Sec.II, we compute the direct amplitude to relative order-v 4 accuracy in the NRQCD factorization formalism.We include the previously known order-α s and order-v 2 corrections and resum leading and next-toleading logarithms of m 2 H /m 2 V to all orders in α s .We compute the indirect amplitude in Sec.III.We provide our numerical results in Sec.IV, and conclude in Sec.V.

II. CALCULATION OF THE DIRECT AMPLITUDE
In this section, we compute the direct amplitude to order-v 4 accuracy in the NRQCD factorization formalism.We work at leading order in α s , but we will include the previously known order-α s v 0 correction in our final results.
We first explain the formalism that we use to compute the direct amplitude in this section.
The creation amplitude of a heavy quarkonium V with polarization vector ǫ(λ) and a photon to relative order v 4 accuracy is given by Here, m is the mass of the heavy quark, g s is the strong coupling, ψ † and χ are Pauli spinor fields that create a heavy quark and an antiquark, respectively, and E i = G i0 and B i = 1 2 ǫ ijk G kj are chromoelectric and chromomagnetic fields, respectively, where G µν is the gluon field-strength tensor.The covariant derivative D = ∇ − ig s A appears in Eq. ( 1) Operators with more than one covariant derivative are defined with The notation T ii δ ij is a shorthand for the symmetric traceless part of a tensor.The short-distance coefficients c n are perturbatively calculable quantities that do not depend on the meson state |V , while the long-distance matrix elements (LDMEs) of NRQCD operators between the vacuum |0 and the meson state |V are nonperturbative quantities.We take the meson state |V to be normalized nonrelativistically.In order to include the polarization vector in the short-distance coefficients c n in Eq. ( 1), we projected the NRQCD operators on the polarization vector of the state V |.
In Eq. ( 1), we included operators that do not contain the chromoelectric or chromomagnetic fields up to dimension 7, and operators that do contain the chromoelectric or chromomagnetic fields up to dimension 6, all of which have definite total angular momentum J = 1, charge conjugation C = −1 and parity P = −1, which are the same as V = J/ψ or Υ(nS).Throughout this paper, we denote the operators that do not contain chromoelectric or chromomagnetic fields as color-singlet operators, and the ones that do contain the chromoelectric or chromomagnetic fields as color-octet operators.Among all possible NRQCD operators, we included in Eq. (1) only the operators whose long-distance matrix elements (LDMEs) contribute to the amplitude up to relative order v 4 , based on the conservative power counting of Refs.[14][15][16][17].In this power counting, the velocity scaling of the LDME of an NRQCD operator is determined by the dimension of the operator, where a power of v is associated to a unit of dimension, and by the contribution to the meson state of the Q Q Fock state created by the operator.For J/ψ or Υ(nS), the leading Fock state contains a Q Q in the color-singlet 3 S 1 state, and this state has the scaling v −3/2 .Subleading Fock states such as the ones that contain Q Q in a color-octet state, or the ones that contain Q Q in a color-singlet D-wave state are suppressed by v and v 2 compared to the leading Fock state, respectively.The color-singlet operator ψ † σ • ǫ(λ)χ is the lowestdimensional operator that creates a Q Q in the leading Fock state ( 3 S 1 ), and so, the LDME V |ψ † σ • ǫ(λ)χ|0 scales like v 3/2 , and contributes to the amplitude at leading order in v. 4 χ scale like v 7/2 and v 11/2 , respectively, because the operators create the Q Q in the leading Fock state and have dimensions that are higher than the lowest-dimensional operator by 2 and 4, respectively.

The LDMEs of the operators ψ
Hence, these LDMEs contribute to the amplitude at relative order v 2 and v 4 , respectively.

The operator ψ
scale like v 9/2 , v 11/2 , and v 11/2 , respectively, and contribute to the amplitude at relative order v 3 , v 4 , and v 4 , respectively.For later use, we also define ratios of LDMEs as follows: There is a color-singlet operator of dimension 7 that does not appear in Eq. ( 1), which is given by 3 D 1 state, its LDME scales like v 15/2 and contributes to the amplitude at relative order v 6 .
Similarly, the color-octet operators of dimension 6 given by ψ † ǫ(λ) )χ do not appear in Eq. (1) because their LDMEs contribute to the amplitude at relative order v 5 and v 6 , respectively.The velocity scalings of these LDMEs can also be determined from the Gremm-Kapustin relations in Eqs.(D2c) and (D2d).
If we follow the power counting of Ref. [9], the color-octet LDMEs except for are suppressed beyond relative order v 4 and do not appear at the current level of accuracy.
We compute the short-distance coefficients c n appearing in Eq. ( 1) at leading order in α s by using the perturbative matching conditions obtained by replacing the meson state V with a perturbative Q Q or a Q Qg state.Since the expression in Eq. ( 1) is only valid to a limited accuracy in v, we expand the perturbative amplitude in powers of the 3-momenta of the Q, Q, and the gluon, and truncate the series to the desired accuracy.We follow a method used in Refs.[18,19], that consists in not projecting to a specific color, spin or orbital angular momentum of the Q Q state, but instead, in only requiring the Q Q or the Q Qg state to have the same J P C = 1 −− as the meson state V .This method has the advantage that fewer matching conditions are required to compute the short-distance coefficients.The caveat is that specific expressions for the matching conditions can be more complicated than when projected to specific color, spin, and orbital angular momentum states.Therefore, this method is suitable for computer-aided, automatized calculations.Also, this method can require including NRQCD operators that have same dimensions as the ones appearing in Eq. ( 1) but have LDMEs that are suppressed beyond relative order v 4 , because the matching conditions obtained in this way do not depend on the probabilities of the Q Q Fock states to be found in the meson state.Hence, in the calculation of the short-distance coefficients, we include all color-singlet operators of dimensions up to 7, and color-octet operators of dimensions up to 6, that have J P C = 1 −− .
If we replace the meson state V with a perturbative Q Q state, the amplitude occurs from order g 0 s , and the color-octet operators do not contribute to the amplitude at this order.We include all color-singlet operators up to dimension 7, which contain at most 4 covariant derivatives.Hence, we must consider the production amplitude of a Q Q state at up to 4th order in the relative momentum of the Q and the Q.We use the kinematical configuration given in Appendix A where the relative 3-momentum between the Q and the Q is given by q.The production amplitude of a Q Q state with J P C = 1 −− and a photon is given at order where we have included all color-singlet operators with J P C = 1 −− up to dimension 7. We take the states |Q and the | Q to be nonrelativistically normalized.We determine the short- , and c D 4 , along with c D 2 D (i D j) that does not appear in Eq. ( 1), by computing the left-and right-hand sides in perturbative QCD and perturbative NRQCD, respectively, and comparing the two sides order by order in the expansion in powers of q up to 4th order.
In order to compute the remaining short-distance coefficients corresponding to the coloroctet LDMEs, we consider the production amplitude of a Q Qg state with J P C = 1 −− which occurs from order g s .We use the kinematical configuration given in Appendix A where the relative 3-momentum between the Q and the Q is given by q 1 and the relative 3-momentum between the Q Q pair and the gluon is given by q 2 .At order g s , the color-octet LDMEs that appear in Eq. ( 1) have matrix elements that are either linear or quadratic in the momenta q 1 or q 2 when the meson state V is replaced by the Q Qg state.We must also include all color-octet operators of dimensions up to 6 with J P C = 1 −− that do not appear in Eq. ( 1), whose matrix elements can also be either linear or quadratic in the momenta q 1 or q 2 .
The color-singlet operators in Eq. ( 1) can also contribute to the Q Qg amplitude at order g s through the gauge fields in the covariant derivatives and through insertions of NRQCD vertices.An NRQCD vertex insertion at order g s involves a heavy-quark propagator, which can produce a factor of 1/|q 2 |.Hence, it is necessary to include all color-singlet operators that contain at most 3 covariant derivatives.Since the lowest dimensional color-singlet operator we consider is of dimension 3, and the highest dimensional color-octet operators are of dimension 6, we need to include NRQCD vertices up to 1/m 3 accuracy.That is, we need to consider two-fermion operators of dimensions up to 7 in the NRQCD Lagrangian, which are given by [20] where c.c. stands for the charge conjugated contribution of the preceding terms.Since we only consider the matching at tree level, we only include the Wilson coefficients at order α 0 s in Eq. ( 5).The production amplitude of a Q Qg state with J P C = 1 −− and a photon at order g s is given by where the left-hand side is calculated in perturbative QCD and is expanded in powers of q 1 and q 2 up to quadratic accuracy.We again take the states |Q and | Q to be nonrelativistically normalized.We determine the short-distance coefficients c B , c DE 0 , and c DE 1 , along with c DE ′ 1 and c DE 2 that do not appear in Eq. ( 1) by computing the left-and right-hand sides in perturbative QCD and NRQCD, respectively, and comparing the two sides order by order in the expansion in powers of q 1 and q 2 up to quadratic order 1 .
In the following sections, we compute the short-distance coefficients c n explicitly.We first calculate the c n in fixed-order perturbation theory, where the QCD amplitudes on the lefthand sides of Eqs. ( 4) and ( 6) are computed at leading order in α s .We obtain corrections 1 The operator ψ † ǫ(λ) χ does not appear in Ref. [16].If we compute the NRQCD matrix elements on the right-hand side of Eq. ( 6) explicitly, the matrix element of this operator is the only matrix element that is quadratic in q 2 .Hence, if we ignore the contribution that is proportional to |q 2 | 2 in both sides of Eq. ( 6), we can ignore this operator from the matching condition without affecting the calculation of the short-distance coefficients in Eq. ( 1).
to the direct amplitude of relative order v 4 which is new in this work, and reproduce the known order-v 2 correction in the fixed-order calculation.We then compute the c n in the light-cone approach, which is valid at leading order in m 2 /m 2 H , that allows us to resum logarithms of m 2 H /m 2 to all orders in α s .We obtain new corrections of relative order v 4 in the light-cone approach, and reproduce the previously calculated order-v 2 correction.We include the order-α s correction to the direct amplitude using the light-cone approach.

A. Fixed-order calculation
At order g 0 s , the direct amplitude for where p γ and ǫ * γ are the momentum and the polarization vector for the photon in the final state.We use the physical gauge for the photon polarization vector, so that ǫ * γ • p γ = 0. Here, e = √ 4πα is the electric charge, e Q is the fractional charge of the heavy quark Q, and is the Yukawa coupling of the Higgs and Q, with G F the Fermi constant.The momenta of the Q and Q are given by p 1 and p 2 , respectively, so that the momentum of the H is where q = 1 2 (p 1 − p 2 ).We choose the heavy-quark mass appearing in y Q to be m(µ), which is the MS mass of the heavy quark Q at scale µ; as we will see in the next section, this choice simplifies the logarithms that appear in the order-α s correction.Since C, P , and T are conserved in the amplitude in Eq. ( 7), the Q Q can only be created with C = −1.We first express the Dirac bilinears ū(p 1 )γ µ v(p 2 ) and ū(p 1 )γ µ γ ν v(p 2 ) in terms of the 3-momenta of the Q and Q in the Q Q rest frame.This can be accomplished by using the method described in Appendix B. Then, we expand the amplitude in powers of q, keeping terms up to relative order q 4 /m 4 .The resulting expression for the amplitude is then a linear combination of the Cartesian tensors built from ξ † ση and q of the form ξ † σ i ηq j • • • q k up to rank 5 and of the form ξ † ηq i q j • • • q k up to rank 4. The contribution from these Cartesian tensors to the total angular momentum J = 1 can be obtained by a reduction method developed in Ref. [21].Finally, the P = −1 contribution is obtained by keeping only the contribution odd in parity, where the parity transform of the Q Q amplitude is given by the replacements q → −q, ξ † ση → −ξ † ση, and ξ † η → −ξ † η.We use the Mathematica package FeynCalc [22,23] and the FeynOnium [24] package to automatize the calculation of the amplitude and the consequent reduction to the J P C = 1 −− contribution.By comparing the J P C = 1 −− contribution of the Q Q amplitude with the right hand side of Eq. ( 4), we obtain the short-distance coefficients where we define r ≡ 4m 2 m 2

H
. The short-distance coefficients c 0 and c D 2 agree with Refs.[4,5], except that our results differ by an overall sign that originates from the sign convention of the J = 1 state employed in Refs.[4,5].We also reproduce the short-distance coefficient c D 4 that can be obtained from Ref. [5].The results for c D (i D j) and c D 2 D (i D j) are new.
The remaining short-distance coefficients corresponding to the color-octet LDMEs are computed from the direct amplitude for H → Q Qg + γ at order g s , which is given by where ǫ * g is the polarization vector of the gluon.The total momentum P of the Q Qg state is given by P = p 1 + p 2 + k g , and the momentum of the H is given by P H = P + p γ , so that in the rest frame of the Q Qg, where , and q 2 = 1 6 (2k g −p 1 −p 2 ).Due to C conservation at this order, the Q Qg can only be produced in a colorsinglet C = −1 state.We use again the Mathematica package FeynCalc and the FeynOnium package to automatize the calculation of the Q Qg amplitude.After expressing the amplitude in terms of the 3-momenta of the Q, Q and g, we expand the amplitude in powers of q 1 and q 2 up to relative order q 2 1 /m 2 , q 2 2 /m 2 , and The resulting expression for the amplitude is then a linear combination of the Cartesian tensors built from ξ † ση, ǫ * g , q 1 , and q 2 up to rank 4. In order to keep only the contribution with negative parity, we keep only the contribution odd in parity, where the parity transform of the Q Qg amplitude is given by the replacements q 1 → −q 1 , q 2 → −q 2 , ǫ * g → −ǫ * g , ξ † ση → −ξ † ση, and ξ † η → −ξ † η.After reducing the Cartesian tensors of odd rank to total angular momentum J = 1, we compare the J P C = 1 −− contribution with the amplitude in the right-hand side of Eq. ( 6) to obtain the short-distance coefficients where The short-distance coefficients in Eqs. ( 9) and ( 12) allow us to compute the direct amplitude to relative order v 4 accuracy.Because the indirect amplitude will be available only at leading order in m 2 /m 2 H (see Sec. III), only the limit m 2 /m 2 H → 0 of the short-distance coefficients in Eqs. ( 9) and ( 12) will be employed in the total decay amplitude.This amounts to setting r = 0.
We note that the method described in this section can be easily applied to compute contributions of the QCD amplitude with different J P C .For example, the J P C = 1 +− contribution can be obtained by keeping parity even terms in the amplitude.The J P C = 1 +− contribution can then be used to obtain the short-distance coefficients for the H → h c + γ amplitude.We show this result in Appendix C.

B. Light-cone calculation
In the fixed-order calculation, the short-distance coefficients contain contributions from the scales m and m H . Since m H is much larger than m, the logarithms of m 2 H /m 2 that appear in corrections of higher orders in α s can potentially spoil the convergence of the perturbation series.If we work at leading power in m 2 /m 2 H , the light-cone approach provides a factorization formula that separates the contribution at the scale m from the contribution at the scale m H [11,12].The light-cone approach also enables us to resum the logarithms of m 2 H /m 2 to all orders in α s by solving an evolution equation.The factorization formula is given by [11,12,25] where T H (x, µ) is the perturbative hard part that contains the contribution at the scale m H , while the contribution at scales m and below are contained in the decay constant f ⊥ V (µ) and the light-cone distribution amplitude (LCDA) φ ⊥ V (x, µ) of the meson V .The decay constant and the LCDA are nonperturbative quantities defined by The nonlocal operator Q α (x) is defined by where P is the momentum of the meson V , and the decay constant f ⊥ V (µ) is defined by integrating Eq. ( 14) over x and considering the normalization of the LCDA, which is given by Here, Q(x) is the QCD quark field.The light-cone vectors n and n are given by n = m H pγ•P H p γ and n = 2 m H P H − n, which satisfy n 2 = n2 = 0 and n • n = 2. Here, P H and p γ are the momenta of the H and the photon, respectively.For any 4-vector a µ , we , where P is the path ordering operator, ensures the gauge invariance of the nonlocal operator.
The hard part T H (x, µ) has been computed to next-to-leading order (NLO) accuracy in α s and is given by [6,25] where , and N c = 3 is the number of colors.The imaginary part in Eq. ( 16) comes from the order-α s correction where the virtual lines can be on shell.The expression for the NLO correction to T H (x, µ) in Eq. ( 16) is valid when the decay constant f ⊥ V (µ) and the LCDA φ ⊥ V (x, µ) are renormalized in the MS scheme, and the heavy-quark mass in the Yukawa coupling y Q is the MS mass at the scale µ.If we use the heavy-quark pole mass m instead of the MS mass at the scale µ, the order-α s correction in T H (x, µ) involves a logarithm of m 2 H /m 2 ; instead, using the MS mass in the Yukawa coupling y Q ensures that T H depends only on the scale m H and removes this logarithmic contribution from it [6].
The nonlocal operator Q α (x) has an anomalous dimension known to NLO in α s that allows us to resum logarithms of m 2 H /m 2 to NLL accuracy [11,[26][27][28][29].In order to resum the logarithms of m 2 H /m 2 that appear in corrections of higher orders in α s to the short-distance coefficients, we apply the factorization formula [Eq.(13)] to the perturbative amplitudes This involves calculating the decay constant and the LCDA with the meson state V replaced by the perturbative Q Q and Q Qg states.Then, each of the short-distance coefficients c n is given by a convolution of the hard part T H and a distribution in x that satisfies the same evolution equation as the nonlocal operator in Eq. ( 14).Equivalently, we can apply the NRQCD factorization formula to Eq. ( 14) to obtain an expression of the decay constant and the LCDA in terms of NRQCD LDMEs, so that [13] where the short-distance coefficients cn (x) are computed from the matching conditions that are similar to Eqs. ( 4) and (6), where the c n on the right-hand sides are replaced by cn (x), and the left-hand sides are replaced by , respectively.In the matching conditions, the scale µ 0 in Eq. ( 17) is the scale where the NRQCD LDMEs on the right-hand side are defined.The factorization formula [Eq.( 13)] implies that It is worth noting that Q Q|Q α (x)|0 and Q Qg|Q α (x)|0 contain contributions from both negative and positive charge conjugation.From the fact that the charge conjugated of the operator Q α (x) is given by −Q α (1 − x), we can see that the contributions of negative charge conjugation to Q Q|Q α (x)|0 and Q Qg|Q α (x)|0 are given by the contribution symmetric If we replace the meson state V with the Q Q state, we obtain We use the same strategy as the fixed-order calculation to obtain the contribution with J P C = 1 −− ; we express the Dirac bilinears in terms of the 3-momenta of the Q and the Q in the Q Q rest frame, and then we expand Eq. ( 19) in powers of q up to relative order q 4 /m 4 .Then, Eq. ( 19) is given by a linear combination of the delta function δ(x − 1/2) and its derivatives.In order to keep only the contribution with C = −1, we ignore the odd derivatives of δ(x − 1/2).The J = 1 contribution is then obtained by reducing the Cartesian tensors of the form ξ † ηq i q j • • • q k up to rank 4 and of the form ξ † σ i ηq j • • • q k up to rank 5 using the reduction method of Ref. [21].We keep only the contribution with negative parity in order to obtain the J P C = 1 −− contribution.From the matching condition (17) we obtain The short-distance coefficients c0 (x) and cD 2 (x) agree with known results in Ref. [5].The agreement with the fixed-order calculation can be easily verified using Eq. ( 18).
To compute the remaining short-distance coefficients corresponding to the color-octet LDMEs, we replace the meson state V with the Q Qg state to obtain After expressing the amplitude in terms of the 3-momenta of the Q, Q and g, we expand the amplitude in powers of q 1 and q 2 up to relative order q 2 1 /m 2 , q 2 2 /m 2 , and |q 1 ||q 2 |/m 2 .The resulting expression for the amplitude is then a linear combination of the Cartesian tensors built from ξ † ση, ǫ * g , q 1 , and q 2 up to rank 4. We obtain the matching condition for Q Qg with J P C = 1 −− by keeping only the contributions symmetric in x ↔ 1 − x, i.e., C = −1, reducing the Cartesian tensors to J = 1, and keeping only the negative parity contributions.
The resulting matching condition leads to the following short-distance coefficients The agreement with the fixed-order calculation can be easily verified using Eq.(18).
By integrating the short-distance coefficients cn (x) over x, we obtain an expression for f ⊥ V valid up to relative order v 4 .If we include the correction of relative order α s v 0 in the MS scheme from Ref. [25], we obtain The corrections at relative order v 2 agree with Ref. [5].The logarithm of µ 2 0 /m 2 in Eq. ( 23) is the remnant of the renormalization of the decay constant in the MS scheme.
Since the short-distance coefficients cn (x) at leading order in α s are linear combinations of δ(x − 1/2), δ (2) (x − 1/2) and δ (4) (x − 1/2), the LCDA can be written as where and φ ⊥ V (αs) (x, µ 0 ) is the order-α s v 0 correction in the MS scheme given by [25] φ where the plus and plus-plus distributions are defined by The logarithm of µ 2 0 /m 2 in Eq. ( 26) is the remnant of the renormalization of the LCDA in the MS scheme.Note that The decay constant and the LCDA computed in Eqs.(23,24) allow us to resum logarithms of m2 H /m 2 that appear in the QCD corrections to the direct amplitude to all orders in α s .This is accomplished by solving the evolution equation where the evolution kernel V T (x, y) is currently known to NLO accuracy in α s [11,[26][27][28][29].We obtain x, µ 0 ) at scale µ 0 by solving this evolution equation.The direct amplitude with logarithms of m 2 H /m 2 resummed to all orders in α s is then given by Eq. ( 13), by setting µ ∼ m H and µ 0 ∼ m, so that T H (x, µ) is free of The accuracy of the resummation is limited by the accuracy of the evolution kernel; since the evolution kernel is known up to NLO accuracy, we can resum the logarithms to NLL accuracy.The resummation can be carried out using the Gegenbauer polynomials C (3/2) n (2x − 1), which are the eigenfunctions of the LO evolution kernel [30].
The convolution in Eq. ( 13) is given by where TH (n, µ) and φ⊥ V (n, µ) are Gegenbauer moments defined by The solution of the evolution equation in terms of Gegenbauer moments leads to [26][27][28][29] φ⊥ Explicit expressions of U f ⊥ V (µ, µ 0 ) and U n 1 n 2 (µ, µ 0 ) can be found in Ref. [7].We note that U f V (µ, µ 0 ) and U n 1 n 2 (µ, µ 0 ) depend only on µ, µ 0 and the evolution kernel, and are inde- x, µ 0 ) contains singular distributions such as the delta function and their derivatives, the sum in Eq. (32a) can converge badly.The authors of Refs.[7,8] developed a method to overcome the problem of nonconvergence by defining the sum over n 2 in Eq. (32a) as an Abel sum, so that The numerical value of the Abel sum is then evaluated by using Padé approximants of the truncated series.We employ this Abel-Padé method [7,8] to compute the convolution in Eq. (32a).
It is more convenient to express the direct amplitude with resummed logarithms in terms of the decay constant and the LCDA convolved with the hard part, i.e.Eq. ( 13), than to resum the logarithms for each of the short-distance coefficients in Eqs. ( 9) and ( 12).We note that the resummation of logarithms can be carried out separately for the decay constant and the individual contributions to the LCDA in Eq. ( 24), so that where

III. CALCULATION OF THE INDIRECT AMPLITUDE
The indirect amplitude proceeds from H → γ * γ, followed by γ * → V .Since we work in the limit m 2 V /m 2 H → 0, the partial amplitude H → γ * γ can be replaced by the decay amplitude of the Higgs to two photons.Then, the indirect amplitude is given by [4] where the factor 16πm H compensates for the factors that are necessary to obtain the decay rate Γ(H → γγ) from the squared amplitude.We take the scale of the QED coupling at the vertices associated with the virtual photon with virtuality m V to be µ 0 , and the scale of the vertex for the real photon to be 0. The factor α(µ 0 )/α(0) replaces one QED coupling constant at scale 0 in the H → γγ amplitude with the QED coupling constant at scale µ 0 .
We ignore the small imaginary part in the H → γγ amplitude [6,7].The decay constant f V is defined by Since f V is a conserved current in QCD, it does not undergo renormalization.Hence, the indirect amplitude is free of logarithms of m 2 H /m 2 in the limit m 2 V /m 2 H → 0 if one ignores the higher-order electroweak corrections to the indirect amplitude.We can express f V in terms of NRQCD LDMEs up to relative order v 4 using the same techniques we employed to compute the direct amplitude and the decay constant f ⊥ V .We obtain where we included the order-α s v 0 correction computed in Refs.[31,32].We note that f V is positive at leading order in α s and v.If we assume f V > 0, an accurate numerical value for f V can be obtained from the leptonic decay rate so that .
Note that iM ind and iM dir have opposite signs, so that when calculating the decay rate, the direct and indirect amplitudes interfere destructively [4].

IV. NUMERICAL RESULTS
We now present our numerical results for the Higgs decay rate into V +γ for V = J/ψ and Υ(nS) for n = 1, 2, and 3 based on our calculation of the direct amplitude to relative order v 4 accuracy.Our expressions for the direct and indirect amplitudes are given in Eqs.(34) and (40), respectively.If we write where Φ is the phase space and normalization factor given by We present our numerical results for A dir and A ind in the following sections.Then, using the resulting values of A dir and A ind , we compute the decay rate Γ(H → V + γ) from the formula in Eq. (41).

A. Indirect amplitude
We compute the numerical values for the indirect amplitude using Eq. ( 40).We compute the decay constant f V from the measured leptonic decay rate: where we take µ 0 = m V .Then We note that the resulting expression for A ind does not depend on α(µ 0 ).
We use the following input parameters to compute A ind numerically.We take the PDG value for the Higgs mass m H = 125.18± 0.16 [33], and the numerical value for the Higgs two-photon decay rate to be Γ(H → γγ) = 9.34 × 10 −6 GeV, which is computed from the tabulated results of the two-photon branching ratio and the total decay rate of the Higgs in Refs. [34,35].We also take the PDG values for the meson masses m V and use the measured partial widths Γ(V → e + e − ) for the leptonic decay rates [33].We list the values for m V and Γ(V → ℓ + ℓ − ) that we use to compute the indirect amplitude in Table .I. The QED coupling constant at scale 0 is taken to be α(0) = 1/137.036.
We consider the following sources of uncertainties in A ind .The uncertainty in the decay rate Γ(H → γγ) is taken to be 0.01 times the central value, as estimated in Ref. [34] from the higher-order corrections to the decay rate.We consider the experimental uncertainties in Γ(V → ℓ + ℓ − ).We estimate the uncalculated correction of relative order m 2 V /m 2 H to be m 2 V /m 2 H of the central value.We ignore the negligibly small uncertainties in m H and m V compared to other sources of uncertainties.We add the uncertainties in quadrature.Our numerical results for A ind are shown in Table .II.We note that the uncertainties in A ind are less than 2% of the central values.
We can compare our results for A ind with a previous calculation in Ref. [8].Our calculation is equivalent to the one in Ref. [8], except that we use an value of the measured Higgs mass from Ref. [33], which has smaller uncertainties than what was employed in Ref. [8].Our results for A ind in Table .II are compatible with those in Ref. [8] within uncertainties.

B. Direct amplitude
We compute the numerical values of the direct amplitude using Eq. ( 34), so that We now discuss our strategy to compute Eq. ( 45) numerically.The decay constant and DE 1 V , see Eq. ( 23).The dependence on the leading-order LDME V |ψ † σ • ǫ(λ)χ|0 can be eliminated by rewriting the decay constant as and by obtaining the numerical value for f V from the measured leptonic decay rate using Eq. ( 43).The DE 0 V , and DE 1 V , see Eq. (25b), while the remaining contributions in Eq. ( 24) depend only on v 4 S V .We can eliminate the ratios B V and DE 0 V in f ⊥ V (µ 0 ) and φ ⊥ V (2) (x, µ) by using the Gremm-Kapustin relations in Eqs.(D2).We obtain and When computing the convolution in Eq. ( 45), the terms and 4 can contain cross terms that go beyond our current level of accuracy.In order to avoid such contributions, we ignore the cross terms that contribute to the direct amplitude beyond relative order α s v 0 and v 4 , so that We use the results in Eqs.(49) to compute the direct amplitude.Similarly, when calculating the convolution in Eq. ( 45), we ignore the order-α s correction to T H (x, µ) for n = α s , 2, and 4.
We first discuss the numerical input parameters necessary for the computation of the direct amplitude.We compute the decay constant f V from Eq. ( 43) using the measured leptonic decay rate, and α(µ 0 ) = 1/132 regardless of the vector meson state V .We set the central values of µ 0 and µ to be m V and m H , respectively.We compute m(µ), α s (µ 0 ) and α s (µ) using RunDec [36][37][38].We take the QED coupling constant at scale µ to be α(µ) = 1/128.We resum the logarithms in µ/µ 0 to NLL accuracy using the Abel-Padé method [7,8].We take the number of active quark flavors in the evolution kernel to be n f = 4 and 5 for scales below and above m b , respectively.The ratios m 2 v 2 S V for V = J/ψ and Υ(nS) have been obtained from potential-model (Cornell potential) calculations in Refs.[39,40]: The first uncertainty comes from the variation of the potential-model parameters, and the second uncertainty is taken to be ±0.3 and ±0.1 times the central value for V = J/ψ and V = Υ(nS) respectively, which comes from the uncalculated corrections of relative order v 2 .
The potential-model calculations in Refs.[39,40] also let us compute the binding energies E V = m V − 2m, which are given by E J/ψ = 0.306 +0.039 −0.041 ± 0.092 GeV, (51a) E Υ(2S) = 0.482 ± 0.032 ± 0.048 GeV, (51c) where the uncertainties are as in Eqs.(50).The uncertainty in E V from the variation of the potential-model parameters are correlated with the uncertainty in the ratio m 2 v 2 S V .We use these values of the binding energies to compute the heavy-quark mass through the The authors of Refs.[39,41] also found the relation , which is valid if the LDMEs are computed in a potential model like the Cornell potential and in dimensional regularization.By using this relation, we take the central value for the ratio v 4 S V to be ( v 2 S V ) 2 and take the uncertainty to be ±0.3 and ±0.1 times the central value of v 4 S V for V = J/ψ and V = Υ(nS), respectively.The ratios v 2 D V and DE 1 V are not known; since these ratios scale as v 4 , we take the central values of the ratios to be 0 and take the uncertainties to be ±0.09for V = J/ψ, and ±0.01 for V = Υ(nS), respectively.
We now list the sources of uncertainties in A dir .We account for the uncertainties in the NRQCD LDMEs as discussed above.We vary the scales µ and µ 0 between 1 2 m H < µ < 2m H and 1 2 m V < µ 0 < 2m V , while we ignore the negligibly small shifts in the QED couplings α(µ 0 ) and α(µ) from scale variations.We also ignore the uncertainties from m H and m V .
We consider the uncertainty in f V that originates from the experimental uncertainties in the leptonic decay rate.We add the uncertainties in quadrature.Our numerical results for A dir are shown in Table .II.The imaginary part of A dir comes from the imaginary part in the order-α s correction to T H (x, µ).Note that the uncertainties in the real and imaginary parts of A dir are correlated.
For J/ψ, the uncertainties from the LDMEs are comparable to the uncertainties from the variations of µ 0 and µ.If we ignore the uncertainties from scale variations, the uncertainty in Re[A dir ] is about 8% of the central value, and the uncertainty in Im[A dir ] is about 10% of the central value.This is comparable to the nominal size of the relative order-v 4 correction.
For Υ(nS), the uncertainties are dominated by scale variations.If we ignore the uncertainties from scale variations, the uncertainties in the real and imaginary parts of A dir are about 1% of the central values, which are comparable to the nominal size of the relative order-v 4 correction.
While the uncertainties from the LDMEs are expected to be reduced when the LDMEs v 2 D V and DE 1 V are constrained, the reduction of the uncertainties from variations of scales would require calculation of the order-α 2 s and order-α s v 2 correction to the decay constant and the LCDA, and the order-α 2 s correction to T H (x, µ).We again compare our results for A dir with the calculation in Ref. [8].Our results for A dir in Table .II are compatible with those in Ref. [8] within uncertainties.For J/ψ, the uncertainty for A dir is smaller than in Ref. [8], owing to the explicit calculation of the relative order-v 4 corrections included in this work.On the other hand, for Υ(nS), the uncertainties for A dir are larger than those in Ref. [8]; the main reason is that in Ref. [8], the uncertainty from variations in the scales µ 0 and µ were not taken into account, and instead, the uncertainties from the uncalculated order-α 2 s and order-α s v 2 corrections to the real part of A dir were estimated to be C   value.On the other hand, for V = Υ(nS), A ind and A dir are comparable.Due to the large cancellation between A ind and A dir for Υ(nS), the uncertainty in Γ[H → Υ(nS) + γ] is sensitive to the uncertainty in A dir .
Our results are compatible with the previous calculation in Ref. [8]  We note that the decay rates Γ(H ]χ|0 that are currently unknown. In our numerical results, we have estimated their sizes according to the conservative power counting in Refs.[14][15][16][17], and included their effects in the uncertainties.For V = J/ψ, the uncertainties from these unknown LDMEs are significant compared to the total uncertainty in Γ(H → J/ψ + γ).Therefore, knowledge of these LDMEs will improve the accuracy in the prediction of the decay rate Γ(H → J/ψ + γ).Also the systematic inclusion of relative order v 2 effects in the determination of the matrix element significantly improve the determination of the H → J/ψ + γ decay rate.A calculation of the LDMEs in potential NRQCD may help constraining these LDMEs [15].
In Ref. [4], the authors estimated that the process H → J/ψ + γ could be measured at the HL-LHC experiment through the leptonic decays of J/ψ into e + e − and µ + µ − .However, a more recent study by the ATLAS Collaboration found that the HL-LHC would only be able to put an upper bound for the decay rate Γ(H → J/ψ + γ) at 95% confidence level that is about 15 times larger the Standard Model value [42].Although the prospect for measuring the process H → J/ψ + γ at the HL-LHC does not look so good at the moment, it is possible that the experimental methods will improve over time; it is also possible where we define (see also Eq. ( 2)) The color-octet matrix elements for the J P C = 1 +− case are obtained from the color-octet matrix elements for the J P C = 1 −− case in Eq. ( 6) by making the replacements g s E → g s B and g s B → g s E.
From the J P C = 1 +− contributions to the Q Q and Q Qg amplitudes in Eqs. ( 7) and (10), we obtain the short-distance coefficients where r = 4m 2 m 2
Recently, a computation of the decay rate Γ(H → h c + γ) in the NRQCD factorization formalism at leading order in v has appeard in Ref. [44].This is equivalent to our calculation of the short-distance coefficient c 1 in Eq. (C4a), which leads to the expression for the decay rate Γ(H → h c + γ) at leading order in v that is given by where the sum is over the polarizaton of the h c and Φ is the phase space and normalizaiton factor given in Eq. (42).Our result in Eq. (C5) agrees with the decay rate computed in Ref. [44].
Appendix D: Gremm-Kapustin relations The Gremm-Kapustin relations [45] are obtained from where O is an NRQCD operator, and H is the NRQCD Hamiltonian.Computing the commutator [O, H] leads to the following relations where, in calculating the commutator [O, H] we included operators in the Hamiltonian up to 1/m 2 accuracy, with the Wilson coefficients at order α 0 s , and kept NRQCD operators up to dimension 7. Hence, the relations in Eqs.(D2) are valid up to order v 4 relative to the leading-order LDME V |ψ † σ • ǫ(λ)χ|0 and at leading order in α s .These relations can also be verified in perturbation theory.
The Gremm-Kapustin relations provide a way to identify the velocity scalings of the LDMEs that are suppressed beyond the conservative power counting of Refs.[14][15][16][17].Since the binding energy m V − 2m scales like mv 2 , the left-hand side of Eq. (D2c) is suppressed by v 6 compared to the leading-order LDME V |ψ † σ • ǫ(λ)χ|0 .Therefore, the LDME V |ψ † ǫ i (λ)σ j ( ← → D (i g s E j) + g s E (i ← → D j) )χ|0 does not contribute to the amplitude in Eq. ( 1) at relative order v 4 accuracy because it is suppressed by at least v 6 compared to the leadingorder LDME, and scales like v 15/2 .Similarly, the left-hand side of Eq. ( D2d) is suppressed by v 5 compared to the leading-order LDME, and hence, the LDME on the right-hand side of Eq. (D2d) scales like v 13/2 and does not contribute to the amplitude in Eq. ( 1) at relative order v 4 accuracy.

FIG. 1 :
FIG. 1: Feynman diagrams for the (a) direct amplitude at order α 0 s and (b) the indirect amplitude for the process H → V + γ.
the positive charge conjugation contribution of the LCDA, which is antisymmetric in x ↔ 1 − x, does not contribute to the decay amplitude, consistently with the conservation of charge conjugation in the amplitude.Hence, in order to keep only contributions of negative charge conjugation in Q Q|Q α (x)|0 and Q Qg|Q α (x)|0 , we just need to keep contributions that are symmetric in x ↔ 1 − x.
that future experiments at the International Linear Collider, the Circular Electron Positron Collider, the Compact Linear Collider, or the Future Circular Collider will be able to probe such processes.In the case of the Υ(nS), the Standard-Model values for the decay rates Γ[H → Υ(nS) + γ] are about three orders of magnitude smaller than Γ(H → J/ψ + γ), owing to large cancellations between the direct and indirect amplitudes.We show the decay rates Γ(H → J/ψ + γ) and Γ[H → Υ(nS) + γ] when the Higgs-charm coupling is rescaled by a factor κ c and the Higgs-bottom coupling is rescaled by a factor κ b compared to the Standard Model in Fig.2.Due to the cancellation between direct and indirect amplitudes that occurs when κ b ≈ 1, the decay rates Γ[H → Υ(nS) + γ] are highly sensitive to κ b , so much so that the combined rate of Γ[H → Υ(nS) + γ] for n = 1, 2, and 3 increases to be larger than one half of the Standard-Model value of Γ(H → J/ψ + γ) for κ b −1 or κ b 3. Therefore, the numerical results presented in this paper may be useful in determining the size and sign of the Higgs-charm and even more the Higgs-bottom couplings when these processes are measured in future experiments.and iM[H → Q Qg(J P C = 1 +− ) + γ] χ|0 scales like v 11/2 and contributes to the amplitude at relative order v 4 .The color-octet operators in Eq. (1) create Q Q in color-octet states where either the orbital or the spin angular momentum is different from that of the leading Fock state by 1.The contributions of such Fock states are suppressed by v compared to the leading Fock state.Hence, the LDMEs

TABLE I :
[33]es for the meson masses m V and the leptonic decay rates Γ(V → ℓ + ℓ − ) used in the numerical calculation of the direct and indirect amplitudes.All values are taken from Ref.[33].
[8]ios Br(H → V + γ) are shown in TableIII.The corrections computed in this work improve the theoretical accuracy in the prediction of the decay rate Γ(H → J/ψ + γ) compared to a previous calculation in Ref.[8].If we would have used the same method to estimate the uncertainties as in Ref.[8], we would have found uncertainties in Γ[H → Υ(nS) + γ] reduced by almost a factor of two compared to our results in TableIII, so that they would become comparable to the results in Ref.[8], as the error of Γ[H → Υ(nS) + γ] would be dominated by the scale uncertainties.However, due to the difference in the method to estimate uncertainties, our results for Γ[H → Υ(nS) + γ] have larger uncertainties than those of Ref.[8].
[8]vious section, if we would use the same estimates for the uncertainties used in Ref.[8], then the uncertainties in A dir would be reduced, leading to uncertainties in Γ[H → Υ(nS)+γ] and Br[H → Υ(nS) + γ] that are smaller than those of Ref.[8].branching