Neutralino annihilation beyond leading order

High-precision measurements of the relic dark matter density and the calculation of dark matter annihilation branching fractions i or the galactic halo today motivate the computation of the neutralino annihilation cross section beyond leading order. We consider annihilation via squark exchange and parameterize the effective annihilation vertex as a dimension-six operator suppressed by two po squark mass and related to the divergence of the axial vector current of the final-state quarks. Since the axial vector current is conse level in the limit of massless quarks, this dimension-six operator contains a suppression by the quark mass. The quark mass suppre lifted in two ways: (1) by corrections to the dimension-six operator involving the anomalous triangle diagram, and (2) by going to dimensi We address the first of these possibilities by evaluating the anomalous triangle diagram, which contributes to neutralino annihilation to g We relate the triangle diagram via the anomaly equation to the decay of a pseudoscalar into two gluons and use the Adler–Bardeen extract the next-to-leading order (NLO) QCD corrections to χχ → gg from the known corrections to pseudoscalar decay. The strong depen of the dominantχχ → qq̄ cross section on the relative velocity of the neutralinos makes these NLO corrections unimportant at χ decoupling but significant today.  2005 Elsevier B.V.Open access under CC BY license. e de ex ed ymo a ar chy nt tio a rticle ll be onured icroun-


Introduction
The presence of non-baryonic dark matter in the universe is compelling evidence for physics beyond the Standard Model. While many very different models have been proposed to explain the dark matter (wimps, axions, wimpzillas, modified gravity, etc.) [1], the thermal production of stable weaklyinteracting particles with weak-scale mass remains as the most attractive and predictive explanation for the observed dark matter relic abundance, and further allows the solution to the dark matter problem to be linked to the solution to the hierarchy problem and tested at current and future collider experiments. Supersymmetry provides an especially attractive explanation with the lightest supersymmetric particle (LSP) as the dark mat- Another important physics application of these higher-order QCD corrections is the calculation of signals from WIMP annihilations in the galactic halo or in the interior of the Sun. The NLO corrections are potentially important in the evaluation of the branching fractions for the observable gamma ray and neutrino signals, respectively. As we shall see, these NLO corrections turn out to be much more important for the calculation of the observable gamma ray and neutrino signals than in the calculation of the relic density because of the strong dependence of the tree-level annihilation cross-section on the relative velocity of the neutralinos.

Neutralino annihilation cross section
The behavior of the annihilation cross section depends on the composition of the neutralino. Throughout this Letter we assume that the LSP is largely gaugino as motivated by mSUGRA models [4]. The processes that contribute to the cross section up to order α 2 s and one loop are shown in Fig. 1. The tree-level diagram is shown in Fig. 1(a). Figs. 1(b)-(d) show the diagrams with t -channel squark exchange, whereas (e)-(j) show the diagrams with s-channel Z, H 0 , h 0 , A 0 exchanges. The gauge and Higgs bosons couple to the Higgsino part of the LSP and thus their contributions are suppressed for a mostly-gaugino neutralino. 1 The corresponding suppression factors for the s-and p-wave terms in the cross section are given in Table 1.

