$Z$-boson decays to a vector quarkonium plus a photon

We compute the decay rates for the processes $Z\to V+\gamma$, where $Z$ is the $Z$ boson, $\gamma$ is the photon, and $V$ is one of the vector quarkonia $J/\psi$ or $\Upsilon(nS)$, with $n=1$, $2$, or $3$. Our computations include corrections through relative orders $\alpha_s$ and $v^2$ and resummations of logarithms of $m_Z^2/m_Q^2$, to all orders in $\alpha_s$, at NLL accuracy. ($v$ is the velocity of the heavy quark $Q$ or the heavy antiquark $\bar{Q}$ in the quarkonium rest frame, and $m_Z$ and $m_Q$ are the masses of $Z$ and $Q$, respectively.) Our calculations are the first to include both the order-$\alpha_s$ correction to the light-cone distributions amplitude and the resummation of logarithms of $m_Z^2/m_Q^2$ and are the first calculations for the $\Upsilon(2S)$ and $\Upsilon(3S)$ final states. The resummations of logarithms of $m_Z^2/m_Q^2$ that are associated with the order-$\alpha_s$ and order-$v^2$ corrections are carried out by making use of the Abel-Pad\'e method. We confirm the analytic result for the order-$v^2$ correction that was presented in a previous publication, and we correct the relative sign of the direct and indirect amplitudes and some choices of scales in that publication. Our branching fractions for $Z\to J/\psi+\gamma$ and $Z\to \Upsilon(1S)+\gamma$ differ by $2.0\,\sigma$ and $-4.0\,\sigma$, respectively, from the branching fractions that are given in the most recent publication on this topic (in units of the uncertainties that are given in that publication). However, we argue that the uncertainties in the rates are underestimated in that publication.


