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

TL;DR

This paper calculates higher-order relativistic corrections to the Higgs decay into a J/psi meson and a photon, improving theoretical precision and exploring sensitivity to Higgs-charm and bottom couplings.

Contribution

It introduces order v^4 corrections and resums logarithms for the decay rate, enhancing the theoretical accuracy of Higgs decay predictions.

Findings

Decay rate predictions are more precise with reduced uncertainties.

Results are consistent with previous studies within error margins.

Decay channels involving Upsilon are highly sensitive to bottom quark coupling.

Abstract

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 boson 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.

Tables

Table 1. Table 1: Values for the meson masses m V subscript 𝑚 𝑉 m_{V} and the leptonic decay rates Γ ​ ( V → ℓ + ​ ℓ − ) Γ → 𝑉 superscript ℓ superscript ℓ \Gamma(V\to\ell^{+}\ell^{-}) used in the numerical calculation of the direct and indirect amplitudes. All values are taken from Ref.  Tanabashi:2018oca .

$V$	$m_V$ (GeV)	$\Gamma \left(V\to {\ell }^{+}{\ell }^{-}\right)$ (keV)
$J/\psi$	$3.0969$	$5.55\pm 0.14\pm 0.02$
${\rm Y}\left(1S\right)$	$9.4603$	$1.340\pm 0.018$
${\rm Y}\left(2S\right)$	$10.02326$	$0.612\pm 0.011$
${\rm Y}\left(3S\right)$	$10.3352$	$0.443\pm 0.008$

Table 2. Table 2: Numerical results for the amplitudes 𝒜 ind subscript 𝒜 ind {\cal A}_{\rm ind} and 𝒜 dir subscript 𝒜 dir {\cal A}_{\rm dir} .

$V$	${𝒜}_{\mathrm{ind}}\times {10}^5$ (GeV1/2)	${𝒜}_{\mathrm{dir}}\times {10}^5$ (GeV1/2)
$J/\psi$	$-\left({11.73}_{-0.16}^{+0.16}\right)$	$\left({0.631}_{-0.080}^{+0.071}\right)+\left({0.065}_{-0.012}^{+0.015}\right)i$
${\rm Y}\left(1S\right)$	$\left({3.288}_{-0.033}^{+0.033}\right)$	$-\left({2.719}_{-0.142}^{+0.136}\right)-\left({0.291}_{-0.040}^{+0.055}\right)i$
${\rm Y}\left(2S\right)$	$\left({2.158}_{-0.026}^{+0.026}\right)$	$-\left({1.896}_{-0.104}^{+0.101}\right)-\left({0.197}_{-0.027}^{+0.037}\right)i$
${\rm Y}\left(3S\right)$	$\left({1.808}_{-0.022}^{+0.022}\right)$	$-\left({1.614}_{-0.093}^{+0.090}\right)-\left({0.164}_{-0.023}^{+0.031}\right)i$

Table 3. Table 3: Numerical results for Γ ​ ( H → V + γ ) Γ → 𝐻 𝑉 𝛾 \Gamma(H\to V+\gamma) and Br ​ ( H → V + γ ) Br → 𝐻 𝑉 𝛾 {\rm Br}(H\to V+\gamma) .

$V$	$\Gamma \left(H\to V+\gamma \right)$ (GeV)	$\mathrm{Br}\left(H\to V+\gamma \right)$
$J/\psi$	$\left({1.231}_{-0.037}^{+0.038}\right)\times {10}^{-8}$	$\left({3.01}_{-0.15}^{+0.15}\right)\times {10}^{-6}$
${\rm Y}\left(1S\right)$	$(4.08{}_{-1.23}{}^{+1.65})\times {10}^{-11}$	$\left({9.97}_{-3.03}^{+4.04}\right)\times {10}^{-9}$
${\rm Y}\left(2S\right)$	$(1.07{}_{-0.37}{}^{+0.57})\times {10}^{-11}$	$\left({2.62}_{-0.91}^{+1.39}\right)\times {10}^{-9}$
${\rm Y}\left(3S\right)$	$(0.77{}_{-0.28}{}^{+0.43})\times {10}^{-11}$	$\left({1.87}_{-0.69}^{+1.05}\right)\times {10}^{-9}$

Equations

$$i\mathcal{M}(H\to V+\gamma )$$

$${\psi }^{†}{\overleftrightarrow{D}}^{i_1}\dots {\overleftrightarrow{D}}^{i_n}\chi$$

$$⟨v_S^2{⟩}_V$$

$$⟨v_S^4{⟩}_V$$

$$⟨v_D^2{⟩}_V$$

$$⟨B{⟩}_V$$

$$⟨D{E}_0{⟩}_V$$

$$⟨D{E}_1{⟩}_V$$

$$i\mathcal{M}[H\to Q\bar{Q}(J^{PC}=1^{--})+\gamma ]$$

$${\mathcal{L}}_{2-\mathrm{f}}$$

$$i\mathcal{M}[H\to Q\bar{Q}g(J^{PC}=1^{--})+\gamma ]$$

$$i{\mathcal{M}}_{\mathrm{dir}}(H\to Q\bar{Q}+\gamma )$$

$$p_{\gamma }^0=|{\mathbf{p}}_{\gamma }|=\frac{m_H^2-(p_1+p_2{)}^2}{2\sqrt{(p_1+p_2{)}^2}}=\frac{m_H^2-4m^2-4{\mathbf{q}}^2}{4\sqrt{m^2+{\mathbf{q}}^2}},$$

$$c_0$$

$$c_{{\mathbf{D}}^2}$$

$$c_{{\mathbf{D}}^{(i}{\mathbf{D}}^{j)}}$$

$$c_{{\mathbf{D}}^4}$$

$$c_{{\mathbf{D}}^2{\mathbf{D}}^{(i}{\mathbf{D}}^{j)}}$$

$$i{\mathcal{M}}_{\mathrm{dir}}(H\to Q\bar{Q}g+\gamma )$$

$$p_{\gamma }^0=|{\mathbf{p}}_{\gamma }|=\frac{m_H^2-P_0^2}{2P_0},$$

$$c_B$$

$$c_{D{E}_0}$$

$$c_{D{E}_1}$$

$$c_{D{E}_1^{\prime }}$$

$$c_{D{E}_2}$$

$$i{\mathcal{M}}_{\mathrm{dir}}(H\to V+\gamma )=\frac{i}{2}e{e}_Q{y}_Q{ϵ}_{\gamma }^{\ast }\cdot {ϵ}^{\ast }(\lambda )f_V^{\perp }(\mu ){\int }_0^{1}dx{T}_H(x,\mu ){\varphi }_V^{\perp }(x,\mu )+O(m^2/m_H^2),$$

$$f_V^{\perp }(\mu ){ϵ}_{\perp }^{\ast }{}^{\alpha }(\lambda ){\varphi }_V^{\perp }(x,\mu )=⟨V|{\mathcal{Q}}^{\alpha }(x)|0⟩.$$

$${\mathcal{Q}}^{\alpha }(x)=\int \frac{d\omega }{2\pi }{e}^{-i(x-1/2)\omega \bar{n}\cdot P}(\bar{Q}{W}_c)(\omega \bar{n}/2)\bar{n}/{\gamma }_{\perp }^{\alpha }(W_c^{†}Q)(-\omega \bar{n}/2),$$

$$T_H(x,\mu )$$

$$f_V^{\perp }({\mu }_0){\varphi }_V^{\perp }(x,{\mu }_0)$$

$$c_n=-\frac{i}{2}e{e}_Q{y}_Q{\mathbf{ϵ}}^{\ast }(\lambda )\cdot {\mathbf{ϵ}}_{\gamma }^{\ast }{\int }_0^{1}dx{T}_H(x,{\mu }_0){\tilde{c}}_n(x)+O(m^2/m_H^2).$$

$$-⟨Q\bar{Q}|{\mathcal{Q}}^{\alpha }(x)|0⟩$$

$${\tilde{c}}_0(x)$$

$${\tilde{c}}_{{\mathbf{D}}^2}(x)$$

$${\tilde{c}}_{{\mathbf{D}}^{(i}{\mathbf{D}}^{j)}}(x)$$

$${\tilde{c}}_{{\mathbf{D}}^4}(x)$$

$${\tilde{c}}_{{\mathbf{D}}^2{\mathbf{D}}^{(i}{\mathbf{D}}^{j)}}(x)$$

$$-⟨Q\bar{Q}g|{\mathcal{Q}}^{\alpha }(x)|0⟩$$

$${\tilde{c}}_B(x)$$

Full text

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 Khachatryan:2016vau , 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 Aaboud:2018txb ; Aaboud:2018fhh . One way is to measure the Higgs decay into $c\bar{c}$, by identifying charm jets in the final state. Another way is to measure the decay of the Higgs boson to a charmonium and a photon Bodwin:2013gca . Compared to $H\to c\bar{c}$, the process $H\to$ charmonium$+\gamma$ has an advantage that the charmonium provides a clean final state through its electromagnetic decays. Higgs decay into charmonium$+\gamma$ also allows a simultaneous measurement of the size and the sign of the Higgs charm coupling. The current upper limits for the branching ratio ${\rm Br}(H\to J/\psi+\gamma)$ and the cross section $\sigma(pp\to ZH)\times{\rm Br}(H\to c\bar{c})$ at 95% confidence level are both about 2 orders of magnitude larger than the Standard-Model predictions Aaboud:2018txb ; Aaboud:2018fhh .