The anomaly equation
The leading contribution to neutralino annihilation via exchange of a squark of massM, shown in Fig. 1(a), can be reduced to an effective vertex described by a dimension-six operator suppressed byM 2 , 1 The Higgsino fraction suppression can be removed at the cost of going to one loop [5].
where c is a dimensionless coefficient. This dimension-six operator corresponds to taking the leading term in the expansion of the squark propagator in powers of 1/M 2 ; in particular, we work in the limit m 2 χ M 2 .
In the static limit, where the relative velocity of the two neutralinos can be neglected, the operator O 6 is related to the divergence of the axial vector current of the quarksqq: In the massless quark limit, m q = 0, the axial vector current is conserved at tree level, ∂ µ (qγ µ γ 5 q) = 0, and all tree amplitudes due to the dimension-six operator vanish; in particular, radiating additional gluons cannot lift the suppression. Even at the loop level, for example, diagrams involving the exchange of a virtual gluon, the suppression is still valid unless the anomalous triangle diagram is involved. This is the well-known partiallyconserved axial current (PCAC) condition.
Indeed, only through the anomalous loop diagrams is the conservation of the axial vector current violated, in the form [8] (3) µνG (a)µν with 1 2G µν = µναβ G αβ denoting the dual color field strength tensor. The simplest such anomalous diagram is shown in Fig. 1(c). Neglecting the mass of the internal quark q and using the anomaly equation, this diagram can be written in the form for m q m χ . In the opposite limit, m q m χ , the very heavy quark decouples; the top quark contribution can be neglected if m χ 100 GeV. The leading-order (LO) calculation using the anomaly equation was first studied correctly by Ref. [7] in the γ γ channel in QED.
The gluonic decay amplitude of a fundamental pseudoscalar, A 0 → gg, is also related to the anomaly equation. This decay proceeds through a quark loop. In the heavy quark limit, m Q m A , the divergence term on the left-hand side of Eq. (3) Table 1 Dependence of the cross section from each diagram on various suppression factors. The diagrams are shown in Fig. 1. The p-wave contribution to the cross section is shown only when it is important due to suppression or absence of the s-wave component. The column labeled "χ mixing" sketches the dependence of the cross section on the neutralino composition. Interference terms between the various diagrams carry a combination of the suppression factors corresponding to each diagram, and are not shown. We note that the Z-pole structure of diagram (i) is canceled by a numerator factor supplied in accordance with Yang's theorem [6,7].
becomes insignificant, leading to Note that because the Yukawa coupling of A 0 to the quark Q is proportional to the quark mass, the m Q dependence of the A 0 → gg partial width drops out in the limit m Q m A . In contrast, the neutralino pair annihilates into two gluons via the anomaly diagram, which is dominated by light quarks, m q m χ . In this limit, the term proportional to m q on the right-hand side of Eq. (3) becomes insignificant, leading to Thus we see that two seemingly different processes, A 0 → gg through a heavy quark loop and χχ → gg through a light quark loop, are related by the anomaly equation to the same gluonic operator. The Adler-Bardeen theorem [9] guarantees that the anomaly equation, Eq. (3), is valid to all orders of α s . One can take advantage of this anomaly property to obtain the higherorder QCD corrections to χχ → gg from the known results for A 0 → gg at next-to-leading order (NLO) from Ref. [10]. Note that, because of the non-Abelian nature of the gauge field, the above gluonic operator also incorporates tri-gluon amplitudes beyond leading order.

Beyond dimension-six
If we include higher terms in the 1/M 2 expansion, the PCAC constraint will be lifted. The dimension-eight operator corresponding to an amplitude proportional to 1/M 4 survives even in the massless quark limit, m q → 0. One can use this dimensioneight amplitude to calculate the rate of χχ → qqg from diagrams such as Fig. 1(b); a full calculation was done in Ref. [11]. However, the contribution to χχ → qqg from the dimensionsix operator due to the anomaly with a virtual gluon turning into a quark pair ( Fig. 1(c)) suffers lessM suppression and will dominate the dimension-eight term forM m χ even though the order in α s is higher.
The effective vertex of the dimension-eight operator for χχ → qqg is where G (c)αβ is the gluonic field strength for color index c and t c is the corresponding SU(3) generator. This operator has exactly the same form as that in Eq. (6) of Ref. [12] for the χχ → ff γ amplitude computed through explicit expansion of the propagators to order 1/M 4 in the limit m q = 0.
It is interesting to note that the amplitude for χχ → qqg from the dimension-six operator with one gluon splitting into qq (shown in Fig. 1(d)) yields an operator of the same form as in Eq. (7). This allows the interference term between this diagram and the dimension-eight process to be easily obtained.