I. INTRODUCTION
The rare decays of the Higgs boson (H) to a vector quarkonium (V ) and a photon (γ) have been proposed as processes with which to measure the Higgs-boson couplings to the charm and bottom quarks [1]. Even at a high-luminosity LHC, observations of these decay processes would be challenging. It has been pointed out in Refs. [2,3] that the decays of the Z boson Z → V + γ could provide means to calibrate the experimental techniques that might be used to measure the H → V + γ decay rates.
As was emphasized in Ref. [1], in the decays H → V + γ, two processes give important contributions to the amplitude: (1) a direct process, in which the Higgs boson decays to a heavy-quark-antiquark pair (QQ), which emits a photon and evolves into a quarkonium; (2) an indirect process, in which the Higgs boson decays through a virtual heavy-quark or W -boson loop into a photon and a virtual photon, with the virtual photon decaying into a heavy quarkonium. In the case of the decays H → V +γ, the indirect process is enhanced for massive particles in the virtual loop because the Higgs-boson coupling to the loop particle is proportional to the mass of the particle.
In a classic paper [4], analytic expressions were given for the direct amplitudes and the corresponding decay rates for Z-boson decays to a photon plus an S-wave or a P -wave quarkonium. These expressions were calculated at leading order (LO) in α s , the QCD running coupling, and at LO in v 2 , where v is the velocity of the heavy quark (Q) or the heavy antiquark (Q) in the heavy-quarkonium rest frame.
Calculations of exclusive quarkonium production processes can be simplified by making use of the light-cone approach [5,6], which yields a systematic expansion in powers of m V /m hard , where m V is the quarkonium mass and m hard is the hard-scattering scale, which is of order the Z-boson mass m Z in the present case. In the light-cone approach, nonperturbative effects in the quarkonium system are parametrized in terms of the quarkonium light-cone distribution amplitudes (LCDAs). A heavy-quarkonium LCDA can, by virtue of nonrelativistic QCD (NRQCD) factorization [7], be written as a sum of products of shortdistance coefficients times NRQCD long-distance matrix elements (LDMEs) [8].
In Ref. [9], calculations of the rates for Z-boson decays to a photon plus η c , J/ψ, χ c0 , χ c1 , χ c2 , or h c were presented. These calculations were based on the direct amplitude and were carried out at LO in α s and v 2 in both the NRQCD and light-cone formalisms.
In Ref. [2], the decay rates for the processes Z → V + γ, where V is the J/ψ or the Υ(1S), were computed in the leading-power light-cone approximation, which is valid up to corrections of order m 2 V /m 2 Z . The calculations were carried out at next-to-leading order (NLO) in α s and v 2 . Reference [2] also gave the first result for the short-distance coefficient of the order-v 2 (relativistic) corrections. The calculations in Ref. [2] included contributions from the indirect production process. These contributions were found to be small, producing effects of less than 1 % on the rates, because, in contrast with the Higgs-boson indirect amplitude, the Z-boson indirect amplitude is not proportional to the mass of the loop particle. The calculation in Ref. [2] did not include the effects of resummation of large logarithms of the ratio m 2 Z /m 2 Q , where m Q is the heavy-quark mass. This resummation was estimated in Ref. [2] to produce a 1.5 % effect in the rate for Z → J/ψ + γ.
In Ref. [3], the decay rates for the processes Z → V +γ, where V is the J/ψ or the Υ(1S), were also computed in the leading-power light-cone approximation at NLO in v 2 and at LO in α s . Logarithms of m 2 Z /m 2 Q were resummed to all orders in α s at leading logarithmic (LL) accuracy. In the case of the order-v 2 correction, the resummation of logarithms of m 2 Z /m 2 Q was carried out by introducing a model for the LCDA whose second moment in the light-cone momentum fraction x (in the narrow-width approximation) matches the second x moment of the order-v 2 term in the nonrelativistic expansion of the LCDA. It was found in Ref. [3] that the resummation effects are much larger than the 1.5 % estimate that was given in Ref. [2].
In principle, one can carry out the resummation of logarithms of m 2 Z /m 2 Q for the order-v 2 and order-α s corrections to the LCDA directly, avoiding the unknown uncertainties that are associated with the introduction of a model LCDA. However, as was pointed out in Refs. [10,11], the standard approach for such calculations, namely, expansion in a series in the LO evolution eigenvectors (Gegenbauer polynomials), fails because the eigenvector series diverges, even though the evolved LCDA itself is well defined. This divergence can be traced to the fact that the order-v 2 and order-α s corrections to the LCDA contain distributions (generalized functions) [11]. A general solution to this problem was given in Ref. [11], where it was shown that the evolved LCDA can be obtained by using the so-called Abel-Padé method to sum the divergent eigenvector series. The Abel-Padé method allows one to compute the resummation of logarithms of m 2 Z /m 2 Q for the order-v 2 and order-α s corrections to the LCDA from first principles.
In the present paper, we compute the decay rates for the processes Z → V + γ, where V is one of the four states J/ψ and Υ(nS), with n = 1, 2, or 3. Our computation is carried out at leading power in the light-cone formalism and through orders α s and v 2 .
Logarithms of m 2 Z /m 2 Q are resummed in the direct amplitude at next-to-leading-logarithmic (NLL) accuracy. The computations of the rates for Z → V + γ in this paper are the first to include both the order-α s corrections to the LCDA and the resummation of logarithms of m 2 Z /m 2 Q . The calculation includes the effects of the indirect process, as well as the effects of the direct process.
In comparison with the central values in Ref. [3], our branching fraction for Z → J/ψ + γ is shifted by about +12 %, which is +2.0 σ in the uncertainties of Ref. [3], and our branching fraction for Z → Υ(1S) + γ is shifted by about −11 %, which is −4.0 σ in the uncertainties of Ref. [3]. We argue that the uncertainties in the rates are underestimated in Ref. [3].
We have also confirmed the result in Ref. [2] for the short-distance coefficient of the orderv 2 correction. Our result for the relative sign between the direct and indirect amplitudes differs from that in Ref. [2], resulting in positive (negative) interference for the J/ψ + γ [Υ(nS) + γ] final state. As the indirect amplitude is small relative to the direct amplitude, the effect of this sign change is much less than the uncertainties in the calculation. We have also corrected some choices of scales in the calculation in Ref. [2]. The effects of these corrections tend to cancel the effects of the resummations of logarithms of m 2 Z /m 2 Q , which are not included in Ref. [2].
The remainder of this paper is organized as follows. In Sec. II we give the expression for the direct amplitude, and in Sec. III we discuss the resummation of logarithms of m 2 Z /m 2 Q in the direct amplitude. In Sec. IV we give the expression for the indirect amplitude. Section V contains a discussion of the numerical calculation of the rates and the uncertainties in that calculation. In Sec. VI, we present our numerical results and compare them with results from previous computations. Finally, in Sec. VII, we summarize and discuss our results.