It is crucial that the decay rate $\Gamma(H\to J/\psi+\gamma)$ 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 Bodwin:2013gca ; Bodwin:2014bpa ; Koenig:2015pha ; Bodwin:2016edd ; Bodwin:2017wdu . Especially, approaches based on nonrelativistic effective field theories allow a systematic improvement of theoretical accuracy Bodwin:2013gca ; Bodwin:2014bpa ; Bodwin:2016edd ; Bodwin:2017wdu . In the nonrelativistic QCD (NRQCD) effective field theory Bodwin:1994jh , decay and production processes involving a heavy quarkonium are given by a double series in $\alpha_{s}$ and $v$, where $v$ is the typical velocity of a heavy quark $Q$ in a heavy quarkonium; for charmonium, $v^{2}\approx 0.3$, and for bottomonium, $v^{2}\approx 0.1$. Currently, the decay rates $\Gamma(H\to V+\gamma)$ for $V=J/\psi$ or $\Upsilon(nS)$ for $n=1,2,$ and $3$ have been computed to relative order $\alpha_{s}v^{0}$ and $v^{2}$ accuracy Shifman:1980dk ; Bodwin:2014bpa ; Bodwin:2016edd ; Bodwin:2017wdu . In Refs. Bodwin:2014bpa ; Bodwin:2016edd ; Bodwin:2017wdu , the large logarithms of $m_{H}^{2}/m_{V}^{2}$ that appear in higher order corrections in $\alpha_{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 $\alpha_{s}$ by combining the NRQCD and the light-cone formalisms Lepage:1980fj ; Chernyak:1983ej ; Jia:2008ep .

In this paper, we improve the accuracy of the Standard Model prediction of the decay rates $\Gamma(H\to V+\gamma)$ for $V=J/\psi$ or $\Upsilon(nS)$ for $n=1,2,$ and $3$ by computing the order-$v^{4}$ correction to the decay rate in the NRQCD factorization formalism. In our numerical analysis, we do not consider the $\psi(2S)$ meson, because, to date, there are no available estimates of the relevant NRQCD matrix elements accounting for open-flavor threshold effects and nonrelativistic corrections in a complete and model-independent way. These effects may be particularly important for this state, as it is just $43$ MeV below the $D\bar{D}$ threshold. We work in the limit $m_{V}^{2}/m_{H}^{2}\to 0$, where the calculation simplifies dramatically. In this limit, $H\to V+\gamma$ 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 $\bar{Q}$ through the Yukawa interaction, and the $Q\bar{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_{H}^{2}/m_{V}^{2}$ to all orders in $\alpha_{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 $\alpha_{s}^{n}\log^{n}(m_{H}^{2}/m_{V}^{2})$ and $\alpha_{s}^{n}\log^{n-1}(m_{H}^{2}/m_{V}^{2})$, respectively, for all orders $n\geq 1$ in $\alpha_{s}$. We note that the light-cone formalism applies only to the leading-order piece in the expansion in powers of $m_{V}^{2}/m_{H}^{2}$ Lepage:1980fj ; Chernyak:1983ej . In the indirect process, the Higgs boson first decays into a $\gamma$ and a $\gamma^{*}$, and the $\gamma^{*}$ evolves into a quarkonium Bodwin:2013gca . 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_{V}^{2}/m_{H}^{2}\to 0$, we find simplifications in the indirect process that let us compute the indirect amplitude accurately from the known calculation of the $H\to\gamma\gamma$ decay amplitude and the leptonic decay rate of the meson $V$. Also, the logarithms of $m_{H}^{2}/m_{V}^{2}$ 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}\approx 0.005$, which indeed makes the direct amplitude numerically smaller than the indirect one. In the case of the bottom quark $y_{b}\approx 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 $\alpha_{s}$ and order $v^{2}$ corrections and resum leading and next-to-leading logarithms of $m_{H}^{2}/m_{V}^{2}$ to all orders in $\alpha_{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 $\alpha_{s}$, but we will include the previously known order $\alpha_{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 $\epsilon(\lambda)$ and a photon to relative order $v^{4}$ accuracy is given by

[TABLE]

Here, $m$ is the mass of the heavy quark, $g_{s}$ is the strong coupling, $\psi^{\dagger}$ and $\chi$ are Pauli spinor fields that create a heavy quark and an antiquark, respectively, and $E^{i}=G^{i0}$ and $B^{i}=\frac{1}{2}\epsilon^{ijk}G^{kj}$ are chromoelectric and chromomagnetic fields, respectively, where $G^{\mu\nu}$ is the gluon field-strength tensor. The covariant derivative $\bm{D}=\bm{\nabla}-ig_{s}\bm{A}$ appears in Eq. (1) in the combination $\psi^{\dagger}\overleftrightarrow{\bm{D}}\chi=\psi^{\dagger}\bm{D}\chi-(\bm{D}\psi)^{\dagger}\chi$. Operators with more than one covariant derivative are defined with

[TABLE]

The notation $T^{(ij)}=\tfrac{1}{2}(T^{ij}+T^{ji})-\tfrac{1}{3}T^{ii}\delta^{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\rangle$, while the long-distance matrix elements (LDMEs) of NRQCD operators between the vacuum $|0\rangle$ and the meson state $|V\rangle$ are nonperturbative quantities. We take the meson state $|V\rangle$ 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 $\langle 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/\psi$ or $\Upsilon(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. Brambilla:2001xy ; Brambilla:2002nu ; Brambilla:2006ph ; Brambilla:2008zg . 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\bar{Q}$ Fock state created by the operator. For $J/\psi$ or $\Upsilon(nS)$, the leading Fock state contains a $Q\bar{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\bar{Q}$ in a color-octet state, or the ones that contain $Q\bar{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 $\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi$ is the lowest-dimensional operator that creates a $Q\bar{Q}$ in the leading Fock state (${}^{3}S_{1}$), and so, the LDME $\langle V|\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle$ scales like $v^{3/2}$, and contributes to the amplitude at leading order in $v$. The LDMEs of the operators $\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)(-\tfrac{i}{2}\overleftrightarrow{\bm{D}})^{2}\chi$ and $\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)(-\tfrac{i}{2}\overleftrightarrow{\bm{D}})^{4}\chi$ scale like $v^{7/2}$ and $v^{11/2}$, respectively, because the operators create the $Q\bar{Q}$ in the leading Fock state and have dimensions that are higher than the lowest-dimensional operator by 2 and 4, respectively. Hence, these LDMEs contribute to the amplitude at relative order $v^{2}$ and $v^{4}$, respectively. The operator $\psi^{\dagger}\epsilon^{i}(\lambda)\sigma^{j}(-\tfrac{i}{2})^{2}\overleftrightarrow{\bm{D}}^{(i}\overleftrightarrow{\bm{D}}^{j)}\chi$ creates a color-singlet $Q\bar{Q}$ in a ${}^{3}D_{1}$ state. Since the $D$-wave Fock state has a contribution to the meson state suppressed by $v^{2}$ compared to the leading Fock state, the LDME $\langle V|\psi^{\dagger}\epsilon^{i}(\lambda)\sigma^{j}(-\tfrac{i}{2})^{2}\overleftrightarrow{\bm{D}}^{(i}\overleftrightarrow{\bm{D}}^{j)}\chi|0\rangle$ scales like $v^{11/2}$ and contributes to the amplitude at relative order $v^{4}$. The color-octet operators in Eq. (1) create $Q\bar{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 $\langle V|\psi^{\dagger}g_{s}\bm{B}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle$, $\langle V|\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\tfrac{1}{3}(\overleftrightarrow{\bm{D}}\cdot g_{s}\bm{E}+g_{s}\bm{E}\cdot\overleftrightarrow{\bm{D}})\chi|0\rangle$, and $\langle V|\psi^{\dagger}\bm{\epsilon}(\lambda)\cdot\tfrac{1}{2}[\bm{\sigma}\times(\overleftrightarrow{\bm{D}}\times g_{s}\bm{E}-g_{s}\bm{E}\times\overleftrightarrow{\bm{D}})]\chi|0\rangle$ 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:

[TABLE]

There is a color-singlet operator of dimension 7 that does not appear in Eq. (1), which is given by $\frac{1}{2}\psi^{\dagger}\epsilon^{i}(\lambda)\sigma^{j}(-\tfrac{i}{2})^{2}\{\overleftrightarrow{\bm{D}}^{(i}\overleftrightarrow{\bm{D}}^{j)},(-\tfrac{i}{2}\overleftrightarrow{\bm{D}})^{2}\}\chi$. Because this operator creates a $Q\bar{Q}$ in a ${}^{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 $\psi^{\dagger}\bm{\epsilon}(\lambda)\cdot\frac{i}{2}(\overleftrightarrow{\bm{D}}\times g_{s}\bm{E}+g_{s}\bm{E}\times\overleftrightarrow{\bm{D}})\chi$ and $\psi^{\dagger}\epsilon^{i}(\lambda)\sigma^{j}(\overleftrightarrow{\bm{D}}^{(i}g_{s}\bm{E}^{j)}+g_{s}\bm{E}^{(i}\overleftrightarrow{\bm{D}}^{j)})\chi$ 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. (65) and (65).

If we follow the power counting of Ref. Bodwin:1994jh , the color-octet LDMEs except for $\langle V|\psi^{\dagger}\tfrac{1}{3}(\overleftrightarrow{\bm{D}}\cdot g_{s}\bm{E}+g_{s}\bm{E}\cdot\overleftrightarrow{\bm{D}})\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle$ 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 $\alpha_{s}$ by using the perturbative matching conditions obtained by replacing the meson state $V$ with a perturbative $Q\bar{Q}$ or a $Q\bar{Q}g$ 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$, $\bar{Q}$, and the gluon, and truncate the series to the desired accuracy. We follow a method used in Refs. Braaten:1996rp ; Brambilla:2017kgw , that consists in not projecting to a specific color, spin or orbital angular momentum of the $Q\bar{Q}$ state, but instead, in only requiring the $Q\bar{Q}$ or the $Q\bar{Q}g$ state to have the same $J^{PC}=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 the 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\bar{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^{PC}=1^{--}$.

If we replace the meson state $V$ with a perturbative $Q\bar{Q}$ state, the amplitude occurs from order $g_{s}^{0}$, 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\bar{Q}$ state at up to fourth order in the relative momentum of the $Q$ and the $\bar{Q}$. We use the kinematical configuration given in Appendix A where the relative 3-momentum between the $Q$ and the $\bar{Q}$ is given by $\bm{q}$. The production amplitude of a $Q\bar{Q}$ state with $J^{PC}=1^{--}$ and a photon is given at order $g_{s}^{0}$ by

[TABLE]

where we have included all color-singlet operators with $J^{PC}=1^{--}$ up to dimension 7. We take the states $|Q\rangle$ and the $|\bar{Q}\rangle$ to be nonrelativistically normalized. We determine the short-distance coefficients $c_{0}$, $c_{\bm{D}^{2}}$, $c_{\bm{D}^{(i}\bm{D}^{j)}}$, and $c_{\bm{D}^{4}}$, along with $c_{\bm{D}^{2}\bm{D}^{(i}\bm{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 $\bm{q}$ up to fourth order.

In order to compute the remaining short-distance coefficients corresponding to the color-octet LDMEs, we consider the production amplitude of a $Q\bar{Q}g$ state with $J^{PC}=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 $\bar{Q}$ is given by $\bm{q}_{1}$ and the relative 3-momentum between the $Q\bar{Q}$ pair and the gluon is given by $\bm{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 $\bm{q}_{1}$ or $\bm{q}_{2}$ when the meson state $V$ is replaced by the $Q\bar{Q}g$ state. We must also include all color-octet operators of dimensions up to 6 with $J^{PC}=1^{--}$ that do not appear in Eq. (1), whose matrix elements can also be either linear or quadratic in the momenta $\bm{q}_{1}$ or $\bm{q}_{2}$. The color-singlet operators in Eq. (1) can also contribute to the $Q\bar{Q}g$ 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/|\bm{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 Manohar:1997qy

[TABLE]

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 $\alpha_{s}^{0}$ in Eq. (5). The production amplitude of a $Q\bar{Q}g$ state with $J^{PC}=1^{--}$ and a photon at order $g_{s}$ is given by

[TABLE]

where the left-hand side is calculated in perturbative QCD and is expanded in powers of $\bm{q}_{1}$ and $\bm{q}_{2}$ up to quadratic accuracy. We again take the states $|Q\rangle$ and $|\bar{Q}\rangle$ to be nonrelativistically normalized. We determine the short-distance coefficients $c_{B}$, $c_{DE_{0}}$, and $c_{DE_{1}}$, along with $c_{DE_{1}^{\prime}}$ 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 $\bm{q}_{1}$ and $\bm{q}_{2}$ up to quadratic order.111The operator $\psi^{\dagger}\bm{\epsilon}(\lambda)\cdot\tfrac{i}{2}(\overleftrightarrow{\bm{D}}\times g_{s}\bm{E}+g_{s}\bm{E}\times\overleftrightarrow{\bm{D}})\chi$ does not appear in Ref. Brambilla:2006ph . 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 $\bm{q}_{2}$. Hence, if we ignore the contribution that is proportional to $|\bm{q}_{2}|^{2}$ on 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).

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 left-hand sides of Eqs. (4) and (6) are computed at leading order in $\alpha_{s}$. We obtain corrections 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_{H}^{2}$, that allows us to resum logarithms of $m_{H}^{2}/m^{2}$ to all orders in $\alpha_{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 $\alpha_{s}$ correction to the direct amplitude using the light-cone approach.

II.1 Fixed-order calculation

At order $g_{s}^{0}$, the direct amplitude for $H\to Q\bar{Q}+\gamma$ is given by

[TABLE]

where $p_{\gamma}$ and $\epsilon_{\gamma}^{*}$ 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 $\epsilon_{\gamma}^{*}\cdot p_{\gamma}=0$. Here, $e=\sqrt{4\pi\alpha}$ is the electric charge, $e_{Q}$ is the fractional charge of the heavy quark $Q$, and $y_{Q}=\overline{m}(\mu)(\sqrt{2}G_{F})^{\frac{1}{2}}$ is the Yukawa coupling of the Higgs boson and $Q$, with $G_{F}$ the Fermi constant. The momenta of the $Q$ and $\bar{Q}$ are given by $p_{1}$ and $p_{2}$, respectively, so that the momentum of the $H$ is $P_{H}=p_{1}+p_{2}+p_{\gamma}$. This implies $m_{H}^{2}=(p_{1}+p_{2}+p_{\gamma})^{2}=2p_{\gamma}\cdot(p_{1}+p_{2})+(p_{1}+p_{2})^{2}$, so that in the rest frame of the $Q\bar{Q}$,

[TABLE]

where $q=\frac{1}{2}(p_{1}-p_{2})$. We choose the heavy-quark mass appearing in $y_{Q}$ to be $\overline{m}(\mu)$, which is the $\overline{\rm MS}$ mass of the heavy quark $Q$ at scale $\mu$; as we will see in the next section, this choice simplifies the logarithms that appear in the order $\alpha_{s}$ correction. Since $C$, $P$, and $T$ are conserved in the amplitude in Eq. (7), the $Q\bar{Q}$ can only be created with $C=-1$. We first express the Dirac bilinears $\bar{u}(p_{1})\gamma^{\mu}v(p_{2})$ and $\bar{u}(p_{1})\gamma^{\mu}\gamma^{\nu}v(p_{2})$ in terms of the 3-momenta of the $Q$ and $\bar{Q}$ in the $Q\bar{Q}$ rest frame. This can be accomplished by using the method described in Appendix B. Then, we expand the amplitude in powers of $\bm{q}$, keeping terms up to relative order $\bm{q}^{4}/m^{4}$. The resulting expression for the amplitude is then a linear combination of the Cartesian tensors built from $\xi^{\dagger}\bm{\sigma}\eta$ and $\bm{q}$ of the form $\xi^{\dagger}\sigma^{i}\eta q^{j}\cdots q^{k}$ up to rank 5 and of the form $\xi^{\dagger}\eta q^{i}q^{j}\cdots 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. CoopeSnider . Finally, the $P=-1$ contribution is obtained by keeping only the contribution odd in parity, where the parity transform of the $Q\bar{Q}$ amplitude is given by the replacements $\bm{q}\to-\bm{q}$, $\xi^{\dagger}\bm{\sigma}\eta\to-\xi^{\dagger}\bm{\sigma}\eta$, and $\xi^{\dagger}\eta\to-\xi^{\dagger}\eta$. We use the Mathematica package FeynCalc Mertig:1990an ; Shtabovenko:2016sxi and the FeynOnium FeynOnium package to automatize the calculation of the amplitude and the consequent reduction to the $J^{PC}=1^{--}$ contribution. By comparing the $J^{PC}=1^{--}$ contribution of the $Q\bar{Q}$ amplitude with the right-hand side of Eq. (4), we obtain the short-distance coefficients

[TABLE]

where we define $r\equiv\frac{4m^{2}}{m_{H}^{2}}$. The short-distance coefficients $c_{0}$ and $c_{\bm{D}^{2}}$ agree with Refs. Bodwin:2013gca ; Bodwin:2014bpa , except that our results differ by an overall sign that originates from the sign convention of the $J=1$ state employed in Refs. Bodwin:2013gca ; Bodwin:2014bpa . We also reproduce the short-distance coefficient $c_{\bm{D}^{4}}$ that can be obtained from Ref. Bodwin:2014bpa . The results for $c_{\bm{D}^{(i}\bm{D}^{j)}}$ and $c_{\bm{D}^{2}\bm{D}^{(i}\bm{D}^{j)}}$ are new.

The remaining short-distance coefficients corresponding to the color-octet LDMEs are computed from the direct amplitude for $H\to Q\bar{Q}g+\gamma$ at order $g_{s}$, which is given by

[TABLE]

where $\epsilon_{g}^{*}$ is the polarization vector of the gluon. The total momentum $P$ of the $Q\bar{Q}g$ state is given by $P=p_{1}+p_{2}+k_{g}$, and the momentum of the $H$ is given by $P_{H}=P+p_{\gamma}$, so that in the rest frame of the $Q\bar{Q}g$,

[TABLE]

where $P_{0}=2|\bm{q}_{2}|+\sqrt{(\bm{q}_{1}+\bm{q}_{2})^{2}+m^{2}}+\sqrt{(\bm{q}_{1}-\bm{q}_{2})^{2}+m^{2}}$, $q_{1}=\frac{1}{2}(p_{1}-p_{2})$, and $q_{2}=\frac{1}{6}(2k_{g}-p_{1}-p_{2})$. Due to $C$ conservation at this order, the $Q\bar{Q}g$ can only be produced in a color-singlet $C=-1$ state. We use again the Mathematica package FeynCalc and the FeynOnium package to automatize the calculation of the $Q\bar{Q}g$ amplitude. After expressing the amplitude in terms of the 3-momenta of the $Q$, $\bar{Q}$ and $g$, we expand the amplitude in powers of $\bm{q}_{1}$ and $\bm{q}_{2}$ up to relative order $\bm{q}_{1}^{2}/m^{2}$, $\bm{q}_{2}^{2}/m^{2}$, and $|\bm{q}_{1}||\bm{q}_{2}|/m^{2}$ The resulting expression for the amplitude is then a linear combination of the Cartesian tensors built from $\xi^{\dagger}\bm{\sigma}\eta$, $\bm{\epsilon}_{g}^{*}$, $\bm{q}_{1}$, and $\bm{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\bar{Q}g$ amplitude is given by the replacements $\bm{q}_{1}\to-\bm{q}_{1}$, $\bm{q}_{2}\to-\bm{q}_{2}$, $\bm{\epsilon}_{g}^{*}\to-\bm{\epsilon}_{g}^{*}$, $\xi^{\dagger}\bm{\sigma}\eta\to-\xi^{\dagger}\bm{\sigma}\eta$, and $\xi^{\dagger}\eta\to-\xi^{\dagger}\eta$. After reducing the Cartesian tensors of odd rank to total angular momentum $J=1$, we compare the $J^{PC}=1^{--}$ contribution with the amplitude on the right-hand side of Eq. (6) to obtain the short-distance coefficients

[TABLE]

where $r=\frac{4m^{2}}{m_{H}^{2}}$.

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_{H}^{2}$ (see Sec. III), only the limit $m^{2}/m_{H}^{2}\to 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^{PC}$. For example, the $J^{PC}=1^{+-}$ contribution can be obtained by keeping parity-even terms in the amplitude. The $J^{PC}=1^{+-}$ contribution can then be used to obtain the short-distance coefficients for the $H\to h_{c}+\gamma$ amplitude. We show this result in Appendix C.

II.2 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_{H}^{2}/m^{2}$ that appear in corrections of higher orders in $\alpha_{s}$ can potentially spoil the convergence of the perturbation series. If we work at leading power in $m^{2}/m_{H}^{2}$, the light-cone approach provides a factorization formula that separates the contribution at the scale $m$ from the contribution at the scale $m_{H}$ Lepage:1980fj ; Chernyak:1983ej . The light-cone approach also enables us to resum the logarithms of $m_{H}^{2}/m^{2}$ to all orders in $\alpha_{s}$ by solving an evolution equation. The factorization formula is given by Lepage:1980fj ; Chernyak:1983ej ; Wang:2013ywc

[TABLE]

where $T_{H}(x,\mu)$ 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}^{\perp}(\mu)$ and the light-cone distribution amplitude (LCDA) $\phi_{V}^{\perp}(x,\mu)$ of the meson $V$. The decay constant and the LCDA are nonperturbative quantities defined by

[TABLE]

The nonlocal operator ${\cal Q}^{\alpha}(x)$ is defined by

[TABLE]

where $P$ is the momentum of the meson $V$, and the decay constant $f_{V}^{\perp}(\mu)$ is defined by integrating Eq. (14) over $x$ and considering the normalization of the LCDA, which is given by $\int_{0}^{1}dx\,\phi_{V}^{\perp}(x,\mu)=1$. Here, $Q(x)$ is the QCD quark field. The light-cone vectors $n$ and $\bar{n}$ are given by $\bar{n}=\frac{m_{H}}{p_{\gamma}\cdot P_{H}}p_{\gamma}$ and $n=\frac{2}{m_{H}}P_{H}-\bar{n}$, which satisfy $n^{2}=\bar{n}^{2}=0$ and $\bar{n}\cdot n=2$. Here, $P_{H}$ and $p_{\gamma}$ are the momenta of the $H$ and the photon, respectively. For any 4-vector $a^{\mu}$, we define $a^{\mu}_{\perp}\equiv a^{\mu}-\frac{n^{\mu}}{2}\bar{n}\cdot a-\frac{\bar{n}^{\mu}}{2}{n}\cdot a$. The Wilson line $W_{c}(x)={\cal P}\exp\left[-ig\int_{-\infty}^{0}ds\,\bar{n}\cdot A(x+s\bar{n})\right]$, where ${\cal P}$ is the path-ordering operator, ensures the gauge invariance of the nonlocal operator. The hard part $T_{H}(x,\mu)$ has been computed to next-to-leading order (NLO) accuracy in $\alpha_{s}$ and is given by Wang:2013ywc ; Koenig:2015pha

[TABLE]

where $C_{F}=\frac{N_{c}^{2}-1}{2N_{c}}$, and $N_{c}=3$ is the number of colors. The imaginary part in Eq. (16) comes from the order $\alpha_{s}$ correction where the virtual lines can be on shell. The expression for the NLO correction to $T_{H}(x,\mu)$ in Eq. (16) is valid when the decay constant $f_{V}^{\perp}(\mu)$ and the LCDA $\phi_{V}^{\perp}(x,\mu)$ are renormalized in the $\overline{\rm MS}$ scheme, and the heavy-quark mass in the Yukawa coupling $y_{Q}$ is the $\overline{\rm MS}$ mass at the scale $\mu$. If we use the heavy-quark pole mass $m$ instead of the $\overline{\rm MS}$ mass at the scale $\mu$, the order $\alpha_{s}$ correction in $T_{H}(x,\mu)$ involves a logarithm of $m_{H}^{2}/m^{2}$; instead, using the $\overline{\rm 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 Koenig:2015pha .

The nonlocal operator ${\cal Q}^{\alpha}(x)$ has an anomalous dimension known to NLO in $\alpha_{s}$ that allows us to resum logarithms of $m_{H}^{2}/m^{2}$ to NLL accuracy Lepage:1980fj ; Mueller:1993hg ; Mueller:1994cn ; Vogelsang:1997ak ; Hayashigaki:1997dn . In order to resum the logarithms of $m_{H}^{2}/m^{2}$ that appear in corrections of higher orders in $\alpha_{s}$ to the short-distance coefficients, we apply the factorization formula [Eq. (13)] to the perturbative amplitudes $H\to Q\bar{Q}+\gamma$ and $H\to Q\bar{Q}g+\gamma$. This involves calculating the decay constant and the LCDA with the meson state $V$ replaced by the perturbative $Q\bar{Q}$ and $Q\bar{Q}g$ 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 Jia:2008ep

[TABLE]

where the short-distance coefficients $\tilde{c}_{n}(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 $\tilde{c}_{n}(x)$, and the left-hand sides are replaced by $-\langle Q\bar{Q}(J^{PC}=1^{--})|{\cal Q}^{\alpha}(x)|0\rangle$ and $-\langle Q\bar{Q}g(J^{PC}=1^{--})|{\cal Q}^{\alpha}(x)|0\rangle$, respectively. In the matching conditions, the scale $\mu_{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

[TABLE]

It is worth noting that $\langle Q\bar{Q}|{\cal Q}^{\alpha}(x)|0\rangle$ and $\langle Q\bar{Q}g|{\cal Q}^{\alpha}(x)|0\rangle$ contain contributions from both negative and positive charge conjugation. From the fact that the charge conjugate of the operator ${\cal Q}^{\alpha}(x)$ is given by $-{\cal Q}^{\alpha}(1-x)$, we can see that the contributions of negative charge conjugation to $\langle Q\bar{Q}|{\cal Q}^{\alpha}(x)|0\rangle$ and $\langle Q\bar{Q}g|{\cal Q}^{\alpha}(x)|0\rangle$ are given by the contribution symmetric in $x\leftrightarrow 1-x$. Since the hard part $T_{H}(x,\mu_{0})$ is symmetric in $x\leftrightarrow 1-x$, the positive charge conjugation contribution of the LCDA, which is antisymmetric in $x\leftrightarrow 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 $\langle Q\bar{Q}|{\cal Q}^{\alpha}(x)|0\rangle$ and $\langle Q\bar{Q}g|{\cal Q}^{\alpha}(x)|0\rangle$, we just need to keep contributions that are symmetric in $x\leftrightarrow 1-x$.

If we replace the meson state $V$ with the $Q\bar{Q}$ state, we obtain

[TABLE]

We use the same strategy as the fixed-order calculation to obtain the contribution with $J^{PC}=1^{--}$; we express the Dirac bilinears in terms of the 3-momenta of the $Q$ and the $\bar{Q}$ in the $Q\bar{Q}$ rest frame, and then we expand Eq. (19) in powers of $\bm{q}$ up to relative order $\bm{q^{4}}/m^{4}$. Then, Eq. (19) is given by a linear combination of the delta function $\delta(x-1/2)$ and its derivatives. In order to keep only the contribution with $C=-1$, we ignore the odd derivatives of $\delta(x-1/2)$. The $J=1$ contribution is then obtained by reducing the Cartesian tensors of the form $\xi^{\dagger}\eta q^{i}q^{j}\cdots q^{k}$ up to rank 4 and of the form $\xi^{\dagger}\sigma^{i}\eta q^{j}\cdots q^{k}$ up to rank 5 using the reduction method of Ref. CoopeSnider . We keep only the contribution with negative parity in order to obtain the $J^{PC}=1^{--}$ contribution. From the matching condition (17) we obtain

[TABLE]

The short-distance coefficients $\tilde{c}_{0}(x)$ and $\tilde{c}_{\bm{D}^{2}}(x)$ agree with known results in Ref. Bodwin:2014bpa . 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\bar{Q}g$ state to obtain

[TABLE]

After expressing the amplitude in terms of the 3-momenta of the $Q$, $\bar{Q}$ and $g$, we expand the amplitude in powers of $\bm{q}_{1}$ and $\bm{q}_{2}$ up to relative order $\bm{q}_{1}^{2}/m^{2}$, $\bm{q}_{2}^{2}/m^{2}$, and $|\bm{q}_{1}||\bm{q}_{2}|/m^{2}$. The resulting expression for the amplitude is then a linear combination of the Cartesian tensors built from $\xi^{\dagger}\bm{\sigma}\eta$, $\bm{\epsilon}_{g}^{*}$, $\bm{q}_{1}$, and $\bm{q}_{2}$ up to rank 4. We obtain the matching condition for $Q\bar{Q}g$ with $J^{PC}=1^{--}$ by keeping only the contributions symmetric in $x\leftrightarrow 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

[TABLE]

The agreement with the fixed-order calculation can be easily verified using Eq. (18).

By integrating the short-distance coefficients $\tilde{c}_{n}(x)$ over $x$, we obtain an expression for $f_{V}^{\perp}$ valid up to relative order $v^{4}$. If we include the correction of relative order $\alpha_{s}v^{0}$ in the $\overline{\rm MS}$ scheme from Ref. Wang:2013ywc , we obtain

[TABLE]

The corrections at relative order $v^{2}$ agree with Ref. Bodwin:2014bpa . The logarithm of $\mu_{0}^{2}/m^{2}$ in Eq. (23) is the remnant of the renormalization of the decay constant in the $\overline{\rm MS}$ scheme.

Since the short-distance coefficients $\tilde{c}_{n}(x)$ at leading order in $\alpha_{s}$ are linear combinations of $\delta(x-1/2)$, $\delta^{(2)}(x-1/2)$ and $\delta^{(4)}(x-1/2)$, the LCDA can be written as

[TABLE]

where

[TABLE]

and $\phi_{V}^{\perp}{}^{(\alpha_{s})}(x,\mu_{0})$ is the order $\alpha_{s}v^{0}$ correction in the $\overline{\rm MS}$ scheme given by Wang:2013ywc

[TABLE]

where the plus and plus-plus distributions are defined by

[TABLE]

The logarithm of $\mu_{0}^{2}/m^{2}$ in Eq. (II.2) is the remnant of the renormalization of the LCDA in the $\overline{\rm MS}$ scheme. Note that

[TABLE]

The decay constant and the LCDA computed in Eqs. (23, 24) allow us to resum logarithms of $m_{H}^{2}/m^{2}$ that appear in the QCD corrections to the direct amplitude to all orders in $\alpha_{s}$. This is accomplished by solving the evolution equation

[TABLE]

where the evolution kernel $V_{T}(x,y)$ is currently known to NLO accuracy in $\alpha_{s}$ Lepage:1980fj ; Mueller:1993hg ; Mueller:1994cn ; Vogelsang:1997ak ; Hayashigaki:1997dn . We obtain $f_{V}^{\perp}(\mu)\phi_{V}^{\perp}(x,\mu)$ at scale $\mu$ from $f_{V}^{\perp}(\mu_{0})\phi_{V}^{\perp}(x,\mu_{0})$ at scale $\mu_{0}$ by solving this evolution equation. The direct amplitude with logarithms of $m_{H}^{2}/m^{2}$ resummed to all orders in $\alpha_{s}$ is then given by Eq. (13), by setting $\mu\sim m_{H}$ and $\mu_{0}\sim m$, so that $T_{H}(x,\mu)$ is free of logarithms of $m_{H}^{2}/m^{2}$.222 Equivalently, one may resum logarithms of $m_{H}^{2}/m^{2}$ in the hard part $T_{H}(x,\mu)$, which may be conceptually closer to the effective field theory logic.

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_{n}^{(3/2)}(2x-1)$, which are the eigenfunctions of the LO evolution kernel Brodsky:1980ny . The convolution in Eq. (13) is given by

[TABLE]

where $\hat{T}_{H}(n,\mu)$ and $\hat{\phi}_{V}^{\perp}(n,\mu)$ are Gegenbauer moments defined by

[TABLE]

The solution of the evolution equation in terms of Gegenbauer moments leads to Mueller:1993hg ; Mueller:1994cn ; Vogelsang:1997ak ; Hayashigaki:1997dn

[TABLE]

Explicit expressions of $U_{f_{V}^{\perp}}(\mu,\mu_{0})$ and $U_{n_{1}n_{2}}(\mu,\mu_{0})$ can be found in Ref. Bodwin:2016edd . We note that $U_{f_{V}}(\mu,\mu_{0})$ and $U_{n_{1}n_{2}}(\mu,\mu_{0})$ depend only on $\mu$, $\mu_{0}$ and the evolution kernel, and are independent of $f_{V}^{\perp}(\mu_{0})$ or $\phi_{V}^{\perp}(x,\mu_{0})$. When $\phi_{V}^{\perp}(x,\mu_{0})$ contains singular distributions such as the delta function and their derivatives, the sum in Eq. (32a) can converge badly. The authors of Refs. Bodwin:2016edd ; Bodwin:2017wdu 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

[TABLE]

The numerical value of the Abel sum is then evaluated by using Padé approximants of the truncated series. We employ this Abel-Padé method Bodwin:2016edd ; Bodwin:2017wdu 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

[TABLE]

where

[TABLE]

III Calculation of the indirect amplitude

The indirect amplitude proceeds from $H\to\gamma^{*}\gamma$, followed by $\gamma^{*}\to V$. Since we work in the limit $m_{V}^{2}/m_{H}^{2}\to 0$, the partial amplitude $H\to\gamma^{*}\gamma$ can be replaced by the decay amplitude of the Higgs boson to two photons. Then, the indirect amplitude is given by Bodwin:2013gca

[TABLE]

where the factor $16\pi m_{H}$ compensates for the factors that are necessary to obtain the decay rate $\Gamma(H\to\gamma\gamma)$ 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 $\mu_{0}$, and the scale of the vertex for the real photon to be [math]. The factor $\alpha(\mu_{0})/\alpha(0)$ replaces one QED coupling constant at scale 0 in the $H\to\gamma\gamma$ amplitude with the QED coupling constant at scale $\mu_{0}$. We ignore the small imaginary part in the $H\to\gamma\gamma$ amplitude Koenig:2015pha ; Bodwin:2016edd . The decay constant $f_{V}^{\parallel}$ is defined by

[TABLE]

Since $f_{V}^{\parallel}$ is a conserved current in QCD, it does not undergo renormalization. Hence, the indirect amplitude is free of logarithms of $m_{H}^{2}/m^{2}$ in the limit $m_{V}^{2}/m_{H}^{2}\to 0$ if one ignores the higher-order electroweak corrections to the indirect amplitude. We can express $f_{V}^{\parallel}$ 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}^{\perp}$. We obtain

[TABLE]

where we included the order $\alpha_{s}v^{0}$ correction computed in Refs. Barbieri:1975ki ; Celmaster:1978yz . We note that $f_{V}^{\parallel}$ is positive at leading order in $\alpha_{s}$ and $v$. If we assume $f_{V}^{\parallel}>0$, an accurate numerical value for $f_{V}^{\parallel}$ can be obtained from the leptonic decay rate

[TABLE]

so that

[TABLE]

Note that $i{\cal M}_{\rm ind}$ and $i{\cal M}_{\rm dir}$ have opposite signs, so that when calculating the decay rate, the direct and indirect amplitudes interfere destructively Bodwin:2013gca .

IV Numerical Results

We now present our numerical results for the Higgs decay rate into $V+\gamma$ for $V=J/\psi$ and $\Upsilon(nS)$ for $n=1,2$, and 3 based on our calculation of the direct amplitude to relative order $v^{4}$ accuracy. As we have mentioned in Sec. I, we do not consider the $\psi(2S)$ state, because there are no available estimates of the relevant NRQCD matrix elements that account for open-flavor threshold effects and nonrelativistic corrections in a complete and model-independent way.

Our expressions for the direct and indirect amplitudes are given in Eqs. (34) and (40), respectively. If we write ${\cal M}_{\rm dir}=\epsilon_{\gamma}^{*}\cdot\epsilon^{*}(\lambda){\cal A}_{\rm dir}/\sqrt{\Phi}$ and ${\cal M}_{\rm ind}=\epsilon_{\gamma}^{*}\cdot\epsilon^{*}(\lambda){\cal A}_{\rm ind}/\sqrt{\Phi}$, the decay rate $\Gamma(H\to V+\gamma)$ is given by

[TABLE]

where $\Phi$ is the phase-space and normalization factor given by

[TABLE]

We present our numerical results for ${\cal A}_{\rm dir}$ and ${\cal A}_{\rm ind}$ in the following sections. Then, using the resulting values of ${\cal A}_{\rm dir}$ and ${\cal A}_{\rm ind}$, we compute the decay rate $\Gamma(H\to V+\gamma)$ from the formula in Eq. (41).

IV.1 Indirect amplitude

We compute the numerical values for the indirect amplitude using Eq. (40). We compute the decay constant $f_{V}^{\parallel}$ from the measured leptonic decay rate:

[TABLE]

where we take $\mu_{0}=m_{V}$. Then

[TABLE]

We note that the resulting expression for ${\cal A}_{\rm ind}$ does not depend on $\alpha(\mu_{0})$.

We use the following input parameters to compute ${\cal A}_{\rm ind}$ numerically. We take the PDG value for the Higgs mass $m_{H}=125.18\pm 0.16$ Tanabashi:2018oca , and the numerical value for the Higgs two-photon decay rate to be $\Gamma(H\to\gamma\gamma)=9.34\times 10^{-6}$ GeV, which is computed from the tabulated results of the two-photon branching ratio and the total decay rate of the Higgs boson in Refs. Dittmaier:2011ti ; Dittmaier:2012vm . We also take the PDG values for the meson masses $m_{V}$ and use the measured partial widths $\Gamma(V\to e^{+}e^{-})$ for the leptonic decay rates Tanabashi:2018oca . We list the values for $m_{V}$ and $\Gamma(V\to\ell^{+}\ell^{-})$ that we use to compute the indirect amplitude in Table. 1. The QED coupling constant at scale 0 is taken to be $\alpha(0)=1/137.036$.

We consider the following sources of uncertainties in ${\cal A}_{\rm ind}$. The uncertainty in the decay rate $\Gamma(H\to\gamma\gamma)$ is taken to be $0.01$ times the central value, as estimated in Ref. Dittmaier:2011ti from the higher-order corrections to the decay rate. We consider the experimental uncertainties in $\Gamma(V\to\ell^{+}\ell^{-})$. We estimate the uncalculated correction of relative order $m_{V}^{2}/m_{H}^{2}$ to be $m_{V}^{2}/m_{H}^{2}$ 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 ${\cal A}_{\rm ind}$ are shown in Table. 2. We note that the uncertainties in ${\cal A}_{\rm ind}$ are less than 2% of the central values.

We can compare our results for ${\cal A}_{\rm ind}$ with a previous calculation in Ref. Bodwin:2017wdu . Our calculation is equivalent to the one in Ref. Bodwin:2017wdu , except that we use an updated value of the measured Higgs mass from Ref. Tanabashi:2018oca , which has smaller uncertainties than what was employed in Ref. Bodwin:2017wdu . Our results for ${\cal A}_{\rm ind}$ in Table. 2 are compatible with those in Ref. Bodwin:2017wdu within uncertainties.

IV.2 Direct amplitude

We compute the numerical values of the direct amplitude using Eq. (34), so that

[TABLE]

We now discuss our strategy to compute Eq. (45) numerically. The decay constant $f_{V}^{\perp}(\mu_{0})$ depends on the LDME $\langle V|\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle$ and the ratios of LDMEs $\langle v_{S}^{2}\rangle_{V}$, $\langle v_{S}^{4}\rangle_{V}$, $\langle v_{D}^{2}\rangle_{V}$, $\langle DE_{0}\rangle_{V}$, $\langle B\rangle_{V}$, and $\langle DE_{1}\rangle_{V}$, see Eq. (23). The dependence on the leading-order LDME $\langle V|\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle$ can be eliminated by rewriting the decay constant as

[TABLE]

and by obtaining the numerical value for $f_{V}^{\parallel}$ from the measured leptonic decay rate using Eq. (43). The $\phi_{V}^{\perp}{}^{(2)}(x,\mu)$ term in the LCDA depends on the ratios $\langle v_{S}^{2}\rangle_{V}$, $\langle v_{S}^{4}\rangle_{V}$, $\langle v_{D}^{2}\rangle_{V}$, $\langle DE_{0}\rangle_{V}$, and $\langle DE_{1}\rangle_{V}$, see Eq. (25b), while the remaining contributions in Eq. (24) depend only on $\langle v_{S}^{4}\rangle_{V}$. We can eliminate the ratios $\langle B\rangle_{V}$ and $\langle DE_{0}\rangle_{V}$ in $f_{V}^{\perp}(\mu_{0})$ and $\phi_{V}^{\perp}{}^{(2)}(x,\mu)$ by using the Gremm-Kapustin relations in Eqs. (65). We obtain

[TABLE]

and

[TABLE]

When computing the convolution in Eq. (45), the terms $f_{V}^{\perp}(\mu_{0})\phi_{V}^{\perp}{}^{(n)}(x,\mu_{0})$ for $n=\alpha_{s}$, 2, 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 $\alpha_{s}v^{0}$ and $v^{4}$, so that

[TABLE]

We use the results in Eqs. (49) to compute the direct amplitude. Similarly, when calculating the convolution in Eq. (45), we ignore the order $\alpha_{s}$ correction to $T_{H}(x,\mu)$ for $n=\alpha_{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}^{\parallel}$ from Eq. (43) using the measured leptonic decay rate, and $\alpha(\mu_{0})=1/132$ regardless of the vector meson state $V$. We set the central values of $\mu_{0}$ and $\mu$ to be $m_{V}$ and $m_{H}$, respectively. We compute $\overline{m}(\mu)$, $\alpha_{s}(\mu_{0})$ and $\alpha_{s}(\mu)$ using RunDec Chetyrkin:2000yt ; Schmidt:2012az ; Herren:2017osy . We take the QED coupling constant at scale $\mu$ to be $\alpha(\mu)=1/128$. We resum the logarithms in $\mu/\mu_{0}$ to NLL accuracy using the Abel-Padé method Bodwin:2016edd ; Bodwin:2017wdu . 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}\langle v_{S}^{2}\rangle_{V}$ for $V=J/\psi$ and $\Upsilon(nS)$ have been obtained from potential-model (Cornell potential) calculations in Refs. Bodwin:2007fz ; Chung:2010vz :

[TABLE]

The first uncertainty comes from the variation of the potential-model parameters, and the second uncertainty is taken to be $\pm 0.3$ and $\pm 0.1$ times the central value for $V=J/\psi$ and $V=\Upsilon(nS)$ respectively, which comes from the uncalculated corrections of relative order $v^{2}$. The potential-model calculations in Refs. Bodwin:2007fz ; Chung:2010vz also let us compute the binding energies $E_{V}=m_{V}-2m$, which are given by

[TABLE]

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}\langle v_{S}^{2}\rangle_{V}$. We use these values of the binding energies to compute the heavy-quark mass through the relation $m=\frac{1}{2}(m_{V}-E_{V})$. The authors of Refs. Bodwin:2007fz ; Bodwin:2006dn also found the relation $\langle v_{S}^{4}\rangle_{V}=\left(\langle v_{S}^{2}\rangle_{V}\right)^{2}[1+O(v^{2})]$, 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 $\langle v_{S}^{4}\rangle_{V}$ to be $\left(\langle v_{S}^{2}\rangle_{V}\right)^{2}$ and take the uncertainty to be $\pm 0.3$ and $\pm 0.1$ times the central value of $\langle v_{S}^{4}\rangle_{V}$ for $V=J/\psi$ and $V=\Upsilon(nS)$, respectively. The ratios $\langle v_{D}^{2}\rangle_{V}$ and $\langle DE_{1}\rangle_{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 $\pm 0.09$ for $V=J/\psi$, and $\pm 0.01$ for $V=\Upsilon(nS)$, respectively.

We now list the sources of uncertainties in ${\cal A}_{\rm dir}$. We account for the uncertainties in the NRQCD LDMEs as discussed above. We vary the scales $\mu$ and $\mu_{0}$ between $\frac{1}{2}m_{H}<\mu<2m_{H}$ and $\frac{1}{2}m_{V}<\mu_{0}<2m_{V}$, while we ignore the negligibly small shifts in the QED couplings $\alpha(\mu_{0})$ and $\alpha(\mu)$ from scale variations. We also ignore the uncertainties from $m_{H}$ and $m_{V}$. We consider the uncertainty in $f_{V}^{\parallel}$ that originates from the experimental uncertainties in the leptonic decay rate. We add the uncertainties in quadrature. Our numerical results for ${\cal A}_{\rm dir}$ are shown in Table. 2. The imaginary part of ${\cal A}_{\rm dir}$ comes from the imaginary part in the order $\alpha_{s}$ correction to $T_{H}(x,\mu)$. Note that the uncertainties in the real and imaginary parts of ${\cal A}_{\rm dir}$ are correlated.

For $J/\psi$, the uncertainties from the LDMEs are comparable to the uncertainties from the variations of $\mu_{0}$ and $\mu$. If we ignore the uncertainties from scale variations, the uncertainty in ${\rm Re}[{\cal A}_{\rm dir}]$ is about 8% of the central value, and the uncertainty in ${\rm Im}[{\cal A}_{\rm dir}]$ is about 10% of the central value. This is comparable to the nominal size of the relative order $v^{4}$ correction.

For $\Upsilon(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 ${\cal A}_{\rm 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 $\langle v_{D}^{2}\rangle_{V}$ and $\langle DE_{1}\rangle_{V}$ are constrained, the reduction of the uncertainties from variations of scales would require calculation of the order $\alpha_{s}^{2}$ and order $\alpha_{s}v^{2}$ correction to the decay constant and the LCDA, and the order $\alpha_{s}^{2}$ correction to $T_{H}(x,\mu)$.

We again compare our results for ${\cal A}_{\rm dir}$ with the calculation in Ref. Bodwin:2017wdu . Our results for ${\cal A}_{\rm dir}$ in Table. 2 are compatible with those in Ref. Bodwin:2017wdu within uncertainties. For $J/\psi$, the uncertainty for ${\cal A}_{\rm dir}$ is smaller than in Ref. Bodwin:2017wdu , owing to the explicit calculation of the relative order $v^{4}$ corrections included in this work. On the other hand, for $\Upsilon(nS)$, the uncertainties for ${\cal A}_{\rm dir}$ are larger than those in Ref. Bodwin:2017wdu ; the main reason is that in Ref. Bodwin:2017wdu , the uncertainty from variations in the scales $\mu_{0}$ and $\mu$ were not taken into account, and instead, the uncertainties from the uncalculated order $\alpha_{s}^{2}$ and order $\alpha_{s}v^{2}$ corrections to the real part of ${\cal A}_{\rm dir}$ were estimated to be $C_{F}C_{A}\alpha_{s}^{2}(m)/\pi^{2}$ and $C_{F}\alpha_{s}(m)v^{2}/\pi$ of the central value, and the uncertainty from the uncalculated order $\alpha_{s}^{2}$ correction to the imaginary part was estimated to be $C_{A}\alpha_{s}(m)/\pi$ of the central value. We note that these estimates lead to smaller uncertainties compared to uncertainties estimated from variations of the scales $\mu_{0}$ and $\mu$. If we would use the same uncertainty estimates used in Ref. Bodwin:2017wdu for $\Upsilon(nS)$, then the uncertainties in ${\rm Re}[{\cal A}_{\rm dir}]$ would reduce by a factor of 2, while the uncertainties in ${\rm Im}[{\cal A}_{\rm dir}]$ would increase slightly.

IV.3 Decay rate

We now compute the total decay rate from the numerical results of ${\cal A}_{\rm ind}$ and ${\cal A}_{\rm dir}$ in Table. 2. Our results for the decay rates $\Gamma(H\to V+\gamma)$ are shown in Table 3. When computing the uncertainties in $\Gamma(H\to V+\gamma)$, we consider the correlation between the uncertainties in the real and imaginary parts of ${\cal A}_{\rm dir}$. We also consider the correlation between the uncertainties in ${\cal A}_{\rm ind}$ and ${\cal A}_{\rm dir}$ that comes from the measured leptonic decay rates $\Gamma(V\to\ell^{+}\ell^{-})$. The uncertainty from uncalculated corrections of relative order $m_{V}^{2}/m_{H}^{2}$ is taken to be $m_{V}^{2}/m_{H}^{2}$ of the central value of the decay rate. We also compute the branching ratios ${\rm Br}(H\to V+\gamma)$ by using the total decay rate of the Higgs $\Gamma_{H}$ computed in Refs. Dittmaier:2011ti ; Dittmaier:2012vm : for $m_{H}=125.18$ GeV, it is $\Gamma_{H}=4.10$ MeV, with uncertainties given by $+4.0$% and $-3.9$% of the central value. Our results for the branching ratios ${\rm Br}(H\to V+\gamma)$ are shown in Table 3.

For $V=J/\psi$, ${\cal A}_{\rm ind}$ is more than an order of magnitude larger than ${\cal A}_{\rm dir}$, and so, the uncertainty in the decay rate $\Gamma(H\to J/\psi+\gamma)$ is dominated by the uncertainty in ${\cal A}_{\rm ind}$. As a result, the uncertainty in the prediction for $\Gamma(H\to J/\psi+\gamma)$ is about 3% of the central value. On the other hand, for $V=\Upsilon(nS)$, ${\cal A}_{\rm ind}$ and ${\cal A}_{\rm dir}$ are comparable. Due to the large cancellation between ${\cal A}_{\rm ind}$ and ${\cal A}_{\rm dir}$ for $\Upsilon(nS)$, the uncertainty in $\Gamma[H\to\Upsilon(nS)+\gamma]$ is sensitive to the uncertainty in ${\cal A}_{\rm dir}$.

Our results are compatible with the previous calculation in Ref. Bodwin:2017wdu within errors. The uncertainties in $\Gamma(H\to J/\psi+\gamma)$ and ${\rm Br}(H\to J/\psi+\gamma)$ are slightly smaller than those of Ref. Bodwin:2017wdu , which is a result of the reduction of the uncertainty in ${\cal A}_{\rm ind}$ resulting from the improved measurement of $m_{H}$. On the other hand, the uncertainties in $\Gamma[H\to\Upsilon(nS)+\gamma]$ and ${\rm Br}[H\to\Upsilon(nS)+\gamma]$ are slightly larger than those of Ref. Bodwin:2017wdu . As we discussed in the previous section, if we would use the same estimates for the uncertainties used in Ref. Bodwin:2017wdu , then the uncertainties in ${\cal A}_{\rm dir}$ would be reduced, leading to uncertainties in $\Gamma[H\to\Upsilon(nS)+\gamma]$ and ${\rm Br}[H\to\Upsilon(nS)+\gamma]$ that are smaller than those of Ref. Bodwin:2017wdu .

We also consider the case where the Yukawa coupling $y_{Q}$ deviates from the Standard Model by a factor of $\kappa_{Q}$. In this case, the Higgs decay rate into $V+\gamma$ is given by $|\kappa_{Q}{\cal A}_{\rm dir}+{\cal A}_{\rm ind}|^{2}$. The decay rates $\Gamma[H\to J/\psi+\gamma]$ and $\Gamma[H\to\Upsilon(nS)+\gamma]$ for $-1.5<\kappa_{Q}<3.5$ are plotted in Fig. 2. For $J/\psi$, the decay rate shows moderate dependence on $\kappa_{c}$, while for $\Upsilon(nS)$, the decay rates are very sensitive to $\kappa_{b}$. In the case of $J/\psi$, a reduction of uncertainties in ${\cal A}_{\rm dir}$ will only have a small effect on the sensitivity of the $H\to J/\psi+\gamma$ rate on $\kappa_{c}$, because the uncertainty in $\Gamma(H\to J/\psi+\gamma)$ is dominated by the uncertainty in ${\cal A}_{\rm ind}$. On the other hand, for $\Upsilon$, the decay rate $\Gamma[H\to\Upsilon(nS)+\gamma]$ will be even more sensitive to $\kappa_{b}$ if the accuracy in ${\cal A}_{\rm dir}$ is improved. The reason for this sensitivity is the large cancellation between ${\cal A}_{\rm dir}$ and ${\cal A}_{\rm ind}$ that only happens for $\kappa_{b}$ close to 1. We note that even though the Standard-Model values of $\Gamma[H\to\Upsilon(nS)+\gamma]$ are very small, the decay rate $\Gamma[H\to\Upsilon(nS)+\gamma]$ can still probe large deviations of the Higgs bottom coupling from its Standard Model value. Indeed $\Gamma[H\to\Upsilon(nS)+\gamma]$ may be two/three orders of magnitude larger than the Standard Model value for $\kappa_{b}\lesssim-1$ or $\kappa_{b}\gtrsim 3$. In particular, the case $\kappa_{b}\lesssim-1$ corresponds to a Higgs bottom coupling whose sign is opposite to the Standard Model Yukawa coupling.

V Summary and Discussion

In this work we computed the order $v^{4}$ correction to the decay rate $\Gamma(H\to V+\gamma)$ where $V=J/\psi$ or $\Upsilon(nS)$ based on the nonrelativistic QCD (NRQCD) factorization formalism. By using the light-cone approach, we resummed the logarithms of $m_{H}^{2}/m_{V}^{2}$ that appear in higher order corrections in $\alpha_{s}$ to all orders in $\alpha_{s}$ at the next-to-leading logarithmic accuracy.

If we consider that $\alpha_{s}\approx v^{2}$ at the scale of the heavy-quark mass, the corrections of order $\alpha_{s}^{2}$ and $\alpha_{s}v^{2}$ would be of the same order as the order $v^{4}$ corrections computed in this paper. The calculation of order $\alpha_{s}^{2}$ and order $\alpha_{s}v^{2}$ corrections to the decay rate $\Gamma(H\to V+\gamma)$ requires calculation of the two-loop correction to $c_{0}$ and the one-loop correction to $c_{\bm{D}^{2}}$ in Eq. (1), respectively. In the light-cone approach, the order $\alpha_{s}^{2}$ correction to the hard part $T_{H}(x,\mu)$ in Eq. (16), as well as the two-loop correction to $\tilde{c}_{0}(x)$ and the one-loop correction to $\tilde{c}_{\bm{D}^{2}}(x)$ are necessary. Since these corrections have not been computed, we include the effects of these higher-order corrections in the uncertainties.

Our numerical results for the Standard Model values of the decay rates $\Gamma(H\to V+\gamma)$ and branching ratios ${\rm Br}(H\to V+\gamma)$ are shown in Table 3. The corrections computed in this work improve the theoretical accuracy in the prediction of the decay rate $\Gamma(H\to J/\psi+\gamma)$ compared to a previous calculation in Ref. Bodwin:2017wdu . If we would have used the same method to estimate the uncertainties as in Ref. Bodwin:2017wdu , we would have found uncertainties in $\Gamma[H\to\Upsilon(nS)+\gamma]$ reduced by almost a factor of two compared to our results in Table 3, so that they would become comparable to the results in Ref. Bodwin:2017wdu , as the error of $\Gamma[H\to\Upsilon(nS)+\gamma]$ would be dominated by the scale uncertainties. However, due to the difference in the method to estimate uncertainties, our results for $\Gamma[H\to\Upsilon(nS)+\gamma]$ have larger uncertainties than those of Ref. Bodwin:2017wdu .

We note that the decay rates $\Gamma(H\to V+\gamma)$ depend on the LDMEs $\langle V|\psi^{\dagger}\epsilon^{i}(\lambda)\sigma^{j}(-\tfrac{i}{2})^{2}$ $\times\overleftrightarrow{\bm{D}}^{(i}\overleftrightarrow{\bm{D}}^{j)}\chi|0\rangle$ and $\langle V|\psi^{\dagger}\bm{\epsilon}(\lambda)\cdot\tfrac{1}{2}[\bm{\sigma}\times(\overleftrightarrow{\bm{D}}\times g_{s}\bm{E}-g_{s}\bm{E}\times\overleftrightarrow{\bm{D}})]\chi|0\rangle$ that are currently unknown. In our numerical results, we have estimated their sizes according to the conservative power counting in Refs. Brambilla:2001xy ; Brambilla:2002nu ; Brambilla:2006ph ; Brambilla:2008zg , and included their effects in the uncertainties. For $V=J/\psi$, the uncertainties from these unknown LDMEs are significant compared to the total uncertainty in $\Gamma(H\to J/\psi+\gamma)$. Therefore, knowledge of these LDMEs will improve the accuracy in the prediction of the decay rate $\Gamma(H\to J/\psi+\gamma)$. Also the systematic inclusion of relative order $v^{2}$ effects in the determination of the matrix element $\langle V|\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)(-\tfrac{i}{2}\overleftrightarrow{\bm{D}})^{2}\chi|0\rangle$ would significantly improve the determination of the $H\to J/\psi+\gamma$ decay rate. A calculation of the LDMEs in potential NRQCD may help constrain these LDMEs Brambilla:2002nu .

In Ref. Bodwin:2013gca , the authors estimated that the process $H\to J/\psi+\gamma$ could be measured at the HL-LHC experiment through the leptonic decays of $J/\psi$ into $e^{+}e^{-}$ and $\mu^{+}\mu^{-}$. 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 $\Gamma(H\to J/\psi+\gamma)$ at 95% confidence level that is about 15 times larger the Standard Model value ATL-PHYS-PUB-2015-043 . Although the prospect for measuring the process $H\to J/\psi+\gamma$ 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 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 $\Upsilon(nS)$, the Standard Model values for the decay rates $\Gamma[H\to\Upsilon(nS)+\gamma]$ are about three orders of magnitude smaller than $\Gamma(H\to J/\psi+\gamma)$, owing to large cancellations between the direct and indirect amplitudes.

We show the decay rates $\Gamma(H\to J/\psi+\gamma)$ and $\Gamma[H\to\Upsilon(nS)+\gamma]$ when the Higgs charm coupling is rescaled by a factor $\kappa_{c}$ and the Higgs bottom coupling is rescaled by a factor $\kappa_{b}$ compared to the Standard Model in Fig. 2. Due to the cancellation between direct and indirect amplitudes that occurs when $\kappa_{b}\approx 1$, the decay rates $\Gamma[H\to\Upsilon(nS)+\gamma]$ are highly sensitive to $\kappa_{b}$, so much so that the combined rate of $\Gamma[H\to\Upsilon(nS)+\gamma]$ for $n=1,2$, and $3$ increases to be larger than one half of the Standard Model value of $\Gamma(H\to J/\psi+\gamma)$ for $\kappa_{b}\lesssim-1$ or $\kappa_{b}\gtrsim 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.

Acknowledgements.

We thank Geoffrey Bodwin for helpful discussions. The research of N. B. is supported by the DFG Grant No. BR 4058/2-2. N. B. and A. V. acknowledge support from the DFG cluster of excellence “Origins” (www.origins-cluster.de). The work of H. S. C. is supported by the Alexander von Humboldt Foundation. The work of V. S. is supported in part by the National Science Foundation of China (No. 11135006, No. 11275168, No. 11422544, No. 11375151, No. 11535002) and the Zhejiang University Fundamental Research Funds for the Central Universities (2017QNA3007).

Appendix A Kinematics

In this section, we present the kinematical conventions that we use for the perturbative $Q\bar{Q}$ and the $Q\bar{Q}g$ states. The kinematical conventions that we use here are identical to the ones used in Ref. Brambilla:2017kgw .

A.1 Two-body kinematics

We let the $Q$ and the $\bar{Q}$ to have the momenta $p_{1}$ and $p_{2}$, respectively. We denote the total momentum of the $Q\bar{Q}$ system as $P=p_{1}+p_{2}$, and the relative momentum of the $Q$ and the $\bar{Q}$ as $q=\frac{1}{2}(p_{1}-p_{2})$. In the rest frame of the $Q\bar{Q}$, $q$ and $P$ are given by $q=(0,\bm{q})$, and $P=(2\sqrt{m^{2}+\bm{q}^{2}},\bm{0})$, which leads to $p_{1}=(\sqrt{m^{2}+\bm{q}^{2}},\bm{q})$ and $p_{2}=(\sqrt{m^{2}+\bm{q}^{2}},-\bm{q})$.

A.2 Three-body kinematics

We again let the $Q$ and the $\bar{Q}$ to have the momenta $p_{1}$ and $p_{2}$, respectively, and the gluon carry momentum $k_{g}$. We set the total momentum of the $Q\bar{Q}g$ system to be $P=p_{1}+p_{2}+k_{g}$. We define $q_{1}=\frac{1}{2}(p_{1}-p_{2})$ and $q_{2}=\frac{1}{6}(2k_{g}-p_{1}-p_{2})$, so that $p_{1}=\frac{1}{3}P+q_{1}-q_{2}$, $p_{2}=\frac{1}{3}P-q_{1}-q_{2}$, and $k_{g}=\frac{1}{3}P+2q_{2}$. In the rest frame of the $Q\bar{Q}g$ system, where $\bm{P}=\bm{p}_{1}+\bm{p}_{2}+\bm{k}_{g}=0$, $p_{1}=(\sqrt{m^{2}+(\bm{q}_{1}-\bm{q}_{2})^{2}},\bm{q}_{1}-\bm{q}_{2})$, $p_{2}=(\sqrt{m^{2}+(\bm{q}_{1}+\bm{q}_{2})^{2}},-\bm{q}_{1}-\bm{q}_{2})$, and $k_{g}=(2|\bm{q}_{2}|,2\bm{q}_{2})$, so that $P^{0}=2|\bm{q}_{2}|+\sqrt{(\bm{q}_{1}+\bm{q}_{2})^{2}+m^{2}}+\sqrt{(\bm{q}_{1}-\bm{q}_{2})^{2}+m^{2}}$.

Appendix B Nonrelativistic expansion of Dirac spinors

In order to obtain the matching conditions, we need to express the QCD amplitudes in terms of the 3-momenta of the particles and the 2-component Pauli spinors in the frame where the NRQCD LDMEs are defined. One way to accomplish this is to compute the QCD amplitudes using explicit forms of the Dirac spinors $u(p_{1},s)$ and $v(p_{2},s^{\prime})$ with definite spin. In this appendix, we introduce a simple way to compute the QCD amplitudes with the explicit Dirac spinors that can be easily used in automated calculations.

For the perturbative $Q\bar{Q}$ or the $Q\bar{Q}g$ states, the frame where the NRQCD LDMEs are defined is the frame where the total momentum $P$ of the state is at rest ($\bm{P}=\bm{0}$), so that $P^{0}=\sqrt{P^{2}}$ and $A^{0}=A\cdot P/\sqrt{P^{2}}$ for an arbitrary 4-vector $A$. In this frame, the nonrelativistically normalized Dirac spinors in Dirac basis are given by Bodwin:1994jh

[TABLE]

where $m_{i}^{2}=p_{i}^{2}$, $E_{i}=p_{i}\cdot P/\sqrt{P^{2}}$, and $N_{i}=1/\sqrt{2E_{i}(E_{i}+m_{i})}$. In general, $m_{1}$ and $m_{2}$ can be different. A QCD amplitude for production of a quark with momentum $p_{1}$ and an antiquark with momentum $p_{2}$ involves

[TABLE]

For example, $\bar{u}(p_{1},s)\Gamma v(p_{2},s^{\prime})={\rm tr}[\Gamma v(p_{2},s^{\prime})\otimes\bar{u}(p_{1},s)]$, where $\Gamma$ is a product of gamma matrices. Since $\eta_{s^{\prime}}\otimes\xi_{s}^{\dagger}$ is a $2\times 2$ matrix, it can be written as a linear combination of a $2\times 2$ identity matrix $I$ and the $\bm{\sigma}$ matrices:

[TABLE]

We note that in the Dirac basis,

[TABLE]

and

[TABLE]

so that

[TABLE]

By using Eq. (58), we can compute the quantities of the form $\bar{u}(p_{1},s)\Gamma v(p_{2},s^{\prime})={\rm tr}[\Gamma v(p_{2},s^{\prime})\otimes\bar{u}(p_{1},s)]$, where $\Gamma$ is any product of gamma matrices, as traces of gamma matrices. This can be easily implemented in automated calculations using FeynCalc. We obtain the nonrelativistic expansion of a QCD amplitude by computing the amplitude using Eq. (58) and expanding in powers of the small 3-momenta of the $Q\bar{Q}$ and the $Q\bar{Q}g$ states. Using this method, we easily reproduce the explicit expressions for $\bar{u}(p_{1},s)\Gamma v(p_{2},s^{\prime})$ for $\Gamma=1$, $\gamma^{\mu}$, $\gamma^{\mu}\gamma^{\nu}-\gamma^{\nu}\gamma^{\mu}$, and $\gamma^{\mu}\gamma^{\nu}\gamma^{\sigma}-\gamma^{\sigma}\gamma^{\nu}\gamma^{\mu}$ found in Ref. Braaten:1996rp . We use this method to compute the QCD amplitudes in Eqs. (7), (10), (19), and (21).

The expression in Eq. (58) may serve to relate the nonrelativistic expansion method that we use in this work with the covariant spin-projector method used in previous calculations of the $H\to J/\psi+\gamma$ process in Refs. Bodwin:2013gca ; Bodwin:2014bpa . In the covariant spin-projector method, spin-singlet and spin-triplet contributions are computed separately. The standard forms of the spin projectors, such as the ones used in Ref. Bodwin:2002hg , can be obtained from Eq. (58) by projecting to a spin-singlet or a spin-triplet state using the Clebsch-Gordan coefficients. For a spin-singlet state, $\sum_{s,s^{\prime}}\langle\tfrac{1}{2}s,\tfrac{1}{2}s^{\prime}|00\rangle\xi_{s}^{\dagger}\bm{\sigma}\eta_{s^{\prime}}/\sqrt{2}$ vanishes, and $\sum_{s,s^{\prime}}\langle\tfrac{1}{2}s,\tfrac{1}{2}s^{\prime}|00\rangle\xi_{s}^{\dagger}\eta_{s^{\prime}}/\sqrt{2}$ is, up to a phase, equal to $1$. Hence, the spin-singlet projector is, up to a phase, given by the contribution in Eq. (58) that is proportional to $\xi^{\dagger}_{s}\eta_{s^{\prime}}$, with $\xi^{\dagger}_{s}\eta_{s^{\prime}}$ replaced by $1$. Similarly, for a spin-triplet state with polarization $\lambda$, $\sum_{s,s^{\prime}}\langle\tfrac{1}{2}s,\tfrac{1}{2}s^{\prime}|1\lambda\rangle\xi_{s}^{\dagger}\eta_{s^{\prime}}/\sqrt{2}$ vanishes for all $\lambda\in\{-1,0,+1\}$, and $\sum_{s,s^{\prime}}\langle\tfrac{1}{2}s,\tfrac{1}{2}s^{\prime}|1\lambda\rangle\xi_{s}^{\dagger}\bm{\sigma}\eta_{s^{\prime}}/\sqrt{2}$ is, up to a phase, equal to the polarization vector of the spin-triplet state. Therefore, the spin-triplet projector is, up to a phase, given by the contribution in Eq. (58) that is proportional to $\xi^{\dagger}_{s}\bm{\sigma}\eta_{s^{\prime}}$, with $\xi^{\dagger}_{s}\bm{\sigma}\eta_{s^{\prime}}$ replaced by the polarization vector of the spin-triplet state. As a result, the spin-singlet (spin-triplet) contribution of a QCD amplitude computed using the covariant spin-projector method is, up to a phase, equivalent to the contribution proportional to $\xi^{\dagger}_{s}\eta_{s^{\prime}}$ ($\xi^{\dagger}_{s}\bm{\sigma}\eta_{s^{\prime}}$) in the amplitude computed in the nonrelativistic expansion method. The phase conventions for the spin-singlet and the spin-triplet projectors depend on the conventions for the Clebsch–Gordan coefficients and the Pauli spinors $\xi$ and $\eta$.

One advantage of using the covariant spin-projector method is that covariant expressions of the QCD amplitudes can be obtained easily, unlike the nonrelativistic expansion method in Ref. Braaten:1996rp . On the other hand, covariant expressions can also be obtained if we use Eq. (58) to carry out the nonrelativistic expansion of spinors. Also, in this work, there is no advantage in computing the spin-singlet and spin-triplet contributions separately, because both contributions appear in the matching condition in Eq. (6) simultaneously. Therefore, we compute the QCD amplitudes in Eqs. (7), (10), (19), and (21) using Eq. (58).

Appendix C Short-distance coefficients for $\bm{H\to h_{c}+\gamma}$

The calculation of the short-distance coefficients presented in Sec. II.1 can be easily applied to production amplitudes of other quarkonium states. For example, projecting onto the $J^{PC}=1^{+-}$ state gives us the short-distance coefficients for Higgs decay into $h_{c}+\gamma$ to relative order $v^{2}$. To achieve this, we first write down the matching condition for the $J^{PC}=1^{+-}$ case as

[TABLE]

and

[TABLE]

where we define [see also Eq. (2)]

[TABLE]

The color-octet matrix elements for the $J^{PC}=1^{+-}$ case are obtained from the color-octet matrix elements for the $J^{PC}=1^{--}$ case in Eq. (6) by making the replacements $g_{s}\bm{E}\to g_{s}\bm{B}$ and $g_{s}\bm{B}\to g_{s}\bm{E}$.

From the $J^{PC}=1^{+-}$ contributions to the $Q\bar{Q}$ and $Q\bar{Q}g$ amplitudes in Eqs. (7) and (10), we obtain the short-distance coefficients

[TABLE]

where $r=\frac{4m^{2}}{m_{H}^{2}}$.

Recently, a computation of the decay rate $\Gamma(H\to h_{c}+\gamma)$ in the NRQCD factorization formalism at leading order in $v$ appeard in Ref. Mao:2019hgg . This is equivalent to our calculation of the short-distance coefficient $c_{1}$ in Eq. (62a), which leads to the expression for the decay rate $\Gamma(H\to h_{c}+\gamma)$ at leading order in $v$ that is given by

[TABLE]

where the sum is over the polarization of the $h_{c}$ and $\Phi$ is the phase-space and normalization factor given in Eq. (42). Our result in Eq. (63) agrees with the decay rate computed in Ref. Mao:2019hgg .

Appendix D Gremm-Kapustin relations

The Gremm-Kapustin relations Gremm:1997dq are obtained from

[TABLE]

where ${\cal O}$ is an NRQCD operator, and $H$ is the NRQCD Hamiltonian. Computing the commutator $[{\cal O},H]$ leads to the following relations:

[TABLE]

where, in calculating the commutator $[{\cal O},H]$ we included operators in the Hamiltonian up to $1/m^{2}$ accuracy, with the Wilson coefficients at order $\alpha_{s}^{0}$, and kept NRQCD operators up to dimension 7. Hence, the relations in Eqs. (65) are valid up to order $v^{4}$ relative to the leading-order LDME $\langle V|\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle$ and at leading order in $\alpha_{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. Brambilla:2001xy ; Brambilla:2002nu ; Brambilla:2006ph ; Brambilla:2008zg . Since the binding energy $m_{V}-2m$ scales like $mv^{2}$, the left-hand side of Eq. (65) is suppressed by $v^{6}$ compared to the leading-order LDME $\langle V|\psi^{\dagger}\bm{\sigma}\cdot\bm{\epsilon}(\lambda)\chi|0\rangle$. Therefore, the LDME $\langle V|\psi^{\dagger}\epsilon^{i}(\lambda)\sigma^{j}(\overleftrightarrow{\bm{D}}^{(i}g_{s}\bm{E}^{j)}+g_{s}\bm{E}^{(i}\overleftrightarrow{\bm{D}}^{j)})\chi|0\rangle$ 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 leading-order LDME, and scales like $v^{15/2}$. Similarly, the left-hand side of Eq. (65) is suppressed by $v^{5}$ compared to the leading-order LDME, and hence, the LDME on the right-hand side of Eq. (65) scales like $v^{13/2}$ and does not contribute to the amplitude in Eq. (1) at relative order $v^{4}$ accuracy.