The anomaly at leading order
Since neutralinos are Majorana in nature, the initial state behaves as a pseudoscalar in the zero-velocity limit [13]. In particular, for v rel = 0 the antisymmetrized neutralino spinors reduce to the projection operators [14] where p 1 and p 2 are the four-momenta of the incoming neutralinos and 2P = p 1 + p 2 .
We work in the limit of zero fermion mass, in which case the off-diagonal terms in the squark mass matrices vanish and the squark mass eigenstates coincide with the electroweak eigenstatesq L andq R . Applying the reduction formula Eq. (8) to the amplitude of the χχ → gg diagram shown in Fig. 1(d) allows us to write the amplitude for the diagram involving squarkq i as where the neutralino-quark-squark couplings are defined for right-and left-handed squarks, respectively, as Here T 3 is the squark isospin, Q is the squark electric charge, and N 11 and N 12 are the bino and wino components of the neutralino as defined in Ref. [15]. We also define P R/L = (1 ± γ 5 )/2 in the usual way as the right-and left-handed projection operators. The external gluon momenta are called k 1 and k 2 , and q is the momentum flowing in the loop. The form factor F (q) µν,ab i for squarkq i and the corresponding internal quark q i is given explicitly by The two terms in F (q) µν,ab i correspond to the two possible directions of fermion flow in Fig. 1(d). The diagrams with the neutralinos crossed are already included through the use of the antisymmetrized spinors in Eq. (8).
After summing over left and right squark states, the amplitude becomes with couplings g r,l given in Eq. (10). After integrating Eq. (12), the piece involving the vector coupling vanishes due to the conservation of vectorial current (CVC): the contracted vector behaves as a divergence, so that the resulting vector coupling to multi-gluon states vanishes. Light quark masses can be neglected in the form factor if m 2 b m 2 χ , and the top quark loop amplitude is suppressed if m 2 χ m 2 t . In these limits, the loop amplitude sums over five massless quarks. Note that in the massless quark limit, the form factor F(q) µν,ab i becomes independent of the quark flavor i. The sum over quarks q i in the loop can then be factorized into a sum over the coupling factors A i times the universal massless form factor F (q) µν,ab i . The gluon production amplitude contains a pseudovector triangle diagram. This can be transformed to a pseudoscalar triangle diagram via the axial anomaly, Eq.
The anomaly is computed by relating it to the decay of a fundamental pseudoscalar Higgs boson A 0 to two gluons via a heavy quark loop. If the mass of the heavy quark Q is sufficiently large, m Q m A 0 , then from Eqs. (5) and (6) we have The amplitude becomes where s = 4m 2 χ and C 0 (0, 0, s, m 2 Q , m 2 Q , m 2 Q ) is a three-point Passarino-Veltman integral [16]. In the limit of heavy quark mass, m 2 Q s, the three-point integral reduces to −1/2m 2 Q . In this limit the dependence on the heavy quark mass drops out and the amplitude becomes where we have used Tr[t a t b ] = (1/2)δ ab . Squaring the amplitude and integrating over phase space gives the leading-order (χχ → gg) annihilation cross section, where in our approximation the sum runs over the five light quarks q i . Our result agrees with that of, e.g., Ref. [11] in the limit m q = 0.

Beyond leading order
As discussed in Section 1.2, we can use the Adler-Bardeen theorem [9] and Eqs. (5) and (6) to obtain the QCD corrections to the χχ → gg annihilation cross section by exploiting the known results for pseudoscalar Higgs decays to gluon pairs, A 0 → gg, beyond leading order.
The NLO QCD corrections to the A 0 → gg partial width were calculated by Spira et al. [10]. In the heavy top quark limit, for which our anomaly relation is valid, the NLO QCD corrections are given by a multiplicative factor [10], times the LO decay rate where µ is the renormalization scale. The integer N f counts the number of quark flavors in the gluon splitting, with N f = 5 for m b m χ m t . The diagrams that contribute to A 0 → gg at LO and NLO are shown in Fig. 2. The NLO final states include gg, ggg, and gqq.   3. χχ → gg diagrams that supply a logarithmic factor to cancel that from Fig. 1(d). The one-loop corrections to the gluon legs include quark and gluon loops.
In χχ → gg at NLO, a divergence occurs for the diagram in Fig. 1(d) when the final-state quarks are soft or collinear, in which case the gluon propagator diverges. This is the source of the logarithmic enhancement factor, log(m 2 χ /m 2 q ), found for this diagram in Ref. [12]. However, this logarithmic term is precisely canceled by the renormalization of the strong coupling due to the quark bubble that appears in the virtual part of the NLO correction, shown as the interference of the diagrams in Fig. 3. This is the familiar cancellation of logarithmic divergences guaranteed by the Kinoshita-Lee-Nauenberg theorem [17]. A similar cancellation occurs for the analogous diagrams in which the soft/collinear quarks are replaced with gluons.
We now invoke the Adler-Bardeen theorem [9] and take over the NLO corrections to A 0 → gg to the χχ → gg process in the zero-velocity limit. In this correspondence, the pseudoscalar mass m A is replaced by the χχ center-of-mass energy, equal to 2m χ in the zero-velocity limit. The NLO correction to the cross section for χχ → gg follows immediately from Eqs. (17) and (18), We give explicitly the result for µ = 2m χ and N f = 5, where we set m χ = 100 GeV. The strong coupling is evaluated based on five-flavor running at the scale µ = 2m χ where it appears both explicitly and within the coefficient A i . We note that the above choice of m χ = 100 GeV is well above the current experimental limit [18].