II. LIGHT-CONE AMPLITUDE FOR THE DIRECT PROCESS
The light-cone amplitude for the direct process of Z → V + γ is given, up to corrections of relative order m 2 V /m 2 Z , by where Here, e(> 0) is the electric charge at momentum scale zero, e Q is the fractional charge of the heavy quark Q, f V is the decay constant of the longitudinally polarized vector quarkonium V , ǫ Z is the Z-boson polarization, ǫ V is the quarkonium polarization, ǫ γ and p γ are the photon polarization and momentum, respectively, µ is the renormalization scale, x is the Q momentum fraction of V , which runs from 0 to 1, and g Z and g A are defined by Here, G F is the Fermi constant, and (T f 3 ) L is the eigenvalue of the weak isospin of the lefthanded fermion f , whose value is +1/2 for f = u, c, t, ν e , ν µ , ν τ , and −1/2 for f = d, s, b, e, µ, τ . We use the convention ǫ 0123 = −1.
The longitudinally polarized LCDA φ V is defined in Refs. [12,13] as where p is the quarkonium momentum, z lies along the plus light-cone direction, and is a gauge link that makes the nonlocal operator gauge invariant. Here, g s = √ 4πα s , A µ a is the gluon field with the color index a = 1, 2, . . . , N 2 c − 1, N c = 3, T a is a generator of color SU(3) in the fundamental representation, and the symbol P denotes path ordering.
Note that we have included a factor (−1) in the definition (3) relative to the definition in Refs. [12,13] in order to obtain a positive value for the decay constant. We note that the definition (3) is equivalent to the definition of φ V in Ref. [14], from which we take the order-α s corrections to φ V .
Setting z to 0 and imposing the normalization condition we obtain which allows one to relate the decay constant f V to the quarkonium electromagnetic decay width Γ(V → e + e − ): Here, α(m V ) is the running electromagnetic coupling at the scale m V .
We expand the LCDA at the low-energy scale µ 0 , which is of order m Q , through order v 2 and through order α s : The quantity v 2 V is proportional to the ratio of the NRQCD LDME of order v 2 to the NRQCD LDME of order v 0 . The general expression for the ratio of the NRQCD LDME of order v 2k (k a nonnegative integer) to the NRQCD LDME of order v 0 is Here, ψ is the two-component (Pauli) spinor field that annihilates a heavy quark, χ † is the two-component spinor field that annihilates a heavy antiquark, σ is a Pauli matrix, |V (ǫ V ) denotes the vector quarkonium state in the quarkonium rest frame with spatial polarization ǫ V , and m Q denotes the quark pole mass. The light-cone functions on the right side of Eq. (8) are given by Here, the + and ++ functions are defined by The order-α s light-cone function φ V (x, µ 0 ) was computed in Ref. [14]. In Eq. (10c), we have replaced the pole mass m Q with m Q , the modified-minimal-subtraction (MS) mass, since the pole mass is ill defined, owing to renormalon ambiguities. This change affects the expression Ref. [2]. It can be inferred from the computation for the quarkonium transverse LCDA in Ref. [10] by using the fact that the relativistic corrections to the LCDA are independent of the quarkonium spin [15]. It can also be inferred from the calculation in Ref. [16] for S-wave B c mesons in the limit m c = m b . We have verified this result by using the NRQCD formalism to compute the complete order-v 2 contribution to the direct amplitude, which includes the order-v 2 contribution to φ V in Eq. (10b) and the order-v 2 contribution to the decay constant f V , and making use of the known order-v 2 contribution to f V [see Eq. (12) below].
The decay constant f V is given by where C F = (N 2 c − 1)/(2N c ) and C A = N c = 3 for color SU (3). We note that f V , as defined in Eq. (7), is scale invariant. Hence, the dependence of the expression in brackets on the right side of Eq. (12) on the scale µ 0 implies that Ψ V (0) depends on µ 0 in such a way as to render the complete expression scale invariant. The quarkonium wave function at the origin Ψ V (0) is related to the LO NRQCD LDME [7]: The LO and order-α s contributions to f V were computed in Ref. [14]. The order-v 2 contribution to f V was computed in Ref. [2]. It can be inferred from the order-v 2 contribution to the quarkonium electromagnetic decay rate in Ref. [7].
In this paper, we will use Eq. (7) to compute f V from the measured electromagnetic decay widths. As was pointed out in Ref. [3], this procedure eliminates the uncertainties in the calculation that arise from the use of Eq. (12) in conjunction with phenomenological determinations of Ψ V (0) and v 2 V . Equation (12) was used in the calculation in Ref. [2]. We defer a discussion of the impact of that choice to Sec. VI.
The hard-scattering kernel T H (x, µ) for the process Z → V + γ, through order α s , is [14] T where