Numerical results
In this section we examine the validity of our assumption of massless quarks in the loop for the LO χχ → gg calculation, show the improvement in the renormalization scale dependence obtained in going to NLO, and compare the χχ → gg annihilation cross section to that of the leading-order tree-level process χχ → qq through t -channel squark exchange. We have assumed in our use of the anomaly equation that the five light quarks running in the loop for χχ → gg were massless, and we neglected the heavy top quark contribution. In Fig. 4, we test this assumption for the LO cross section by comparing our approximation to the full cross section including the quark mass dependence (we continue to use the 1/M 2 approximation for the squark propagator). We plot the full cross section normalized to our five-massless-quark approximation as a function of m χ . The full formula differs from our approximation by less than 10% for 6 GeV < m χ < 110 GeV. For heavier neutralinos, the top quark loop starts to have a significant effect, with destructive interference occurring between the top loop and the lighter quark loops. Two of the top quarks in the loop go on shell at m χ = m t , leading to the large dip in the cross section. For neutralinos lighter than about 6 GeV, the nonzero mass of the bottom quark begins to play a significant role, leading to the dip at lower masses.
In an all-orders calculation, physical observables cannot depend on the renormalization scale µ. The µ dependence of our predictions is an artifact of computing to a finite order in perturbation theory. The µ dependence can then be used to estimate the size of the uncomputed higher-order corrections. The dependence of the χχ → gg annihilation cross section on the renormalization scale at LO and NLO is shown in Fig. 5(a) for M = 200 GeV, m χ = 100 GeV, and a pure bino neutralino, N 11 = 1, N 1j = 0 (j = 1). 2 Varying µ by a factor of two in either direction from the central value µ = 2m χ yields a scale dependence of ±16% at LO and ±9% at NLO. The corresponding annihilation cross 2 We note here that the renormalization scale dependence of the tree-level χχ → qq cross section arises only from the running quark mass in the s-wave contribution [11]. Since the cross section for this process is dominated by the quark-mass-independent p-wave part during freeze-out in the early universe, the scale dependence is negligible at tree level. sections including scale uncertainty are shown in Fig. 5(b) as a function of m χ .
Exact cross section formulae for all tree-level two-to-two neutralino annihilation processes are given in Ref. [19]. Using the expansion of the thermally averaged cross section in terms of x = T /m χ , one can compare the leading tree-level process, χχ → qq via t -channel squark exchange, to our results for χχ → gg at NLO given in Eq. (20) both for annihilation during freezeout in the early universe, x ∼ 1/20, and for annihilation in the galactic halo today, x ∼ 0 (corresponding to v/c ∼ 10 −3 ). The cross sections for the two annihilation processes are shown in Fig. 6(a) as a function of m χ forM = 200 GeV and a pure bino neutralino, with x = 1/20 and x = 0. For χχ → qq we sum the cross section over the five light final-state quark flavors and neglect χχ annihilation into lepton pairs. We use the running quark mass m q (µ), which serves to resum the leading logarithmic QCD corrections to this process from soft gluon radiation [11], and take the renormalization scale µ = 2m χ . The annihilation cross section for χχ → gg is dominated by the s-channel component, and so we show only one curve for x = 1/20 and x = 0. The tree-level cross section for χχ → qq depends strongly on the relative velocity of the neutralinos, because the s-channel cross section is suppressed by the final-state quark mass. Thus the annihilation cross section at x = 0, which comes only from the s-wave component, is quite small and is comparable to that from χχ → gg. At x = 1/20, on the other hand, the χχ → qq cross section is dominated by the p-wave component and is larger than χχ → gg by almost a factor of 100.
In Fig. 6(b) we show the corresponding K-factors, defined as the ratio of the total cross section, χχ → qq + gg, to the treelevel χχ → qq cross section. This gives a measure of the relative importance of the χχ → gg component of the total annihilation cross section. We see that during freeze-out, x ∼ 1/20, the χχ → gg contribution is quite small and K = 1.01-1.02 over the range of m χ considered. In the present epoch, however, with x ∼ 0, the K-factor is considerably larger, K = 1.3-30 depending on the mass of the neutralino. Such a large K-factor will impact annihilation branching fractions today, changing the gamma ray flux from the galactic halo and the neutrino flux from inside the Sun [1,11]. We note also that because both the χχ → qq and χχ → gg annihilation cross sections have the same leading 1/M 4 dependence on the squark mass, these K-factors will not depend significantly on the common squark mass scale.
We now exhibit the dependence of the annihilation cross section on the neutralino composition. Because we have worked in the zero-quark-mass limit in our calculation of χχ → gg, we have neglected the couplings of Higgsinos to the internal quark loop, which are proportional to the quark mass. We thus consider only mixed bino-wino neutralinos. We can then parameterize the mixing coefficients N 1j in terms of a bino-wino   In Fig. 7 we again compare the leading tree-level process, χχ → qq via t-channel squark exchange, to our results for χχ → gg at NLO. The cross sections for the two annihilation processes are shown in Fig. 7(a) as a function of the binowino mixing angle φ for m χ = 50 GeV, andM = 200 GeV, with x = 1/20 and x = 0. Again we take the renormaliza-tion scale µ = 2m χ . There is a large enhancement of both annihilation cross sections if the neutralino has a large wino component, φ ∼ 90 • , due to the stronger coupling of the wino to a quark-squark pair. In Fig. 7(b) we show the corresponding K-factors. The nontrivial dependence of the Kfactors on φ arises from the interference among the five light quark loop diagrams that contribute to χχ → gg. By contrast, there is no interference between the amplitudes for χχ → qq with different flavor squarks in the t -channel be-cause the internal squark flavor is fixed by the external quark flavor. Finally, we note that the next-to-next-to-leading order (NNLO) corrections to A 0 → gg have been computed in Ref. [20]. One may be tempted to take this correction over to χχ → gg in the same way as the NLO correction. However, at NNLO the A 0 → gg decay receives a contribution from the interference between the G (a) µνG (a)µν operator and a ∂ µq γ µ γ 5 q operator generated at two-loop level. Once effective operators other than G (a) µνG (a)µν appear, our use of the anomaly equation no longer applies. A proper treatment of χχ → gg at NNLO would thus require a new calculation of the operator matching conditions and renormalization.