A. The Gegenbauer expansion of the amplitude
The scale evolution of the LCDA can be computed conveniently by expanding the LCDA in Gegenbauer polynomials, which are the eigenvectors of the LO evolution kernel. The Gegenbauer expansion of the LCDA is where φ n is the nth Gegenbauer moment of φ V , which can be found by making use of the orthogonality property of the Gegenbauer polynomials: where N n = 4(2n + 3) (n + 1)(n + 2) .
Note that φ n (µ) vanishes for odd n because φ V (x, µ) is symmetric about x = 1/2. We define the LO, order-v 2 , and order-α s Gegenbauer moments of φ V as follows: The moments φ n (µ) can be written in terms of the moments φ n (µ 0 ) and an evolution matrix U nk (µ, µ 0 ): The LL and NLL expressions for U nk (µ, µ 0 ) can be found in the Appendix.
The Gegenbauer expansion of the hard-scattering kernel is given by where T n is the nth Gegenbauer moment of T H , which can be found by making use of the orthogonality property of the Gegenbauer polynomials: We define the LO and order-α s contributions to T n as follows: Making use of the orthogonality property of the Gegenbauer polynomials again, we can write the light-cone amplitude as We also define a decomposition of M into the LO, order-v 2 , and order-α s contributions: where By choosing the scale µ in M(µ) to be m Z , we guarantee that T n (µ) contains no large logarithms. We choose the initial scale of the LCDAs to be µ 0 = m Q . Then, logarithms of m 2 Z /m 2 Q are resummed by the evolution of φ n from the scale µ 0 = m Q to the scale µ = m Z . Using Eq. (24), we find that the resummed direct amplitude is given by We use this expression in our numerical calculations. We make use of expressions for the evolution through NLL accuracy, which are given in the Appendix.

B. The Abel-Padé method
The sum over n in Eq. (25) diverges for M (0,v 2 ) and M (0,1) [10,11]. As was explained in Ref. [11], such divergences can arise because the light-cone distributions contain generalized functions (distributions), rather than ordinary functions. In Ref. [11], it was shown that one can define the generalized functions as a limit of ordinary functions, which leads one to compute M (i,j) as follows: The expression in Eq. (27) is the Abel summation of the eigenfunction series for φ V (x, µ 0 ).
In Ref. [11], the Abel summation was erroneously applied to φ V (x, µ). (See Ref. [17].) We have corrected that error here. The correction amounts to the replacement of z m with z n in Eq. (27).
One can improve upon the convergence of the series in Eq. (27) in the limit z → 1, by constructing a Padé approximant for the nth partial sum before taking the limit z → 1. The use of the Padé approximant is effective in improving the convergence of the series because it provides an approximate analytic continuation for the function of z that is represented by the series. That analytic continuation is valid beyond the radius of convergence of the series, which is typically |z| = 1. The Abel-Padé method was tested extensively against known analytic results for M (i,j) in Ref. [11], and it converged rapidly to the correct value in all cases. We will use it throughout this paper to evaluate M (i,j) .

IV. AMPLITUDE FOR THE INDIRECT PROCESS
The amplitude for the indirect decay amplitude contains the axial-vector-vector triangle diagram as a subdiagram. The amplitude for the axial-vector-vector triangle diagram is given in Ref. [18]. In that paper, the conventions for γ 5 and the completely antisymmetric tensor ǫ ξµνρ are not specified. We fix the overall sign of the triangle amplitude in Ref. [18] in our conventions by requiring that it give the correct axial-vector anomaly. Then, we find that the indirect amplitude for the decay of a Z boson to a photon plus a virtual photon is given by where f denotes any fermion that can appear in the loop in the triangle diagram and Here, ǫ ξ , ǫ * µ , and ǫ * ν are the polarizations of the Z boson, real photon, and virtual photon, respectively, and p γ is the momentum of the real photon (p 2 γ = 0). 1 Then, following Refs. [1,10], we obtain the indirect amplitude for process Z → V + γ: where Here, g V γ is given by The relative sign between the direct amplitude in Eq. (1) and the indirect amplitude in Eq. (30) disagrees with the relative sign that was found in Ref. [2]. That is, we find that the direct and indirect amplitudes interfere constructively for the process Z → J/ψ + γ and interfere destructively for the processes Z → Υ(nS) + γ.

V. COMPUTATION OF THE DECAY RATES
A. Decay rate The rate for the decay of a Z boson into a vector quarkonium plus a photon is easily seen to be where A dir is given in Eq. (26), A ind is given in Eq. (30b), and we have dropped terms of order m 2 V /m 2 Z . In evaluating the expression for A dir in Eq. (26), we take the hard-scattering scale µ to be m Z , and we take the initial scale µ 0 to be the heavy-quark MS mass m Q . The typical momentum scale of loop corrections to the LCDA and to f V is the pole mass, and, so, the pole mass would be a natural choice for µ 0 . However, the pole mass is ill defined, as we have already mentioned, owing to renormalon ambiguities, and the presence of pole-mass renormalons could impact the convergence of the perturbation series unfavorably in higher orders. Therefore, we choose to take µ 0 = m Q . In applying the Abel-Padé method to the expression for A dir in Eq. (26), we take 100 terms in the eigenfunction expansion and use a 50 × 50 Padé approximant. As we have mentioned, in order to minimize uncertainties in f V , we follow Ref. [3] and compute f V from the leptonic width of the quarkonium, using Eq. (7), instead of using the perturbative expression in Eq. (12).  Table I. We do   not use the values for |Ψ V (0)| 2 in our calculations, but we include them here for purposes of later comparison with the calculations in Ref. [2]. We use the values for |Ψ V (0)| 2 and v 2 V from Refs. [19,20], except in the cases of v 2 Υ(1S) and v 2 Υ(2S) . As was explained in Ref. [11], the uncertainties for v 2 Υ(1S) and v 2 Υ(2S) were probably underestimated in Ref. [20]. We use the larger uncertainties for these quantities that are given in Ref. [11].

C. Sources of uncertainties
In calculating the decay rates, we take into account uncertainties in both the direct and indirect amplitudes, as is described below. We also include the uncertainty in the Z-boson total width in computing branching fractions. We compute the overall uncertainties in the rates by making use of the method that is described in Sec. VIE of Ref. [11]. That is, we find the extrema of the rate for values of the input parameters that lie within a hyperellipse that is centered at the central values of the input parameters and whose semimajor axes have lengths that are equal to the uncertainties in the input parameters.

Direct amplitude
In the direct amplitude, we include the uncertainties that arise from the uncertainties in f V and v 2 V . We also include the uncertainties that arise from uncalculated corrections of order α 2 s , order α s v 2 , and order v 4 . We estimate the uncertainties from these uncalculated corrections, relative to the lowest nontrivial order in the direct amplitude, to be for the real part of the direct amplitude and {[C A α s (m Q )/π] 2 +[v 2 ] 2 } 1/2 for the imaginary part of the direct amplitude. (Note that the real part of the direct amplitude starts in absolute order α 0 s and the imaginary part of the direct amplitude starts in absolute order α s .) The coefficient 1/5 in the v 4 uncertainty in the direct amplitude is the known short-distance coefficient for the order-v 4 correction, which arises from the expression [15] for the 2kth x moment of the LCDA x 2k in terms of the order-v 2k LDME ratio v 2k [see Eq. (9)]: We take v 2 = 0.3 for the J/ψ and v 2 = 0.1 for the Υ(nS) states. We also include an uncertainty of m 2 V /m 2 Z in order to account for uncalculated corrections of order m 2 V /m 2 Z .

Indirect amplitude
In indirect amplitude, we include uncertainties that arise from the uncertainties in the leptonic-decay widths of the quarkonia. We assume that the uncertainties in the leptonicdecay widths are 2.5 % for the J/ψ, 1.3 % for the Υ(1S), and 1.8 % for the Υ(2S) and Υ(3S) states. Again, we include an uncertainty of m 2 V /m 2 Z in order to account for uncalculated corrections of order m 2 V /m 2 Z .