Conclusions
We reviewed the dependence of the main neutralino annihilation processes on various suppression factors-the Higgsino fraction, the quark and squark masses, and the relative neutralino velocity-and identified the dominant χχ → qq annihilation process for a gaugino-like neutralino as due to a dimension-six operator in the zero-velocity limit. This dimension-six operator contains the divergence of the axial vector current of the quarks qq, which leads to the wellknown quark mass suppression of the annihilation cross section. This quark mass suppression can be lifted in two ways: either through corrections to the dimension-six operator involving the anomalous triangle diagram, or by going to dimensioneight. We focused on the anomalous triangle diagram, which describes neutralino annihilation to gluon pairs. In the approximation of massless quarks running in the loop, the anomaly equation relates χχ → gg to the seemingly unrelated process of pseudoscalar decay to gluon pairs via a very heavy quark loop. We used this relation to compute χχ → gg in terms of the decay process A 0 → gg. Further, taking advantage of the Adler-Bardeen theorem which guarantees that the anomaly equation is valid to all orders in α s , we extracted the NLO QCD corrections to χχ → gg from the known corresponding results for A 0 → gg and wrote them as a simple multiplicative factor that can be easily inserted into numerical neutralino annihilation codes. For m χ = 100 GeV and µ = 2m χ , the NLO QCD corrections increase the annihilation cross section to gg by 62%. The NLO corrections also reduce the residual renormalization scale dependence of the χχ → gg annihilation cross section from ±16% to ±9%.
Our NLO results were computed in the approximation m b m χ m t . This approximation yields a LO χχ → gg cross section within 10% of the exact result for 6 GeV < m χ < 110 GeV. We finally compared our results for χχ → gg at NLO to the dominant χχ → qq cross section in this neutralino mass range. For neutralino annihilation during freeze-out in the early universe, our results for χχ → gg at NLO constitute only 1-2% of the dominant cross section for a gaugino-like neutralino and are thus of little importance for computing the relic neutralino abundance. However, for neutralino annihilation at the present time the relative neutralino velocity is much lower, leading to a much smaller tree-level χχ → qq cross section. In this situation, the χχ → gg cross section can be as large or larger than χχ → qq, so that NLO corrections can have a significant impact on the computation of gamma ray and neutrino fluxes from neutralino annihilation in the galactic halo and inside the Sun, respectively.