A. Results
Our results for the branching fractions of the Z boson into J/ψ + γ and Υ(nS) + γ are given in Table II. For purposes of comparison, we also show the branching fractions from  Refs. [2] and [3].
Our results for the branching fractions differ considerably from the results in Refs. [2] and [3], in both the central values and in the uncertainties. We now discuss in detail the reasons for those differences.
These differences arise from several sources: (1) we have corrected the value of the scale of Ψ V (0) that was used in Ref. [2]; (2) we have corrected the value of the scale of α s in the order-α s corrections to f V that was used in Ref. [2]; (3) in the direct amplitude, we have absorbed the order-α s and order-v 2 NRQCD corrections to f V in Eq. (12) into an overall factor f V that is determined from the quarkonium electronic decay width, whereas these corrections were computed from the NRQCD expansion and incorporated additively into the direct amplitude in Ref. [2]; (4) we have found a relative sign between the indirect and direct amplitudes that is opposite to the sign that was given in Ref. [2]; (5) Table III, the effects on the branching fractions of the corrections that correspond to these differences are shown. The fractional change in the branching fraction from each correction depends on the order in which the corrections are incorporated into the calculation. In Table III, the fractional changes are computed by incorporating the corrections in the order (1), (2), (3), (4), (5), (6). For each quarkonium state, the product of fractional changes gives the fractional change between our result and that of Ref. [2].
As can be seen from Table III, the effects of corrections (1), (2), (3), and (5) are quite large. However, they tend to cancel each other, and, consequently, our results for branching fractions do not differ so greatly from those in Ref. [2]. We now discuss the corrections to the calculation in Ref. [2] in detail.
In Ref. [2], the decay constant f V was computed by making use of the perturbative expression in Eq. (12). As we have mentioned, this results in greater uncertainties in the calculations. As implemented in Ref. [2], it also leads to shifts in the central values. The  reason for this is that the value for Ψ V (0) that was used in Ref. [2] was extracted from Ref. [21] at the scale m V , while the initial scale µ 0 in Ref. [2] was taken to be m Q . Therefore, the value of Ψ V (0) from Ref. [21] should have been corrected as follows in order to account for the change in the initial scale: The fraction on the right side of Eq. (34) gives correction (1), which produces a correction of +28 % in the rate of Z → J/ψ + γ and a correction of +8 % in the rate of Z → Υ(1S) + γ.
In the expression for the direct amplitude in Ref. [2], there are contributions that are proportional to −8α s (m Z )C F /(4π) − v 2 V /6. These contributions arise when one expresses f V in terms of Ψ V (0), as in Eq. (12). However, the argument of α s should be m Q , rather than m Z . 2 This change of scale accounts for correction (2), which produces a correction of −35 % in the rate of Z → J/ψ + γ and a correction of −17 % in the rate of Z → Υ(1S) + γ.
In the direct amplitude, one can absorb the order-α s and order-v 2 contributions in the NRQCD expansion of f V in Eq. (12) into an overall factor. In our calculation, we express the direct amplitude in terms of the value of f V that one obtains directly from the electronic width of the quarkonium [see Eq. (1)]. As we have mentioned, this approach reduces the size of the uncertainty in the direct amplitude. The effect of absorbing the order-α s and the orderv 2 contributions in the NRQCD expansion of f V into an overall factor f V that is computed from the quarkonium electronic decay rate corresponds to correction (3). Correction (3) changes the rate for Z → J/ψ + γ by −13 % and changes the rate for Z → Υ(1S) + γ by −2 %.
As we have mentioned, our result for the relative sign between the direct and indirect amplitudes disagrees with that in Ref. [2]. Correction (4) accounts for the effects of this change in the relative sign of the indirect amplitude. The numerical effect of correction (4) is very small, changing the rates by only about 2 %, and is insignificant in comparison with the uncertainties in the rates.
In Ref. [2], the resummation of logarithms of m 2 Z /m 2 Q to all orders in α s was estimated to produce a 1.5 % effect in the rate for Z → J/ψ + γ. However, we find a much larger effect, namely, +18 %. We find that the effect of the resummation in the rate for Z → Υ(1S) + γ is +11 %. Correction (5) accounts for these resummation corrections.
In Ref. [2] the initial scale µ 0 = m Q was chosen. As we have explained, we have taken µ 0 = m Q in order to avoid renormalon ambiguities. We have also replaced m Q with m Q in the expression for φ (1) V in Eq. (10c). These differences affect the rate for Z → J/ψ + γ by only +2 % and affect the rate for Z → Υ(1S) + γ by only +1 %. Correction (6) accounts for these differences.
It was claimed in Ref. [2] that only the contributions of the charm-quark, bottom-quark, and τ -lepton loops are important in the indirect amplitude. However, we find that these contributions yield −43 % of the real part of the indirect amplitude in the case of Z → J/ψ+γ and 8 % of the real part of the indirect amplitude in the case of Z → Υ(1S) + γ.
Our uncertainties are considerably smaller than those in Ref. [2]. The differences in uncertainties arise from two principal sources: (1) we have calculated f V from the leptonic width of the quarkonium, using Eq. (7), instead of using the perturbative expression in Eq. (12); and (2) we have taken into account the known short-distance coefficient 1/5 for the order-v 4 corrections in estimating the size of these uncalculated corrections.
The differences between our results for the central values of the branching fractions and those of Ref. [3] arise primarily because our calculations differ from the calculations in Ref. [3] in the following respects: (1) we have included the nonlogarithmic part of the orderα s correction to the LCDA; (2) we have taken µ 0 = m Q for the initial scale, instead of µ 0 = 1 GeV, and we have replaced m Q with m Q in the expression for φ The effects of these differences on the branching fractions are tabulated in Table IV.
As was the case for the corrections to the calculations in Ref. [2], the fractional change in the branching fraction from each correction depends on the order in which the corrections are incorporated into the calculation. In Table IV, the fractional changes are computed by incorporating the corrections in the order (1), (2), (3), (4), (5), (6). For each quarkonium state, the product of fractional changes gives the fractional change between our result and that of Ref. [3], aside from some differences of less than 0.4 % that arise from small differences in the values that are used for the Fermi constant, the heavy-quark pole masses, and the decay constants. As can be seen from Table IV,  The uncertainties in the rates that are given in Ref. [3] are much smaller than the uncertainties that we find. In Ref. [3], uncertainties from uncalculated order-α s corrections are estimated by varying the hard-scattering scale µ. This approach does not take into account uncertainties from uncalculated QCD corrections to the LCDA at the initial scale µ 0 of orders α s (µ 0 ), α 2 s (µ 0 ), and α s (µ 0 )v 2 . We estimate the relative uncertainties from the last two of these sources using the formula to an uncertainty of 8 % in the case of Z → J/ψ + γ and an uncertainty of 2.3 % in the case of Z → Υ(1S) + γ. Our calculation shows that the nonlogarithmic correction to the LCDA of order α s , which is not included in Ref. [3], shifts the rate for Z → J/ψ + γ by about 12 % and shifts the rate for Z → Υ(1S) + γ by about −4 %. In Ref. [3], an uncertainty of about 6 % is given for the rate for Z → J/ψ + γ and an uncertainty of about 3 % is given for the rate for Z → Υ(1S) + γ. Given the uncertainties from uncalculated corrections of order α 2 s (µ 0 ) and α s (µ 0 )v 2 and the shifts from the known corrections of order α s (µ 0 ), we believe that the uncertainties that are given in Ref. [3] are underestimates, especially in the case of the rate for Z → J/ψ + γ.
In Ref. [3], the order-v 2 correction was computed through the use of a model LCDA whose second x moment is adjusted to match the second x moment of the actual order-v 2 correction. The use of a model LCDA circumvents the difficulties of divergent eigenvector series that appear in the resummation of logarithms m 2 Z /m 2 Q . However, the choice of the functional form in the model introduces new uncertainties into the calculation that are not present in a first-principles calculation, such as the calculation in the present paper. In Ref. [22], a model LCDA with the same functional form as the model LCDA in Ref. [3] was used to compute both the order-α s and the order-v 2 correction to the LCDA for the process of Higgs-boson decay to a vector quarkonium plus a photon. It was noted in Ref. [11], that, in this case, the model LCDA does not reproduce the results of the first-principles calculations of the order-α s and the order-v 2 corrections accurately. However, we find that, in the case of the process Z → V + γ, the model LCDA does reproduce the results of firstprinciples calculation of the order-v 2 correction to the LCDA reasonably well. The model LCDA result for the order-v 2 correction differs from the first-principles result by −1.1 % in the case of Z → J/ψ + γ and by +0.8 % in the case of Z → Υ(nS) + γ. This suggests that the difficulties with the model LCDA that were noted in Ref. [11] may arise because of the incorporation of order-α s correction to the LCDA into the model LCDA. We note that the model LCDA contains contributions of order v 4 and higher. As was pointed out in Ref. [11], these contributions are incompatible with the relation between the x moments of the LCDA and the NRQCD LDMEs that is given in Eq. (33). Apparently, the (incorrect) higher-order contributions that are contained in the model LCDA are not numerically significant at the present level of accuracy.

VII. SUMMARY AND DISCUSSION
We have presented a calculation of decay rates for the processes Z → V + γ, where V is one of the vector quarkonia J/ψ or Υ(nS), with n = 1, 2, or 3. Our results for the branching fractions for Z → V + γ are given in Table II. Our calculations contain corrections through relative orders α s and v 2 , as well as logarithms of m 2 Z /m 2 Q , resummed at NLL accuracy to all orders in α s . The use of the Abel-Padé method [11] allows us to compute for the first time the resummation effects for the order-α s corrections to the quarkonium LCDA and to compute from first principles the resummation effects for the order-v 2 corrections to the quarkonium LCDA. The rates for Z → J/ψ + γ and Z → Υ(1S) + γ have been computed previously at lower levels of accuracy [2,3]. Our computations of the rates for the decays Z → Υ(2S) + γ and Z → Υ(3S) + γ are new. We have also verified the expressions for the order-v 2 corrections to the decay rate that are given in Ref. [2].
Our central values for the branching fractions differ from those in Ref. [2] by −10% for the decay Z → J/ψ + γ and by −3% for the decay Z → Υ(1S) + γ. These differences arise principally for the following reasons: (1) we have corrected the value for scale of the quarkonium wave function at the origin that was used in Ref. [2]; (2) we have corrected the value for the scale of α s in the order-α s corrections to the quarkonium decay constant that was used in Ref. [2]; (3) in the direct amplitude, we have replaced the nonrelativistic expansion of f V [in terms of Ψ V (0), α s , and v 2 ] that was used in Ref. [2] with an overall factor f V that is determined from the quarkonium electronic decay rate; (4) we have included resummations of logarithms of m 2 Z /m 2 Q in the direct amplitude, whereas such resummations were not included in the direct amplitude in Ref. [2]. The individual corrections (1)-(4) are quite large, but they tend to cancel each other in the rate. We have also found that the sign of the indirect amplitude, relative to the direct amplitude, is opposite to the sign that is reported in Ref. [2]. The numerical consequences of this change in sign are small.
Our central values for the decay rates differ from those in Ref. [3] by +12 % for the decay Z → J/ψ + γ and by −11 % for the decay Z → Υ(1S) + γ. In the case of the decay Z → J/ψ + γ, most of the shift in the central value occurs because our calculation includes nonlogarithmic corrections to the LCDA of order α s , while the calculation in Ref. [3] does not. In the case of the decay Z → Υ(1S) + γ, the largest difference between our decay rate and that of Ref. [3] occurs because we take the value of v 2 Υ(1S) from the potential-model calculation in Ref. [20], while the calculations in Ref. [3] make use of an estimate v 2 Υ(1S) = 0.1. Other small differences between the results of our calculations and those of Ref. [3] arise for the following reasons: (1) we take the initial scale of the LCDA to be the heavy-quark MS mass, rather than 1 GeV; (2) we include the order-α 2 s contribution to the rate that comes from the absolute square of the order-α s correction to the hard-scattering kernel; (3) we resum logarithms of m 2 Z /m 2 Q at NLL accuracy, rather than LL accuracy; and (4) we include the indirect decay amplitude. We argue that the choice of the heavy-quark mass as the initial scale of the LCDA is more appropriate than the choice 1 GeV because the heavy-quark mass is the typical scale of perturbative loop corrections to the LCDA.
It is argued in Ref. [3] that the value of v 2 Υ(1S) in Ref. [20] cannot be correct because it is negative. However, the minimal-subtraction expression for v 2 Υ(1S) is obtained by subtracting a power divergence. Hence, there is no reason that v 2 Υ(1S) must be nonnegative. One can see that this is so by computing, for example, the minimal-subtraction expression for v 2 for positronium. In the case of positronium, a full calculation, including binding effects, can be carried out reliably in perturbation theory. That computation results in a negative value for v 2 .
The uncertainties in our decay rates are considerably larger than those in Ref. [3]. In Ref. [3], uncertainties that arise from uncalculated corrections of higher orders in α s were estimated by varying the hard-scattering scale µ ∼ m Z . This procedure does not take into account QCD corrections to the LCDA, which reside at a scale µ 0 ∼ m Q and which were not included in the expression for the amplitude in Ref. [3]. Therefore, we believe that the procedure in Ref. [3] underestimates that uncertainties in the rates.
In Ref. [3], the order-v 2 correction to the LCDA were computed by making use of a model for the LCDA whose second x moment, in the narrow-width approximation, agrees with the second x moment of the order-v 2 correction to the LCDA. Such a procedure obviates the use of the Abel-Padé method. However, it introduces model uncertainties that may not be quantifiable. In Ref. [11], it was found that the use of such a model LCDA for both the order-α s and the order-v 2 corrections to the LCDA does not produce accurate results.
However, we have found that, when the model LCDA is used to account only for the order-v 2 correction to the LCDA, it leads to results that differ from our first-principles calculation only by amounts that are, numerically, of order v 4 .
The calculations of the decay rate for Z → V + γ in the present paper improve upon the accuracy of previous theoretical predictions for those rates and give, we believe, more realistic estimates of the theoretical uncertainties. Measurements of the decays Z → V + γ are interesting in their own right as tests of the standard model and as tests of our understanding of the formation of quarkonium bound states in hard-scattering processes.
However, such measurements are also important because they can lead to a better understanding of the experimental difficulties in the observation of quarkonium-plus-photon final states. That understanding may facilitate the observation of the rare decays of the Higgs boson to quarkonium-plus-photon final states, which could yield a first measurement of the Higgs-boson-charm-quark coupling and alternative measurements of the Higgs-bosonbottom-quark coupling.