text:
# Higgs boson production in bottom quark fusion at next-to-next-to-leading order

## I Introduction

The search for the Higgs boson will be a top priority of the CERN Large Hadron Collider (LHC). The LHC’s discovery potential for the Standard Model Higgs boson fully covers the mass range from the experimental lower bound established by the LEP experiments ($`M_H114`$ GeV) up to $`M_H1`$ TeV, beyond which the concept of the Higgs boson as an elementary particle becomes questionable. In addition, the Fermilab Tevatron could find evidence for or even discover the Higgs boson if $`M_H180`$ GeV and if sufficient luminosity lumimoni  can be collected.

The theoretical description of the signal processes for Standard Model Higgs boson production is under good control. For a review, see Ref. Carena:2002es . The dominant production mode is gluon fusion, for which the next-to-next-to-leading order (NNLO) corrections are now available Harlander:2002wh ; Anastasiou:2002yz  and have recently been reconfirmed by Ref. Ravindran:2003um . The radiative corrections for the weak boson fusion channel Han:1992hr  and associated production with a weak gauge boson Han:1991ia  have been known for several years, rendering the theoretical uncertainty in these processes very small. Recently, next-to-leading order (NLO) corrections have also been evaluated for Higgs boson production in association with top quarks Reina:2001sf ; Beenakker:2001rj ; Dawson:2002tg ; Beenakker:2002nc , resulting in a drastic reduction of the scale uncertainty.

These results can be used for supersymmetric Higgs boson production as well. However, because of the enriched particle spectrum in supersymmetric extensions of the Standard Model, they provide only a part of the full production rate in general. Additional contributions arise through intermediate supersymmetric partners Dawson:1996xz  and modified couplings of the Standard Model particles. In order to avoid unnecessary generalizations, we will focus on the Minimal Supersymmetric Standard Model (MSSM) for the rest of this paper (see, e.g., Ref. Martin:1997ns  for an outline of the MSSM). The extent to which our results can be transferred to other models should be clear from this discussion.

The MSSM contains two Higgs doublets, one giving mass to up-type quarks and the other to down-type quarks. The associated vacuum expectation values are labeled $`v_u`$ and $`v_d`$, respectively, and they determine the MSSM parameter $`\mathrm{tan}\beta v_u/v_d`$. After spontaneous symmetry breaking, there are five physical Higgs bosons, whose mass eigenstates are denoted by $`h`$ (“light scalar”), $`H`$ (“heavy scalar”), $`H^\pm `$ (“charged scalars”), and $`A`$ (“pseudoscalar”). One interesting consequence of this more complicated Higgs sector is that, compared to the Standard Model, the bottom quark Yukawa coupling can be enhanced with respect to the top quark Yukawa coupling. In the Standard Model, the ratio of the $`t\overline{t}H`$ and $`b\overline{b}H`$ couplings is given at the tree-level by $`\lambda _t^{\mathrm{SM}}/\lambda _b^{\mathrm{SM}}=m_t/m_b35`$. In contrast, in the MSSM, it depends on the value of $`\mathrm{tan}\beta `$. At leading order,

$$\frac{\lambda _t^{\mathrm{MSSM}}}{\lambda _b^{\mathrm{MSSM}}}=f_\varphi (\alpha )\frac{1}{\mathrm{tan}\beta }\frac{m_t}{m_b},$$
(1)
with

$$f_\varphi (\alpha )=\{\begin{array}{cc}\mathrm{cot}\alpha ,\hfill & \varphi =h,\hfill \\ \mathrm{tan}\alpha ,\hfill & \varphi =H,\hfill \\ \mathrm{cot}\beta ,\hfill & \varphi =A,\hfill \end{array}$$

where $`\alpha `$ is the mixing angle between the weak and the mass eigenstates of the neutral scalars. A value of $`\mathrm{tan}\beta `$ as large as $`3040`$ could be accommodated fairly naturally in the MSSM. Such an enhancement would have (at least) two important consequences. The first is that in the gluon fusion mode it is no longer sufficient to consider top quark loops as the only mediators between the Higgs boson and the gluons; one must also include the effects of bottom quark loops (see Fig. 1). Since one cannot justify the use of an effective field theory in which the bottom quarks are integrated out, this involves computing massive multi-loop diagrams. While the massive NLO calculation (including massive two-loop virtual diagrams) was performed some time ago Spira:1995rr , the NNLO result (requiring up to three loops for the virtual correction) is still beyond the limits of current calculational technology.

The second consequence is that Higgs boson production in association with bottom quarks can become an important channel: $`p\stackrel{~}{p}b\overline{b}\varphi `$ ($`\stackrel{~}{p}\{p,\overline{p}\}`$ and $`\varphi \{h,H,A\}`$ here and in what follows). At first sight, the evaluation of the corresponding cross section is in close analogy to the process $`p\stackrel{~}{p}t\overline{t}\varphi `$. But this is only true if the bottom quarks are observed in the detector, and are thus restricted to large transverse momenta. If at least one of the bottom quarks escapes detection, the production rate must be integrated over all transverse momenta of this bottom quark. Since the Higgs boson is much heavier than the bottom quark, this integration leads to collinear logarithms, $`\mathrm{ln}(m_b/M_\varphi )`$, which require a more careful analysis than in the case of $`t\overline{t}\varphi `$ production.

The subject of this paper is the precise prediction of the total cross section for Higgs boson production in association with bottom quarks, where neither bottom quark need be detected. This requires integrating over the transverse momenta of both final state bottom quarks. Each integration gives rise to collinear logarithms of the kind mentioned above. Since the bottom quarks may remain undetected, it is more appropriate to view our result as a part of the total inclusive Higgs production rate $`\sigma (p\stackrel{~}{p}\varphi +X)`$. In order to emphasize this point, we shall henceforth denote the fully inclusive process mediated through bottom–antibottom annihilation as $`p\stackrel{~}{p}(b\overline{b})\varphi +X`$.

Our calculation is based on the approach of Refs. Dicus:1989cx ; Dicus:1998hs ; Maltoni:2003pn , where the leading-order (LO) partonic process is taken to be $`b\overline{b}\varphi `$. We will refer to this as the variable flavor number scheme (VFS) approach in what follows, as opposed to the fixed flavor number scheme (FFS) approach, where the tree-level process $`ggb\overline{b}\varphi `$ is taken as the lowest-order contribution and bottom quarks cannot appear in the initial state. The initial state bottom quarks in the VFS approach arise (predominantly) from gluon splitting in the proton, parametrized in terms of bottom quark parton distributions Barnett:1988jw ; Olness:1988ep ; Aivazis:1994pi ; Thorne:1998ga ; Kramer:2000hn ; Buza:1996ie ; Buza:1998wv ; Chuvakin:1999nx ; Chuvakin:2000qc ; Chuvakin:2000jm ; Chuvakin:2000zj . In this way, the large collinear logarithms that arise due to the fact that the colliding gluons carry a momentum of the order of $`M_\varphi /2m_b`$ can be resummed through DGLAP evolution. The convolution of these bottom quark densities with the partonic cross section leads to the hadronic cross section $`\sigma \left[p\stackrel{~}{p}(b\overline{b})\varphi +X\right]`$.

This process has been a subject of interest for some time. It is currently known up to NLO in the VFS approach Eichten:1984eu ; Gunion:1987pe ; Olness:1988ep ; Dicus:1989cx ; Dicus:1998hs ; Maltoni:2003pn . In the FFS approach, the calculation is analogous to $`t\overline{t}\varphi `$ production, which is also known to NLO Raitio:1979pt ; Ng:1984jm ; Kunszt:1984ri ; Reina:2001sf ; Beenakker:2001rj ; Dawson:2002tg ; Beenakker:2002nc . The case where one bottom quark is tagged has been computed at NLO in Ref. Campbell:2002zm ; in that case, the LO process in the VFS approach is $`bgb\varphi `$.

There has been an ongoing discussion as to the relative merits of the VFS and FFS approaches Dicus:1998hs ; Balazs:1998sb ; Campbell:2002zm ; Maltoni:2003pn ; Rainwater:2002hm ; Spira:2002rd , especially because the results of the two approaches disagree by more than an order of magnitude for the scale choice $`\mu _F=\mu _R=M_\varphi `$, where $`\mu _F`$ and $`\mu _R`$ are the factorization and the renormalization scales, respectively. Recently, it has been argued Maltoni:2003pn , that the proper factorization scale for this process should be $`\mu _FM_\varphi /4`$ instead of $`M_\varphi `$. Indeed, for this choice, the disagreement between the VFS and the FFS approach is significantly reduced. The result of our paper demonstrates the self-consistency of the VFS approach and confirms the proposed factorization scale of Ref. Maltoni:2003pn  as the appropriate choice.

Given the considerations above, the motivations for a NNLO calculation are manyfold. One is to examine the assertion that $`\mu _FM_\varphi /4`$ is the proper choice for this process at NLO Plehn:2002vy ; Maltoni:2003pn . If the higher order corrections are minimal at that scale, this would be a strong argument in favor of the validity of the VFS approach to $`(b\overline{b})\varphi `$ production. A second motivation, as will be discussed in more detail below, is that the NNLO terms play an exceptional role in the VFS approach to $`(b\overline{b})\varphi `$ production due to the fact that they are the first to consistently include the “parent” process, $`ggb\overline{b}\varphi `$, and thereby sample the same range of bottom quark transverse momenta as the LO FFS approach. A third and perhaps dominant motivation for the NNLO calculation is to reduce the sensitivity of the calculation to the unphysical scale parameters $`\mu _F`$ and $`\mu _R`$, thereby removing a significant source of uncertainty from the theoretical prediction.

In this paper we will present results for the process $`p\stackrel{~}{p}(b\overline{b})\varphi +X`$ at NNLO. As will be shown, they nicely meet all expectations concerning their dependence on the renormalization and factorization scales, thus providing a solid prediction for the total cross section at the LHC and the Tevatron. The inclusive production cross section could have phenomenological implications for the observation of the supersymmetric $`H`$ and $`A`$ bosons, for example, in the $`H/A\mu ^+\mu ^{}`$ decay mode.

The organization of the paper is as follows: In Sec. II we discuss the VFS approach to computing the $`p\stackrel{~}{p}(b\overline{b})\varphi +X`$ process and its motivations. In Sec. III we describe the actual calculation and in Sec. IV we present our numerical results. Analytic results for the partonic cross sections are presented in the appendix.

## II Theoretical description of the production rate

In the FFS approach, the LO partonic process for the production of a Higgs boson in association with a bottom quark pair is of order $`\alpha _s^2`$. A few typical diagrams are shown in Fig. 2. Because of the large mass difference between the bottom quark and the Higgs boson, the total cross section contains large logarithms of the form

$$\begin{array}{c}\hfill l_b\mathrm{ln}(m_b^2/\mu _\varphi ^2),\end{array}$$
(2)
where $`\mu _\varphi `$ is of the order of $`M_\varphi `$. More precisely, every on-shell gluon that splits into a $`b\overline{b}`$ pair with an on-shell bottom quark generates one power of that logarithm. Thus, Figs. 2(a) and 2(b) generate two and one power of $`l_b`$, respectively, while Figs. 2(c) and 2(d) do not generate any $`l_b`$ terms. Furthermore, each higher order in perturbation theory brings in another power of $`l_b`$ due to the radiation of gluons from bottom quarks.

Because $`l_b\mathrm{ln}(m_b^2/M_\varphi ^2)`$ is rather large, $`\alpha _sl_b`$ is not a good expansion parameter. It would be better to re-organize the perturbative series such that terms like $`(\alpha _sl_b)^n`$ are resummed to all orders in $`n`$. This resummation can be achieved by introducing bottom quark parton distribution functions which contain all the collinear terms arising from the splitting of gluons into $`b\overline{b}`$ pairs Barnett:1988jw ; Olness:1988ep ; Aivazis:1994pi ; Thorne:1998ga ; Kramer:2000hn ; Buza:1996ie ; Buza:1998wv ; Chuvakin:1999nx ; Chuvakin:2000qc ; Chuvakin:2000jm ; Chuvakin:2000zj . This constitutes the motivation for using the VFS approach.

Convolving the tree-level process of Fig. 3(a) with these bottom quark distributions will resum the leading logarithms of the form $`(\alpha _sl_b)^2(\alpha _sl_b)^n`$, $`n0`$. In order to retain subleading logarithms, one has to compute higher orders. For example, including the NLO contributions with all the relevant subprocesses ($`b\overline{b}h`$, $`b\overline{b}hg`$, $`gbhb`$, and $`g\overline{b}h\overline{b}`$), resums the terms of order $`\alpha _s^2l_b(\alpha _sl_b)^n`$, $`n0`$. In order to retain all powers of $`l_b`$ at order $`\alpha _s^2_n(\alpha _sl_b)^n`$, it is necessary to evaluate the cross section up to NNLO.

Let us briefly review the idea of the VFS in its simplest form. Assume $`n_{\mathrm{}}=n_f1`$ massless quark flavors and one massive quark flavor of mass $`m_h`$. First, one defines parton densities $`f_i^{(n_{\mathrm{}})}(x,Q^2)`$ for the gluon ($`i=g`$) and the massless flavors ($`i=1,\mathrm{},n_{\mathrm{}}`$) in the standard way, obeying DGLAP evolution with $`n_{\mathrm{}}`$ active flavors. The heavy quark density $`f_{n_f}^{(n_{\mathrm{}})}(x,Q^2)`$ is assumed to vanish. Partonic processes involving the heavy quark should be evaluated by keeping the heavy quark mass. This is called the $`n_{\mathrm{}}`$-flavor scheme.

At a certain scale $`\mu _h^2`$, one relates the $`n_{\mathrm{}}`$\- to the $`n_f`$-flavor scheme by defining initial conditions for new parton densities $`f_i^{(n_f)}`$ in terms of the $`f_i^{(n_{\mathrm{}})}`$:

$$\begin{array}{cc}\hfill f_i^{(n_f)}(x,Q^2=\mu _h^2)& =\underset{j}{}C_{ij}(\mu _h/m_h)f_j^{(n_{\mathrm{}})}(x,Q^2=\mu _h^2),\hfill \\ \hfill i& =g,1,\mathrm{},n_fj=g,1,\mathrm{},n_{\mathrm{}}.\hfill \end{array}$$
(3)
The $`C_{ij}`$ are determined by the requirement that physical quantities are the same (up to higher orders in $`\alpha _s`$) in both the $`n_{\mathrm{}}`$\- and the $`n_f`$-flavor scheme. (This requirement may be implemented asymptotically or using mass dependent terms Aivazis:1994pi ; Thorne:1998ga ; Buza:1998wv ; Chuvakin:1999nx .) Above the matching scale, the DGLAP evolution of the $`f_i^{(n_f)}(x,Q^2)`$ ($`i=g,1,\mathrm{},n_f`$) is performed with $`n_f`$ active flavors.

In general, one assumes the $`n_{\mathrm{}}`$-flavor scheme at scales $`Q^2m_h^2`$ and switches to the $`n_f`$-flavor scheme at larger values of $`Q^2`$. It is also convenient to choose the matching scale in Eq. (3) as $`\mu _h^2=m_h^2`$, which avoids the occurrence of logarithms of $`m_h/\mu _h`$.

For our purposes, $`m_h^2/Q^2m_b^2/M_\varphi ^20.003`$, so threshold effects from the matching prescription should be minimal. For the same reason, it is justified to neglect the bottom quark mass in the partonic process (apart from Yukawa couplings, of course). Indeed, masses must be neglected for initial state bottom quarks in order to avoid violation of the Bloch-Nordsieck theorem (and the consequent failure to fully cancel infrared divergences) at NNLO and beyond Catani:1988xy .

In order to make the following discussion more transparent, let us write the fully inclusive $`(b\overline{b})\varphi `$ production rate in the VFS approach schematically as follows:

$$\begin{array}{cc}\hfill \sigma (p\stackrel{~}{p}& (b\overline{b})\varphi +X)=\hfill \\ & \underset{n=0}{\overset{\mathrm{}}{}}\left(\alpha _sl_b\right)^n\{\alpha _s^2[c_{n0}l_b^2+c_{n1}l_b+c_{n2}]\hfill \\ & +\alpha _s^3d_{n3}+\alpha _s^4d_{n4}+\alpha _s^5d_{n5}+\mathrm{}\}.\hfill \end{array}$$
(4)
The sum over $`n`$ is implicit in the parton densities. A LO calculation determines the coefficients $`c_{n0}`$, while NLO adds the coefficients $`c_{n1}`$. Note that the subprocess $`bgb\varphi `$ does not fully determine the coefficients $`c_{n1}`$; in order to obtain the correct resummation at $`\alpha _s^nl_b^{n1}`$ ($`n2`$), one has to include the real and virtual corrections to the $`b\overline{b}\varphi `$ subprocess as well. In the same way, the sum of all subprocesses that contribute at second order determines $`c_{n2}`$ ($`n0`$), and thus all terms associated with the order $`\alpha _s^2(\alpha _sl_b)^n`$. The NNLO result is thus the first to include all terms of order $`\alpha _s^2`$ (as well as higher order terms resummed in the parton distribution functions). Higher orders in perturbation theory correspond to the coefficients $`d_{nk}`$; their $`l_b`$ terms are — formally — completely contained in the parton densities. This illustrates once more the exceptional role of the NNLO corrections in this approach.

The leading order terms were evaluated by Eichten et al. Eichten:1984eu . The leading $`bg`$ and $`gg`$ initiated processes (Figs. 3(c) and 2(a)) were subsequently added by Dicus and Willenbrock Dicus:1989cx . Ten years later, Dicus et al. Dicus:1998hs  (see also Ref. Balazs:1998sb ) computed the full NLO contribution to $`b\overline{b}\varphi `$ (and related subprocesses), which leads to the single logarithmically suppressed term $`c_{n1}`$ for all $`n0`$. Setting the renormalization and the factorization scale equal to $`M_\varphi `$ ($`\mu _R=\mu _F=M_\varphi `$), they found that the $`𝒪(\alpha _s)`$-corrections to the $`b\overline{b}\varphi `$ subprocess and the contribution from the tree level $`bg\varphi b`$ subprocess are quite large, but of opposite sign. This leads to large cancellations which are particularly drastic at the LHC. They also observed that the contribution from $`bg\varphi b`$ becomes especially large at Higgs boson masses below $`150`$ GeV, meaning that logarithmically suppressed terms become important in this region.

Recently, Maltoni et al. Maltoni:2003pn  revisited the NLO calculation in the light of Ref. Plehn:2002vy , which gives an argument for the proper choice of the factorization scale when using the bottom quark density approach. Following that argument, they determined the factorization scale for the $`(b\overline{b})\varphi `$ process to be $`\mu _FM_\varphi /4`$. With this choice, both the NLO corrections from the $`gb`$ and the $`b\overline{b}`$ initiated process turn out to be very well behaved.

As we will show in Sect. IV, the behavior of the NNLO corrections confirms this choice of scale at NLO, in the sense that the perturbation theory up to NNLO is very well behaved for this choice. This supports the method of Refs. Plehn:2002vy ; Maltoni:2003pn  (see also Ref. Boos:2003yi ) for determining the factorization scale in the VFS approach at lower orders. For the process under consideration, it turns out that the dependence on the unphysical scales of the NNLO result is so weak that the discussion on the proper scale choice becomes irrelevant. The overall conclusion is that the prediction for Higgs boson production in bottom quark fusion is now under very good control.

## III Outline of the calculation

As discussed before, we will neglect the bottom quark mass everywhere except in the Yukawa couplings. The calculation is thus completely analogous to, say, Drell-Yan production of virtual photons Hamberg:1991np ; Harlander:2002wh : One evaluates virtual and real corrections to Higgs production in $`b\overline{b}`$, $`gb`$, $`gg`$, $`bb`$, $`qb`$ and $`q\overline{q}`$ scattering (and the charge conjugated processes) and then performs ultraviolet renormalization and mass factorization.

The subprocesses to be evaluated at the partonic level are given as follows ($`q\{u,d,c,s\}`$):
$``$ up to two loops: $`b\overline{b}\varphi `$ \[Fig. 3(a)\]  $``$ up to one loop: $`b\overline{b}\varphi g`$, $`gb\varphi b`$ \[Figs. 3(b), 3(c)\]  $``$ at tree level:   $`b\overline{b}\varphi gg`$, $`b\overline{b}\varphi q\overline{q}`$, $`b\overline{b}\varphi b\overline{b}`$, $`gb\varphi gb`$, $`bb\varphi bb`$, $`bq\varphi bq`$ \[Figs. 4(a)–4(f)\] $`gg\varphi b\overline{b}`$ \[Figs. 2(a)–2(c)\], $`q\overline{q}\varphi b\overline{b}`$ \[Fig. 2(d)\]    We compute the two-loop virtual terms by employing the method of Refs. Baikov:2000jg ; Harlander:2000mg , which maps them onto three-loop two-point functions. In this way, they can be reduced to a single master integral, using the reduction formulas given in Refs. Chetyrkin:1981qh ; Larin:1991fz . The master integral has been computed in Ref. Gonsalves:1983nq . The pole parts of this “$`b\overline{b}\varphi `$ form factor” can be compared to the general formula of Refs. Catani:1998bh ; Sterman:2002qn  which provides a welcome check. For the generation of the diagrams we use QGRAF Nogueira:1993ex  as embedded in the automated system GEFICOM Geficom ; Steinhauser:2002rq  (see also Ref. Harlander:1998dq ).

The one-loop single emission processes are obtained by computing analytically the full one-loop amplitudes, which are then interfered with the amplitudes of the tree-level processes. The two-particle phase space integrals are also computed analytically.

For the tree-level double emission processes, we express the matrix elements and phase space measures in terms of the variable $`x=M_\varphi ^2/\widehat{s}`$, where $`\widehat{s}`$ is the center of mass energy. Then we expand the integrands in terms of $`(1x)`$ Harlander:2002wh . This leaves us with only one nontrivial phase space integral, independent of the order of the expansion. The regular integrands and the finite integration region ensure the validity of the interchange of integration and expansion. Keeping of the order of ten terms in the expansion in $`(1x)`$ leads to a hadronic result that is already phenomenologically equivalent to the analytic result. By evaluating the expansion up to sufficiently high orders, however, one can invert the series Kilgore:2002sk  by mapping the expansion onto a set of basis functions. The latter can be deduced from the known NNLO Drell-Yan result Hamberg:1991np ; Harlander:2002wh .

All algebraic manipulations are performed with the help of the program FORM Vermaseren:2000nd .

For a consistent treatment of the NNLO process, it is not sufficient to evaluate only the partonic cross section at NNLO. Another ingredient is the proper parton densities, obeying NNLO Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution. At present, only approximate evolution kernels are known, derived from moments of the structure functions Retey:2000nq ; Larin:1994vu ; Larin:1997wd . On this basis, approximate NNLO parton distribution sets have been evaluated Martin:2002dr . We use this set in all of our numerical analyses below. Once parton distributions that use exact NNLO evolution become available, it is a straightforward task to update the analysis using the partonic results presented in Appendix A.

Let us now turn to the underlying interaction and the renormalization of the partonic results. We ignore the bottom quark mass and the electroweak interactions, so for our purposes, the Lagrangian is:

$$\begin{array}{cc}\hfill _{b\overline{b}\varphi }& =\frac{1}{4}F_{\mu \nu }^aF^{a\mu \nu }+\underset{q}{}\overline{q}i\text{/}Dq+\overline{b}i\text{/}Db\lambda _b^\mathrm{B}\overline{b}\varphi b,\hfill \end{array}$$
(5)
where $`F_{\mu \nu }^a`$ is the gluon field strength tensor, $`D_\mu `$ is the QCD covariant derivative, and the sum runs over the quarks $`u,d,s,c`$. $`\lambda _b^\mathrm{B}`$ is a bare bottom Yukawa coupling constant. In the modified minimal subtraction $`\left(\overline{\text{MS}}\right)`$ scheme, the scalar coupling is renormalized such that<sup>1</sup><sup>1</sup>1We refrain from quoting terms proportional to $`\gamma _\mathrm{E}`$ and $`\mathrm{ln}4\pi `$ that drop out of $`\overline{\text{MS}}`$-renormalized quantities.

$$\begin{array}{cc}\hfill \lambda _b^\mathrm{B}& \lambda _bZ_m(\alpha _s),\left[\varphi =h,H\right]\hfill \\ \hfill Z_m(\alpha _s)& =1\frac{\alpha _s}{\pi }\frac{1}{ϵ}\hfill \\ \hfill +& \left(\frac{\alpha _s}{\pi }\right)^2\left[\frac{1}{ϵ^2}\left(\frac{15}{8}\frac{n_f}{12}\right)+\frac{1}{ϵ}\left(\frac{101}{48}+\frac{5}{72}n_f\right)\right]\hfill \\ & +𝒪(\alpha _s^3),\hfill \end{array}$$
(6)
where $`ϵ=(4D)/2`$ and $`D`$ is the number of space-time dimensions in which we evaluate the loop (and phase-space) integrals. $`Z_m(\alpha _s)`$ is identical to the quark mass renormalization constant of QCD Chetyrkin:1997dh ; Vermaseren:1997fq . Here and in the following, we use the short hand notations $`\lambda _b\lambda _b^{(n_f)}(\mu _R)`$ and $`\alpha _s\alpha _s^{(n_f)}(\mu _R)`$ for the $`\overline{\text{MS}}`$-renormalized Yukawa and strong coupling constants, respectively. $`\mu _R`$ is the renormalization scale, and $`n_f`$ is the number of “active” quark flavors. We will set $`n_f=5`$ in our numerical analyses.

There are (at least) two methods of obtaining the result for pseudoscalar production. The first is to replace the Yukawa interaction term in Eq. (5) with a pseudoscalar interaction,

$$\lambda _b^\mathrm{B}\overline{b}\varphi bi\lambda _b^\mathrm{B}\overline{b}\varphi \gamma _5b,$$
(7)
and proceed by direct calculation.

The second method is to exploit the chiral symmetry of the bottom quarks in Eq. (5), which implies that we are free to perform independent left-handed and right-handed phase rotations of the bottom quarks. If we perform the rotation

$$b_Rib_R^{}b_Lb_L^{},$$
(8)
the Lagrangian becomes

$$\begin{array}{cc}\hfill _{b\overline{b}\varphi }& \frac{1}{4}F_{\mu \nu }^aF^{a\mu \nu }+\underset{q}{}\overline{q}i\text{/}Dq+\overline{b}^{}i\text{/}Db^{}i\lambda _b^\mathrm{B}\overline{b}^{}\varphi \gamma _5b^{},\hfill \end{array}$$
(9)
and we find the same interaction Lagrangian as in Eq. (7). This implies that the cross section for pseudoscalar Higgs boson production, written in terms of the Yukawa coupling $`\lambda _b`$, is identical to the cross section for scalar Higgs boson production to all orders in $`\alpha _s`$.

Following the prescription of Larin Larin:1993tq <sup>2</sup><sup>2</sup>2Note that only the e-print of Ref. Larin:1993tq  discusses the renormalization of the pseudoscalar current. for the treatment of $`\gamma _5`$ in dimensional regularization, we have performed the direct calculation through NNLO and find that this is indeed the case.

Even in the direct calculation, one can see that this identity will hold to all orders in $`\alpha _s`$ with the following argument. If we square the amplitude before computing loop integrals, all fermion lines are closed loops. The fact that we set the bottom quark mass to zero means that both Higgs boson vertices (in both the scalar and pseudoscalar cases) must appear on the same fermion line. If only one Higgs vertex were to appear on a fermion line, there would be an odd number of $`\gamma `$ matrices in the fermion trace which would therefore vanish. In the pseudoscalar case, this means that nonvanishing fermion traces must contain either zero or two $`\gamma _5`$ matrices. The prescription of Larin Larin:1993tq  allows one to assume anticommutativity of $`\gamma _5`$ and identify $`\gamma _5^2=1`$ when two $`\gamma _5`$-matrices are on the same fermion line. Thus, the $`\gamma _5`$-matrices can be eliminated and we see that the calculation for pseudoscalar Higgs boson production is identical, diagram by diagram of the squared amplitude, to that for scalar Higgs boson production, apart from the different Yukawa couplings.

For the sake of completeness, let us remark that the Standard Model value for the coupling constant is given by $`\lambda _b=\sqrt{2}m_b/v`$, where $`v246`$ GeV is the vacuum expectation value for the Higgs boson field, and $`m_b`$ is the running $`\overline{\text{MS}}`$ mass of the bottom quark, $`m_b(\mu _R)`$, evaluated at the renormalization scale $`\mu _R`$. In the MSSM we have

$$\begin{array}{cc}\hfill \lambda _b& =\{\begin{array}{cc}\sqrt{2}\frac{m_b}{v}\frac{\mathrm{sin}\alpha }{\mathrm{cos}\beta },\hfill & \varphi =h,\hfill \\ \sqrt{2}\frac{m_b}{v}\frac{\mathrm{cos}\alpha }{\mathrm{cos}\beta },\hfill & \varphi =H,\hfill \\ \sqrt{2}\frac{m_b}{v}\mathrm{tan}\beta ,\hfill & \varphi =A.\hfill \end{array}\hfill \end{array}$$
(10)
The renormalized partonic results have a dependence on the unphysical scales $`\mu _F`$ and $`\mu _R`$, both explicitly in terms of logarithms, and implicitly through the parameters $`\alpha _s(\mu _R)`$ and $`\lambda _b(\mu _R)`$. The variation of $`\alpha _s`$ and $`\lambda _b`$ with $`\mu _R`$ is governed by the renormalization group equations (RGEs)

$$\mu _R^2\frac{\mathrm{d}}{\mathrm{d}\mu _R^2}a_s=\beta (a_s)a_s,\mu _R^2\frac{\mathrm{d}}{\mathrm{d}\mu _R^2}\lambda _b=\gamma ^m(a_s)\lambda _b,a_s\frac{\alpha _s}{\pi },$$
(11)
where

$$\begin{array}{cc}\hfill \beta (a_s)& =a_s\beta _0a_s^2\beta _1a_s^3\beta _2+𝒪(a_s^4),\hfill \\ \hfill \beta _0& =\frac{11}{4}\frac{1}{6}n_f,\hfill \\ \hfill \beta _1& =\frac{51}{8}\frac{19}{24}n_f,\hfill \\ \hfill \beta _2& =\frac{2857}{128}\frac{5033}{1152}n_f+\frac{325}{3456}n_f^2,\hfill \\ \hfill \gamma ^m(a_s)& =a_s\gamma _0^ma_s^2\gamma _1^ma_s^3\gamma _2^m+𝒪(a_s^4),\hfill \\ \hfill \gamma _0^m& =1,\hfill \\ \hfill \gamma _1^m& =\frac{101}{24}\frac{5}{36}n_f,\hfill \\ \hfill \gamma _2^m& =\frac{1249}{64}\left(\frac{277}{216}+\frac{5}{6}\zeta _3\right)n_f\frac{35}{1296}n_f^2.\hfill \end{array}$$
(12)
Here, $`\zeta _n\zeta (n)`$ is Riemann’s $`\zeta `$-function ($`\zeta _31.20206`$). In order to evaluate $`\alpha _s(\mu _R)`$ from the initial value<sup>3</sup><sup>3</sup>3The numerical value of $`\alpha _s(M_Z)`$ has to be set in accordance with the parton sets that are used, see below. $`\alpha _s(M_Z)`$, $`\beta (a_s)`$ is expanded up to $`\alpha _s^{\mathrm{}}`$, with $`\mathrm{}=1`$ at LO, $`\mathrm{}=2`$ at NLO, and $`\mathrm{}=3`$ at NNLO. The resulting differential equation of Eq. (11) is solved numerically.

In order to evaluate $`\lambda _b(\mu )`$ from its initial value $`\lambda (\mu _0)`$, we combine the two RGEs of Eq. (11) to obtain

$$\begin{array}{c}\hfill \lambda _b(\mu )=\lambda _b(\mu _0)\frac{c(a_s(\mu ))}{c(a_s(\mu _0))},\end{array}$$
(13)
with

$$\begin{array}{cc}\hfill c(a)=a^{c_0}\{& 1+(c_1b_1c_0)a\hfill \\ & +\frac{1}{2}\left[(c_1b_1c_0)^2+c_2b_1c_1+b_1^2c_0b_2c_0\right]a^2\hfill \\ & +𝒪(a^3)\},c_i\frac{\gamma _i^m}{\beta _0},b_i\frac{\beta _i}{\beta _0}.\hfill \end{array}$$
(14)
$`\beta _i^m`$ and $`\gamma _i`$ have been defined in Eq. (12). Both $`a_s(\mu )`$ and $`a_s(\mu _0)`$ are calculated from $`\alpha _s(M_Z)`$ using the procedure described above. Working at LO (NLO, NNLO), we truncate the term in braces at order $`a^0`$ ($`a^1`$, $`a^2`$).

Convolution of the partonic cross section with the parton densities cancels the $`\mu _F`$ dependence up to higher orders and results in the physical hadronic cross section. The variation of the hadronic cross section with $`\mu _F`$ and $`\mu _R`$ is thus an indication of the size of higher order effects.

## IV Results

The analytic expressions for the partonic cross section are quite voluminous and will be deferred to the appendix. In this section, we study the behavior of the NNLO result with respect to variations of the input parameters, in particular the Higgs boson mass and the collider type (LHC and Tevatron). Special emphasis is placed on the variation of the results with the renormalization and factorization scale, from which we estimate the theoretical uncertainty of the prediction for Higgs boson production in $`b\overline{b}`$ annihilation.

Because the cross sections for the neutral Higgs bosons in $`b\overline{b}`$ annihilation differ only in the magnitudes of the Yukawa couplings (within our approximations), we will restrict our discussion to the production of a Standard Model Higgs boson. In the limit that supersymmetric partners are heavy, their virtual contributions are insignificant and the predictions for supersymmetric Higgs bosons can be obtained from the Standard Model values by rescaling them with the proper coupling constants (cf. Eq. (10)).

All the numerical results have been obtained using Martin-Roberts-Stirling-Thorne (MRST) parton distributions. In particular, we use the MRST2001 sets Martin:2001es  at LO $`\left[\alpha _s(M_Z)=0.130\right]`$ and NLO $`\left[\alpha _s(M_Z)=0.119\right]`$, and MRSTNNLO Martin:2002dr  at NNLO $`\left[\alpha _s(M_Z)=0.1155\right]`$.

In order to obtain an overall picture of the renormalization and factorization scale dependence of the cross section, we plot $`\sigma (pp(b\overline{b})H+X)`$ for $`\sqrt{s}=14`$ TeV as a function of the two parameters $`\mu _F`$ and $`\mu _R`$ in Fig. 5. The corresponding plot for the Tevatron, i.e., $`\sigma (p\overline{p}(b\overline{b})H+X)`$ for $`\sqrt{s}=1.96`$ TeV, is shown in Fig. 6. The Higgs boson mass is fixed to $`M_H=120`$ GeV. Subpanels (a), (b), and (c) show the LO, NLO, and NNLO prediction, respectively. Note the extremely large variation of the scales, by a factor of 10 above and below $`M_H`$. Apart from the region of large $`\mu _F`$ and small $`\mu _R`$, one observes a clear reduction of the scale dependence with increasing order of perturbation theory, both for $`\mu _F`$ and $`\mu _R`$. Notably, we find minimal radiative corrections and a particularly weak dependence on the renormalization and factorization scales in the vicinity of $`(\mu _R,\mu _F)=(M_\varphi ,0.25M_\varphi )(\overline{\mu }_R,\overline{\mu }_F)`$. This agrees with the observation of Ref. Maltoni:2003pn  that the proper factorization scale for this process should be around $`\mu _F=M_\varphi /4`$.

To illustrate this observation, we display separately the $`\mu _R`$\- and $`\mu _F`$-variation of the cross section at the LHC in Fig. 7 for $`M_H=120`$ GeV, and in Fig. 8 for $`M_H=300`$ GeV. In subpanels (a), the factorization scale is fixed to $`\mu _F=\overline{\mu }_F=0.25M_H`$, and the renormalization scale is varied within $`0.1\mu _R/M_H10`$. In subpanels (b), the renormalization scale is fixed to $`\mu _R=\overline{\mu }_R=M_H`$, and the factorization scale is varied within $`0.1\mu _F/M_H10`$. The reduction in the scale dependence with increasing order of perturbation theory is clearly visible. As opposed to the lower order curves which have a monotonic dependence on $`\mu _{R/F}`$ within the displayed range, the NNLO curves develop a maximum, so that it is possible to define a “point of least sensitivity” for them. In all cases, this falls nicely into a region where the radiative corrections are small. Note also that the central values for the NNLO curves are perfectly consistent between panels (a) and (b). These observations confirm that $`\overline{\mu }_F=0.25M_H`$ and $`\overline{\mu }_R=M_H`$ are indeed the appropriate scale choices for this process.

The corresponding curves for Run II at the Tevatron are shown in Fig. 9 (we only show results at the Tevatron for $`M_H=120`$ GeV). As opposed to the LHC, the reduction of the renormalization scale dependence with increasing order of perturbation theory is less drastic. One even observes a slight increase in the $`\mu _R`$ dependence between NLO and NNLO. However, the absolute variation is very small. The factorization scale dependence is quite similar to that observed for the LHC. Again, the central values for the NNLO curves of subpanels (a) and (b) coincide nicely. Note that the cross section at the Tevatron is typically about two orders of magnitude smaller than at the LHC.

Fig. 10(a) shows the LO, NLO, and NNLO predictions for the cross section $`\sigma (pp(b\overline{b})H+X)`$ at the LHC as a function of the Higgs boson mass. The two curves at each order correspond to two different choices of the factorization scale, $`\mu _F=0.1M_H`$ and $`\mu _F=0.7M_H`$. From subpanels (b) of Figs. 7 and 8 one can see that this roughly defines the maximal $`\mu _F`$-variation at NNLO between $`0.1M_H`$ and $`10M_H`$. Since the renormalization scale dependence is very weak \[cf. subpanels (a) of Figs. 7 and 8\], we fix $`\mu _R=\overline{\mu }_R=M_H`$. Taking the width of these bands as an indication of the theoretical uncertainty, we observe an improvement of the accuracy of the prediction for $`M_H=120`$ GeV from $`70\%`$ at LO to $`40\%`$ at NLO to $`15\%`$ at NNLO. At larger Higgs boson masses, the scale uncertainty is smaller, amounting to $`40\%`$ at LO, $`17\%`$ at NLO, and $`5\%`$ at NNLO for $`M_H=300`$ GeV.

The cross section for the Tevatron at $`\sqrt{s}=1.96`$ TeV center-of-mass energy is shown in Fig. 10(b). Here the renormalization scale dependence within the range $`0.1\mu _R/M_H10`$ at NNLO is larger than the factorization scale dependence (cf. Fig. 9). Nevertheless, we apply the same prescription as for the LHC and plot the LO, NLO, and NNLO cross section at $`(\mu _R,\mu _F)=(1,0.1)M_H`$ and $`(\mu _R,\mu _F)=(1,0.7)M_H`$. This is justified since $`\mu _R`$-variation on absolute scales is still very small, in particular if it is restricted to a more reasonable range of about a factor of five above and below $`M_H`$. We obtain a reduction of the scale uncertainty at the Tevatron for $`M_H=120`$ GeV from around $`60\%`$ at LO, to $`30\%`$ at NLO, to $`10\%`$ at NNLO.

## V Conclusions

We have computed the total cross section for Higgs boson production in $`b\overline{b}`$ fusion at NNLO in QCD. We have argued that the NNLO plays an exceptional role in this process, as it incorporates all subleading logarithms at order $`\alpha _s^2`$. The results are very stable with respect to changes of the renormalization and factorization scales. We find that the radiative corrections are particularly small at factorization scales of around $`\mu _F=M_\varphi /4`$, in agreement with the arguments of Refs. Plehn:2002vy ; Maltoni:2003pn .

We conclude that the inclusive cross section for Higgs boson production in bottom quark annihilation at hadron colliders is under good theoretical control.

#### Acknowledgments.

We would like to thank Scott Willenbrock, Fabio Maltoni, and Zack Sullivan for their encouragement and valuable comments. We are grateful to Tilman Plehn and Werner Vogelsang for enlightening discussions on bottom quark parton densities. R.V.H. would like to thank André Turcot for emphasizing the importance of this process, and Kostia Chetyrkin for useful comments. The work of W.B.K. was supported by the U. S. Department of Energy under Contract No. DE-AC02-98CH10886.

## Appendix A Partonic results

It is convenient to write the partonic cross section in the following way:

$$\begin{array}{cc}\hfill \widehat{\sigma }_{ij}(x)& =\sigma ^0\mathrm{\Delta }_{ij}(x),i,j\{b,\overline{b},g,q,\overline{q}\},\hfill \end{array}$$
(15)
where $`\widehat{\sigma }_{ij}`$ is the cross section for the process $`ij\varphi +X`$. $`i`$ and $`j`$ label the partons in the initial state, $`\varphi `$ means either a scalar or pseudoscalar Higgs boson, and $`X`$ denotes any number of quarks or gluons in the final state. Here and in what follows, $`q`$ denotes any of the light quarks $`u,d,s,c`$. The normalization factor, $`\sigma _0`$, is

$$\sigma _0=\frac{\pi }{12}\frac{\lambda _b^2}{M_\varphi ^2}.$$
(16)
The correction terms are written as a perturbative expansion:

$$\mathrm{\Delta }_{ij}(x)=\mathrm{\Delta }_{ij}^{(0)}(x)+\frac{\alpha _s}{\pi }\mathrm{\Delta }_{ij}^{(1)}(x)+\left(\frac{\alpha _s}{\pi }\right)^2\mathrm{\Delta }_{ij}^{(2)}(x)+𝒪(\alpha _s^3).$$
(17)
Explicit dependence on the number of active flavors $`n_f`$ appears only at NNLO. Because the terms are large and cumbersome, it is convenient to write

$$\mathrm{\Delta }_{ij}^{(2)}(x)=\mathrm{\Delta }_{ij}^{(2)A}(x)+n_f\mathrm{\Delta }_{ij}^{(2)F}(x).$$
(18)
All results will be presented for the scale choices $`\mu _F=\mu _R=M_\varphi `$. The corresponding expressions for general values of $`\mu _F`$ and $`\mu _R`$ can be reconstructed from renormalization scale invariance of the partonic, and factorization invariance of the hadronic cross section.<sup>4</sup><sup>4</sup>4 The analytic results including all scale dependences can also be obtained from the authors upon request.

### A.1 The $`b\overline{b}`$ subprocess

In the VFS approach, the tree-level $`b\overline{b}`$ annihilation term is the LO contribution. Thus, this is the only term for which $`\mathrm{\Delta }_{ij}^{(0)}(x)`$ does not vanish. The LO contribution to $`b\overline{b}\varphi +X`$ is

$$\mathrm{\Delta }_{b\overline{b}}^{(0)}(x)=\delta (1x).$$
(19)
The NLO contribution is

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{b\overline{b}}^{(1)}(x)=& \frac{48\zeta _2}{3}\delta (1x)+\frac{16}{3}𝒟_1(1x)\frac{16+8x+8x^2}{3}\mathrm{ln}\left(1x\right)+\frac{4x4x^2}{3}\frac{8}{3}\frac{\mathrm{ln}\left(x\right)}{1x}+\frac{8+4x+4x^2}{3}\mathrm{ln}\left(x\right),\hfill \end{array}$$
(20)
where $`𝒟_n(1x)\left[{\displaystyle \frac{\mathrm{ln}^n(1x)}{1x}}\right]_+`$, and $`\zeta _2\pi ^2/61.64493`$, $`\zeta _31.20206`$.

At NNLO, the contributions are ($`\zeta _4\pi ^4/901.08232`$)

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{b\overline{b}}^{(2)A}=& \frac{115+116\zeta _2156\zeta _319\zeta _4}{18}\delta \left(1x\right)\frac{404396\zeta _21146\zeta _3}{27}𝒟_0\left(1x\right)\hfill \\ & +\frac{204200\zeta _2}{9}𝒟_1\left(1x\right)\frac{44}{3}𝒟_2\left(1x\right)+\frac{128}{9}𝒟_3\left(1x\right)\hfill \\ & \frac{128+64x+64x^2}{9}\mathrm{ln}^3\left(1x\right)+\frac{140+40x+92x^28x^3}{9}\mathrm{ln}^2\left(1x\right)\hfill \\ & \frac{248}{9}\frac{\mathrm{ln}^2\left(1x\right)\mathrm{ln}\left(x\right)}{1x}+\frac{248+168x+168x^2}{9}\mathrm{ln}^2\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & \frac{604+138x+423x^2+44x^3}{27}\mathrm{ln}\left(1x\right)+\frac{200+100x+100x^2}{9}\zeta _2\mathrm{ln}\left(1x\right)\hfill \\ & +24\frac{\mathrm{ln}\left(1x\right)\mathrm{ln}\left(x\right)}{1x}\frac{216+110x+160x^224x^3}{9}\mathrm{ln}\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{148}{9}\frac{\mathrm{ln}\left(1x\right)\mathrm{ln}^2\left(x\right)}{1x}\frac{148+110x+110x^2}{9}\mathrm{ln}\left(1x\right)\mathrm{ln}^2\left(x\right)\hfill \\ & +\frac{20}{9}\frac{\mathrm{ln}\left(1x\right)Li_2\left(1x\right)}{1x}\frac{2078x78x^2}{9}\mathrm{ln}\left(1x\right)Li_2\left(1x\right)\hfill \\ & +\frac{4640+1017x+2958x^2+721x^3}{324}\frac{14015x+147x^28x^3}{9}\zeta _2\hfill \\ & \frac{382+191x+191x^2}{9}\zeta _3\frac{146}{9}\frac{\mathrm{ln}\left(x\right)}{1x}+\frac{164}{9}\frac{\zeta _2\mathrm{ln}\left(x\right)}{1x}\frac{23}{3}\frac{\mathrm{ln}^2\left(x\right)}{1x}\frac{44}{27}\frac{\mathrm{ln}^3\left(x\right)}{1x}\hfill \\ & +\frac{876+249x+444x^2+38x^3}{54}\mathrm{ln}\left(x\right)\frac{164+126x+126x^2}{9}\zeta _2\mathrm{ln}\left(x\right)\hfill \\ & +\frac{138+73x+115x^212x^3}{18}\mathrm{ln}^2\left(x\right)+\frac{44+21x+21x^2+4x^3}{27}\mathrm{ln}^3\left(x\right)\hfill \\ & +\frac{Li_2\left(1x\right)}{1x}+\frac{58}{9}\frac{Li_2\left(1x\right)\mathrm{ln}\left(x\right)}{1x}\frac{142}{9}\frac{Li_3\left(1x\right)}{1x}\frac{64}{9\left(1x\right)}Li_3\left(\frac{1x}{x}\right)\hfill \\ & +\frac{751x10x^2+10x^3}{9}Li_2\left(1x\right)\frac{58+88x+88x^2+2x^3}{9}Li_2\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & \frac{x2x^22x^3}{9}Li_2\left(1x^2\right)+\frac{x^3}{3}Li_2\left(1x^2\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{142+37x+37x^2+6x^3}{9}Li_3\left(1x\right)+\frac{64+94x+94x^26x^3}{9}Li_3\left(\frac{1x}{x}\right)\hfill \\ & \frac{7x^3}{18}Li_3\left(1x^2\right)\frac{x^3}{18}Li_3\left(\frac{1x^2}{x^2}\right)\frac{2x^3}{3}\left[Li_3\left(\frac{1x}{1+x}\right)Li_3\left(\frac{1x}{1+x}\right)\right],\hfill \end{array}$$
(21)
$$\begin{array}{cc}\hfill \mathrm{\Delta }_{b\overline{b}}^{(2)F}=& \frac{210\zeta _2+18\zeta _3}{27}\delta \left(1x\right)+\frac{5672\zeta _2}{81}𝒟_0\left(1x\right)\frac{40}{27}𝒟_1\left(1x\right)+\frac{8}{9}𝒟_2\left(1x\right)\hfill \\ & \frac{8+4x+4x^2}{9}\mathrm{ln}^2\left(1x\right)+\frac{40+8x+32x^2}{27}\mathrm{ln}\left(1x\right)\frac{16}{9}\frac{\mathrm{ln}\left(1x\right)\mathrm{ln}\left(x\right)}{1x}\hfill \\ & +\frac{16+8x+8x^2}{9}\mathrm{ln}\left(1x\right)\mathrm{ln}\left(x\right)+\frac{10}{9}\frac{\mathrm{ln}\left(x\right)}{1x}+\frac{2}{3}\frac{\mathrm{ln}^2\left(x\right)}{1x}\frac{2}{9}\frac{Li_2\left(1x\right)}{1x}\frac{56+x+55x^2}{81}\hfill \\ & +\frac{8+4x+4x^2}{9}\zeta _2\frac{10+3x+7x^2}{9}\mathrm{ln}\left(x\right)\frac{12+7x+7x^2}{18}\mathrm{ln}^2\left(x\right)+\frac{2}{9}Li_2\left(1x\right).\hfill \end{array}$$
(22)
Note that the $`𝒟_n`$ terms in this result could also be derived by other methods Kramer:1998iq ; Vogt:2000ci ; Harlander:2001is ; Catani:2001ic ; Kidonakis:2003tx ; Magnea:1991qg ; Contopanagos:1997nh .

### A.2 The $`bg`$ subprocess

The $`bg\varphi +X`$ subprocess first enters at NLO, where the contribution is

$$\begin{array}{cc}& \mathrm{\Delta }_{bg}^{(1)}=\mathrm{\Delta }_{\overline{b}g}^{(1)}=\frac{x2x^2+2x^3}{2}\mathrm{ln}\left(1x\right)\frac{3x10x^2+7x^3}{8}\frac{x2x^2+2x^3}{4}\mathrm{ln}\left(x\right).\hfill \end{array}$$
(23)
At NNLO, the contribution is

$$\begin{array}{cc}& \mathrm{\Delta }_{bg}^{(2)A}=\mathrm{\Delta }_{\overline{b}g}^{(2)A}=\hfill \\ & \frac{257x514x^2+514x^3}{144}\mathrm{ln}^3\left(1x\right)\hfill \\ & +\frac{1659x+272x^2237x^3}{16}\mathrm{ln}^2\left(1x\right)\frac{11x94x^2+62x^3}{8}\mathrm{ln}^2\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{16+28x731x^2+726x^3}{48}\mathrm{ln}\left(1x\right)\frac{35x70x^2+70x^3}{12}\zeta _2\mathrm{ln}\left(1x\right)\hfill \\ & +\frac{65x508x^2+774x^3}{24}\mathrm{ln}\left(1x\right)\mathrm{ln}\left(x\right)\frac{3x+174x^298x^3}{24}\mathrm{ln}\left(1x\right)\mathrm{ln}^2\left(x\right)\hfill \\ & +\frac{77x+134x^286x^3}{24}\mathrm{ln}\left(1x\right)Li_2\left(1x\right)+\frac{3x+6x^2+6x^3}{4}\mathrm{ln}\left(1x\right)Li_2\left(1x^2\right)\hfill \\ & \frac{208411x1350x^2+1781x^3}{288}\frac{1631x+176x^2169x^3}{16}\zeta _2+\frac{161x322x^2+322x^3}{48}\zeta _3\hfill \\ & +\frac{32x+536x^2993x^3}{48}\mathrm{ln}\left(x\right)+\frac{7x230x^2+134x^3}{24}\zeta _2\mathrm{ln}\left(x\right)\hfill \\ & \frac{47x604x^2+1028x^3}{96}\mathrm{ln}^2\left(x\right)\frac{35x+146x^2+12x^3}{144}\mathrm{ln}^3\left(x\right)\hfill \\ & +\frac{4843x+152x^2+236x^3}{24}Li_2\left(1x\right)\frac{10x+34x^217x^3}{6}Li_2\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & \frac{19x+40x^2+21x^3}{48}Li_2\left(1x^2\right)\frac{9x+18x^2+10x^3}{12}Li_2\left(1x^2\right)\mathrm{ln}\left(x\right)\hfill \\ & \frac{14x+8x^2+9x^3}{3}Li_3\left(1x\right)+\frac{13x+118x^218x^3}{12}Li_3\left(\frac{1x}{x}\right)\hfill \\ & +\frac{9x+18x^214x^3}{48}Li_3\left(1x^2\right)+\frac{9x+18x^214x^3}{48}Li_3\left(\frac{1x^2}{x^2}\right)\hfill \\ & +\frac{3x+6x^2+6x^3}{2}\left[Li_3\left(\frac{1x}{1+x}\right)Li_3\left(\frac{1x}{1+x}\right)\right],\hfill \\ \hfill \mathrm{\Delta }_{bg}^{(2)F}=& \mathrm{\Delta }_{\overline{b}g}^{(2)F}=0.\hfill \end{array}$$
(24)
### A.3 The $`bq`$ subprocess

All remaining components contribute only at NNLO and beyond. The contribution to $`bq\varphi +X`$, where $`q`$ is a light ($`u`$,$`d`$,$`s`$ or $`c`$) quark, is:

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{bq}^{(2)A}=& \mathrm{\Delta }_{\overline{b}q}^{(2)A}=\mathrm{\Delta }_{b\overline{q}}^{(2)A}=\mathrm{\Delta }_{\overline{b}\overline{q}}^{(2)A}=\frac{4+3x3x^24x^3}{9}\mathrm{ln}^2\left(1x\right)\hfill \\ & +\frac{2x+2x^2}{3}\mathrm{ln}^2\left(1x\right)\mathrm{ln}\left(x\right)+\frac{457x+75x^222x^3}{27}\mathrm{ln}\left(1x\right)\hfill \\ & \frac{x4x^24x^3}{3}\mathrm{ln}\left(1x\right)\mathrm{ln}\left(x\right)\frac{2x+2x^2}{3}\mathrm{ln}\left(1x\right)\mathrm{ln}^2\left(x\right)\hfill \\ & +\frac{4x+4x^2}{3}\mathrm{ln}\left(1x\right)Li_2\left(1x\right)\frac{208915x+1410x^2703x^3}{648}\hfill \\ & \frac{4+3x3x^24x^3}{9}\zeta _2+\frac{93x264x^2+20x^3}{108}\mathrm{ln}\left(x\right)\frac{2x+2x^2}{3}\zeta _2\mathrm{ln}\left(x\right)\hfill \\ & \frac{3x+15x^2+40x^3}{72}\mathrm{ln}^2\left(x\right)+\frac{x+x^2}{36}\mathrm{ln}^3\left(x\right)+\frac{163x+21x^2+8x^3}{18}Li_2\left(1x\right)\hfill \\ & \left(x+x^2\right)Li_2\left(1x\right)\mathrm{ln}\left(x\right)\frac{2x+2x^2}{3}Li_3\left(1x\right)+\frac{2x+2x^2}{3}Li_3\left(\frac{1x}{x}\right),\hfill \\ \hfill \mathrm{\Delta }_{bq}^{(2)F}=& \mathrm{\Delta }_{\overline{b}q}^{(2)F}=\mathrm{\Delta }_{b\overline{q}}^{(2)F}=\mathrm{\Delta }_{\overline{b}\overline{q}}^{(2)F}=0.\hfill \end{array}$$
(25)
### A.4 The $`gg`$ subprocess

The contribution to $`gg\varphi +X`$ is:

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{gg}^{(2)A}=& \left(x+2x^23x^3\right)\mathrm{ln}^2\left(1x\right)\frac{x+4x^2+4x^3}{2}\mathrm{ln}^2\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{23x+52x^275x^3}{8}\mathrm{ln}\left(1x\right)+\frac{5x+16x^24x^3}{4}\mathrm{ln}\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{x+4x^2+4x^3}{4}\mathrm{ln}\left(1x\right)\mathrm{ln}^2\left(x\right)\left(x+4x^2+4x^3\right)\mathrm{ln}\left(1x\right)Li_2\left(1x\right)\hfill \\ & \frac{163x+1528x^21691x^3}{128}+\left(x+2x^23x^3\right)\zeta _2\hfill \\ & \frac{54x+312x^2223x^3}{64}\mathrm{ln}\left(x\right)+\frac{x+4x^2+4x^3}{2}\zeta _2\mathrm{ln}\left(x\right)\hfill \\ & \frac{16x+111x^243x^3}{64}\mathrm{ln}^2\left(x\right)+\frac{7x+25x^2+34x^3}{48}\mathrm{ln}^3\left(x\right)\hfill \\ & \frac{4x15x^262x^3}{16}Li_2\left(1x\right)+\frac{11x+44x^2+30x^3}{16}Li_2\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{x^26x^3}{32}Li_2\left(1x^2\right)+\frac{3x+6x^2+38x^3}{64}Li_2\left(1x^2\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{x+3x^2+18x^3}{8}Li_3\left(1x\right)\frac{15x+60x^2+30x^3}{16}Li_3\left(\frac{1x}{x}\right)\hfill \\ & \frac{5x+10x^2+74x^3}{128}Li_3\left(1x^2\right)\frac{3x+6x^2+70x^3}{128}Li_3\left(\frac{1x^2}{x^2}\right)\hfill \\ & \frac{x+2x^2+2x^3}{32}\left[Li_3\left(\frac{1x}{1+x}\right)Li_3\left(\frac{1x}{1+x}\right)\right],\hfill \\ \hfill \mathrm{\Delta }_{gg}^{(2)F}=& 0.\hfill \end{array}$$
(26)
### A.5 The $`bb`$ subprocess

The contribution to $`bb\varphi +X`$ is:

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{bb}^{(2)A}=& \mathrm{\Delta }_{\overline{b}\overline{b}}^{(2)A}=\frac{8+6x6x^28x^3}{9}\mathrm{ln}^2\left(1x\right)+\frac{4x+4x^2}{3}\mathrm{ln}^2\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{8138x+174x^244x^3}{27}\mathrm{ln}\left(1x\right)\frac{10x20x^224x^3}{9}\mathrm{ln}\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{4}{9}\frac{\mathrm{ln}\left(1x\right)\mathrm{ln}^2\left(x\right)}{1+x}\frac{16}{9}\frac{\mathrm{ln}\left(1x\right)Li_2\left(1x\right)}{1+x}+\frac{8}{9}\frac{\mathrm{ln}\left(1x\right)Li_2\left(1x^2\right)}{1+x}\hfill \\ & \frac{4+10x+14x^2}{9}\mathrm{ln}\left(1x\right)\mathrm{ln}^2\left(x\right)+\frac{16+16x+32x^2}{9}\mathrm{ln}\left(1x\right)Li_2\left(1x\right)\hfill \\ & \frac{84x+4x^2}{9}\mathrm{ln}\left(1x\right)Li_2\left(1x^2\right)\frac{52357x+510x^2205x^3}{81}\hfill \\ & \frac{8+6x6x^28x^3}{9}\zeta _2+\frac{117x279x^2+20x^3}{54}\mathrm{ln}\left(x\right)\frac{4x+4x^2}{3}\zeta _2\mathrm{ln}(x)\hfill \\ & \frac{x+11x^2+34x^3}{36}\mathrm{ln}^2(x)+\frac{2812x+17x^2}{54}\mathrm{ln}^3\left(x\right)\frac{14}{27}\frac{\mathrm{ln}^3\left(x\right)}{1+x}\hfill \\ & +\frac{4}{3}\frac{Li_2\left(1x\right)\mathrm{ln}\left(x\right)}{1+x}\frac{8}{9}\frac{Li_2\left(1x^2\right)\mathrm{ln}\left(x\right)}{1+x}\frac{4}{3}\frac{Li_3\left(1x\right)}{1+x}+\frac{4}{3\left(1+x\right)}Li_3\left(\frac{1x}{x}\right)\hfill \\ & +\frac{1}{3}\frac{Li_3\left(1x^2\right)}{1+x}+\frac{1}{9\left(1+x\right)}Li_3\left(\frac{1x^2}{x^2}\right)+\frac{20}{9\left(1+x\right)}\left[Li_3\left(\frac{1x}{1+x}\right)Li_3\left(\frac{1x}{1+x}\right)\right]\hfill \\ & +\frac{1612x+16x^2+11x^3}{9}Li_2\left(1x\right)\frac{12+8x+29x^2}{9}Li_2\left(1x\right)\mathrm{ln}\left(x\right)\hfill \\ & +\frac{x}{9}Li_2\left(1x^2\right)+\frac{1611x+11x^2}{18}Li_2\left(1x^2\right)\mathrm{ln}\left(x\right)+\frac{1227x+5x^2}{9}Li_3\left(1x\right)\hfill \\ & \frac{1216x9x^22x^3}{9}Li_3\left(\frac{1x}{x}\right)\frac{1213x+13x^2}{36}Li_3\left(1x^2\right)\hfill \\ & \frac{43x+3x^2}{36}Li_3\left(\frac{1x^2}{x^2}\right)\frac{2013x+13x^2}{9}\left[Li_3\left(\frac{1x}{1+x}\right)Li_3\left(\frac{1x}{1+x}\right)\right],\hfill \\ \hfill \mathrm{\Delta }_{bb}^{(2)F}=& \mathrm{\Delta }_{\overline{b}\overline{b}}^{(2)F}=0.\hfill \end{array}$$
(27)
### A.6 The $`q\overline{q}`$ subprocess

The contribution to $`q\overline{q}\varphi +X`$ is ($`q\{u,d,s,c\}`$):

$$\begin{array}{cc}\hfill \mathrm{\Delta }_{q\overline{q}}^{(2)A}=& \frac{2x8x^2+6x^3}{9}\frac{x2x^23x^3}{9}\mathrm{ln}\left(x\right)+\frac{x^3}{9}\mathrm{ln}^2\left(x\right)\hfill \\ & \frac{4x^3}{9}Li_2\left(1x\right)+\frac{2x^3}{9}Li_2\left(1x^2\right),\hfill \\ \hfill \mathrm{\Delta }_{q\overline{q}}^{(2)F}=& 0.\hfill \end{array}$$
(28)


format: text

text:
# Dark Matter and Dark Energy: Summary and Future Directions

## 1 The Density Budget of the Universe

It is convenient to express the mean densities $`\rho _i`$ of various quantities in the Universe in terms of their fractions relative to the critical density: $`\mathrm{\Omega }_i\rho _i/\rho _{crit}`$. The theory of cosmological inflation strongly suggests that the the total density should be very close to the critical one: $`\mathrm{\Omega }_{tot}1`$, and this is supported by the available data on the cosmic microwave background radiation (CMB) (Bond 2003). The fluctuations observed in the CMB at a level $`10^5`$ in amplitude exhibit a peak at a partial wave $`\mathrm{}200`$, as would be produced by acoustic oscillations in a flat Universe with $`\mathrm{\Omega }_{tot}1`$. At lower partial waves, $`\mathrm{}200`$, the CMB fluctuations are believed to be dominated by the Sachs-Wolfe effect due to the gravitational potential, and more acoustic oscillations are expected at larger $`\mathrm{}>200`$, whose relative heights depend on the baryon density $`\mathrm{\Omega }_b`$. At even larger $`\mathrm{}>1000`$, these oscillations should be progressively damped away.

Fig. 1 compares measurements of CMB fluctuations made before WMAP (Bond 2003) with the WMAP data themselves (Bennett et al. 2003; Hinshaw et al. 2003), that were released shortly after this Meeting. The position of the first acoustic peak indeed corresponds to a flat Universe with $`\mathrm{\Omega }_{tot}1`$: in particular, now WMAP finds $`\mathrm{\Omega }_{tot}=1.02\pm 0.02`$ (Spergel et al. 2003), and two more acoustic peaks are established with high significance, providing a new determination of $`\mathrm{\Omega }_bh^2=0.0224\pm 0.0009`$, where $`h0.7`$ is the present Hubble expansion rate $`H`$, measured in units of $`100`$km/s/Mpc. The likelihood functions for various cosmological parameters are shown in Fig. 2. Remarkably, there is excellent consistency between the estimate of the present-day Hubble constant $`H72`$km/s/Mpc from WMAP (Spergel et al. 2003) with that inferred from the local distance ladder based, e.g., on Cepheid variables.

As seen in Fig 3, the combination of CMB data with those on high-redshift Type-Ia supernovae (Perlmutter 2003; Perlmutter & Schmidt 2003) and on large-scale structure (Peacock 2003$`a`$, $`b`$) favour strongly a flat Universe with about 30 % of (mainly dark) matter and 70 % of vacuum (dark) energy. Type-Ia supernovae probe the geometry of the Universe at redshifts $`z<1`$. They disagree with a flat $`\mathrm{\Omega }_{tot}=1`$ Universe that has no vacuum energy, and also with an open $`\mathrm{\Omega }_m0.3`$ Universe (Perlmutter 2003; Perlmutter & Schmidt 2003). They appear to be adequate standard candles, and two observed supernovae with $`z>1`$ argue strongly against dust or evolution effects that would be sufficient to cloud their geometrical interpretation. The supernovae indicate that the expansion of the Universe is currently accelerating, though it had been decelerating when $`z`$ was $`>1`$. There are good prospects for improving substantially the accuracy of the supernova data, by a combination of continued ground-based and subsequent space observations using the SNAP satellite project (Perlmutter 2003; Perlmutter & Schmidt 2003).

It is impressive that the baryon density inferred from WMAP data (Spergel et al. 2003) is in good agreement with the value calculated previously on the basis of Big-Bang nucleosynthesis (BBN), which depends on completely different (nuclear) physics. Fig. 4 compares the abundances of light elements calculated using the WMAP value of $`\mathrm{\Omega }_bh^2`$ with those inferred from astrophysical data (Cyburt et al. 2003). Depending on the astrophysical assumptions that are made in extracting the light-element abundances from astrophysical data, there is respectable overlap.

As we heard at this Meeting, several pillars of inflation theory have now been verified by WMAP and other CMB data (Bond 2003): the Sachs-Wolfe effect due to fluctuations in the large-scale gravitational potential were first seen by the COBE satellite, the first acoustic peak was seen in the CMB spectrum at $`\mathrm{}210`$ and this has been followed by two more peaks and the intervening dips, the damping tail of the fluctuation spectrum expected at $`\mathrm{}>1000`$ has been seen, polarization has been observed, and the primary anisotropies are predominantly Gaussian. WMAP has, additionally, measured the thickness of the last scattering surface and observed the reionization of the Universe when $`z20`$ by the first generation of stars (Kogut et al. 2003). Remaining to be established are secondary anisotropies, due, e.g., to the Sunyaev-Zeldovich effect, weak lensing and inhomogeneous reionization, and tensor perturbations induced by gravity waves.

As we also heard at this meeting, the values of $`\mathrm{\Omega }_{CDM}`$ inferred from X-ray studies of gas in rich clusters using the Chandra satellite (Rees 2003), which indicate $`\mathrm{\Omega }_{CDM}=0.325\pm 0.34`$, gravitational lensing (Schneider 2003) and data on large-scale structure, e.g., from the 2dF galaxy redshift survey (Peacock 2003$`a`$, $`b`$), are very consistent with that inferred by combining CMB and supernova data. The WMAP data confirm this concordance with higher precision: $`\mathrm{\Omega }_{CDM}h^2=0.111\pm 0.009`$ (Spergel et al. 2003).

The 2dF galaxy survey has examined two wedges through the Universe. Significant structures are seen at low redshifts, which die away at larger redshifts where the Universe becomes more homogeneous and isotropic. The perturbation power spectrum at these large scales matches nicely with that seen in the CMB data, whilst the structures seen at small scales would not be present in a baryon-dominated Universe, or one with a significant fraction of hot dark matter. Indeed, the 2dF data were used to infer an upper limit on the sum of the neutrino masses of 1.8 eV (Elgaroy et al. 2002), which has recently been improved using WMAP data (Spergel et al. 2003) to

$$\mathrm{\Sigma }_{\nu _i}m_{\nu _i}<0.7\mathrm{eV},$$
(1)
as seen in Fig. 5. This impressive upper limit is substantially better than even the most stringent direct laboratory upper limit on an individual neutrino mass, as discussed in the next Section. Thw WMAP data also provide (Crotty et al. 2003) a new limit on the effective number of light neutrino species, beyond the three within the Standard Model:

$$1.5<\mathrm{\Delta }N_\nu ^{eff}<\mathrm{\hspace{0.33em}4.2}.$$
(2)
This limit is not as stringent as that from LEP, but applies to additional light degrees of freedom that might not be produced in $`Z`$ decay.

## 2 What is it?

As discussed here by Kolb (2003), particle candidates for dark matter range from the axion with a mass $`>10^{15}`$ GeV (van Bibber 2003) to cryptons with masses $`<10^{+15}`$ GeV (Ellis et al. 1990; Benakli et al. 1999), via neutrinos with masses $`<10^{10}`$ GeV, the gravitino and the lightest supersymmetric particle with a mass $`>10^2`$ GeV (Ellis et al. 1984; Goldberg 1983). In recent years, there has been considerable experimental progress in understanding neutrino masses, so I start with them, even though cosmology now disfavours the hot dark matter they would provide (Spergel et al. 2003). All the others are candidates for cold dark matter, except for the gravitino, which might constitute warm dark matter, another possibility now disfavoured by the WMAP evidence for reionization when $`z20`$ (Kogut et al. 2003).

### 2.1 Neutrinos

Particle theorists expect particles to have masses that vanish exactly only if they are protected by some unbroken gauge symmetry, much as the photon is massless because of the U(1) gauge symmetry of electromagnetism, that is associated with the conservation of electric charge. There is no corresponding exact gauge symmetry to protect lepton number, so we expect it to be violated and neutrinos to acquire masses. This is indeed the accepted interpretation of the observed oscillations between different types of neutrinos, which are made possible by mixing into non-degenerate mass eigenstates (Wark 2003); Pakvasa & Valle 2003).

Neutrino masses could arise even within the Standard Model of particle physics, without adding any new particles, at the expense of introducing a interaction between two neutrino fields and two Higgs fields (Barbieri et al. 1980):

$$\frac{1}{M}\nu H\nu Hm_\nu =\frac{0|H0^2}{M}.$$
(3)
However, such an interaction would be non-renormalizable, and therefore is not thought to be fundamental. The (presumably large) mass scale $`M`$ appearing the denominator of (3) is generally thought to originate from the exchange of some massive fermionic particle that mixes with the light neutrino (Gell-Mann et al. 1979; Yanagida 1979; Mohapatra & Senjanovic 1980):

$`(\nu _L,N)\left(\begin{array}{cc}0& M_D\\ M_D^T& M\end{array}\right)\left(\begin{array}{c}\nu _L\\ N\end{array}\right),`$ (8)
Diagonalization of this matrix naturally yields small neutrino masses, since we expect that the Dirac mass term $`m_D`$ is of the same order as quark and lepton masses, and $`Mm_W`$.

We have the following direct experimental upper limits on neutrino masses. From measurements of the end-point in Tritium $`\beta `$ decay, we know that (Weinheimer et al. 1999; Lobashov et al. 1999):

$$m_{\nu _e}<\mathrm{\hspace{0.33em}2.5}\mathrm{eV},$$
(9)
and there are prospects to improve this limit down to about 0.5 eV with the proposed KATRIN experiment (Osipowicz et al. 2001). From measurements of $`\pi \mu \nu `$ decay, we know that (Hagiwara et al. 2002):

$$m_{\nu _\mu }<\mathrm{\hspace{0.33em}190}\mathrm{KeV},$$
(10)
and there are prospects to improve this limit by a factor $`20`$. From measurements of $`\tau n\pi \nu `$ decay, we know that (Hagiwara et al. 2002):

$$m_{\nu _\tau }<\mathrm{\hspace{0.33em}18.2}\mathrm{MeV},$$
(11)
and there are prospects to improve this limit to $`5`$ MeV.

However, the most stringent laboratory limit on neutrino masses may come from searches for neutrinoless double-$`\beta `$ decay, which constrain the sum of the neutrinos’ Majorana masses weighted by their couplings to electrons (Klapdor-Kleingrothaus et al. 2001):

$$m_\nu _e|\mathrm{\Sigma }_{\nu _i}m_{\nu _i}U_{ei}^2|<0.35\mathrm{eV}$$
(12)
and there are prospects to improve this limit to $`0.01`$ eV in a future round of experiments. The impact of the limit (12) in relation to the cosmological upper limit (1) is discussed below, after we have gathered further experimental input from neutrino-oscillation experiments.

The neutrino mass matrix (8) should be regarded also as a matrix in flavour space. When it is diagonalized, the neutrino mass eigenstates will not, in general, coincide with the flavour eigenstates that partner the mass eigenstates of the charged leptons. The mixing matrix between them (Maki et al. 1962) may be written in the form

$`V=\left(\begin{array}{ccc}c_{12}& s_{12}& 0\\ s_{12}& c_{12}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& c_{23}& s_{23}\\ 0& s_{23}& c_{23}\end{array}\right)\left(\begin{array}{ccc}c_{13}& 0& s_{13}\\ 0& 1& 0\\ s_{13}e^{i\delta }& 0& c_{13}e^{i\delta }\end{array}\right).`$ (22)
where the symbols $`c,s_{ij}`$ denote the standard trigonometric functions of the three real ‘Euler’ mixing angles $`\theta _{12,23,31}`$, and $`\delta `$ is a CP-violating phase that can in principle also be observed in neutrino-oscillation experiments (De Rújula et al. 1999). Additionally, there are two CP-violating phases $`\varphi _{1,2}`$ that appear in the double-$`\beta `$ observable (12), but do not affect neutrino oscillations.

The pioneering Super-Kamiokande and other experiments have shown that atmospheric neutrinos oscillate, with the following difference in squared masses and mixing angle (Fukuda et al. 1998):

$$\delta m^2\mathrm{\hspace{0.33em}2.4}\times 10^3\mathrm{eV}^2,\mathrm{sin}^22\theta _{23}1.0,$$
(23)
which is very consistent with the K2K reactor neutrino experiment (Ahn et al. 2002), as seen in the left panel of Fig. 6. A flurry of recent solar neutrino experiments, most notably SNO (Ahmad et al. 2002$`a`$, $`b`$), have established beyond any doubt that they also oscillate, with

$$\delta m^2\mathrm{\hspace{0.33em}6}\times 10^5\mathrm{eV}^2,\mathrm{tan}^2\theta _{31}0.5.$$
(24)
Most recently, the KamLAND experiment has reported a deficit of electron antineutrinos from nuclear power reactors, leading to a very similar set of preferred parameters, as seen in the right panel of Fig. 6 (Eguchi et al. 2003).

Using the range of $`\theta _{12}`$ allowed by the solar and KamLAND data, one can establish a correlation between the relic neutrino density $`\mathrm{\Omega }_\nu h^2`$ and the neutrinoless double-$`\beta `$ decay observable $`m_\nu _e`$, as seen in Fig. 7 (Minakata & Sugiyama 2002). Pre-WMAP, the experimental limit on $`m_\nu _e`$ could be used to set the bound (Minakata & Sugiyama 2002)

$$10^3<\mathrm{\Omega }_\nu h^2<\mathrm{\hspace{0.33em}10}^1.$$
(25)
Alternatively, now that WMAP has set a tighter upper bound $`\mathrm{\Omega }_\nu h^2<0.0076`$ (1), one can use this correlation to set an upper bound:

$$<m_\nu >_e<\mathrm{\hspace{0.33em}0.1}\mathrm{eV},$$
(26)
which is difficult to reconcile with the signal reported in (Klapdor-Kleingrothaus et al. 2002).

The ‘Holy Grail’ of neutrino physics is CP violation in neutrino oscillations, which would manifest itself as a difference between the oscillation probabilities for neutrinos and antineutrinos (De Rújula et al. 1999):

$`P\left(\nu _e\nu _\mu \right)P\left(\overline{\nu }_e\overline{\nu }_\mu \right)`$ $`=\mathrm{\hspace{0.33em}16}s_{12}c_{12}s_{13}c_{13}^2s_{23}c_{23}\mathrm{sin}\delta `$
$`\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Delta }m_{12}^2}{4E}}L\right)\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Delta }m_{13}^2}{4E}}L\right)\mathrm{sin}\left({\displaystyle \frac{\mathrm{\Delta }m_{23}^2}{4E}}L\right),`$
For this to be observable, $`\mathrm{\Delta }m_{12}`$ and $`\theta _{12}`$ have to be large, as SNO and KamLAND have shown to be the case, and also $`\theta _{13}`$ has to be large enough - which remains to be seen.

In fact, even the minimal seesaw model contains many additional parameters that are not observable in neutrino oscillations (Casas & Ibarra 2001). In addition to the three masses of the charged leptons, there are three light-neutrino masses, three light-neutrino mixing angles and three CP-violating phases in the light-neutrino sector: the oscillation phase $`\delta `$ and the two Majorana phases that are relevant to neutrinoless double-$`\beta `$ decay experiments. As well, there are three heavy singlet-neutrino masses, three more mixing angles and three more CP-violating phases that become observable in the heavy-neutrino sector, making a total of 18 parameters. Out of all these parameters, so far just four light-neutrino parameters are known, two differences in masses squared and two real mixing angles.

As discussed later, it is often thought that there may be a connection between CP violation in the neutrino sector and the baryon density in the Universe, via leptogenesis. Unfortunately, this connection is somewhat indirect, since the three CP-violating phases measurable in the light-neutrino sector to not contribute to leptogenesis, which is controlled by the other phases that are not observable directly at low energies (Ellis & Raidal 2002). However, if the seesaw model is combined with supersymmetry, these extra phases contribute to the renormalization of soft superymmetry-breaking parameters at low energies, and hence may have indirect observable effects (Ellis et al. 2002$`a`$, $`b`$).

### 2.2 Problems with Cold Dark Matter?

A recurring question is whether the cold dark matter paradigm for structure formation is compatible with all the observational data (Navarro 2003; Abadi et al. 2002). One of the issues is the mass profile at the core of a galactic halo. As we heard at this meeting, density profiles clearly differ from naive power laws, and are much shallower than simple isothermal models. However, the haloes of $`\gamma <2.5`$ seem to be in reasonable agreement with CDM sumulations, though there are still problems with $`\gamma >2.5`$ galaxies. The theoretical predictions are not yet conclusive, though, with questions such as triaxiality, departures from equilibrium and time dependence remaining to be resolved. Another issue is halo substructure and the abundance of Milky Way satellites. However, the latest news is that there is apparently better agreement between the number of Milky Way satellites and the number of massive halo substructures predicted in CDM simulations. It has been argued that a Milky Way-like stellar disk would be thickened if there were substructures in the halo, but recent simulations do not display any significant such effect (Navarro 2003; Abadi et al. 2002).

Thus, the latest advice from the simulators seems to be that there is no showstopper for the CDM paradigm, so let us examine some of the candidates for the CDM.

### 2.3 Axions

As we heard here from van Bibber (van Bibber 2003), axions were invented in order to conserve CP in the strong interactions: in the picturesque analogy of Sikivie (Sikivie 1996), to explain why the strong-interaction pool table is apparently horizontal. Axion-like particles may be characterized by a two-dimensional parameter space, consisting of the axion mass $`m_a`$ and its coupling to pairs of photons, $`g_{a\gamma }`$. Many areas of this parameter space are excluded by laser experiments, telescopes, searches for solar axions currently being extended by the CAST experiment at CERN (Irastorza et al. 2002), astrophysical constraints (Hagiwara et al. 2002) and searches for halo axions using microwave cavities (Asztalos et al. 2001), as seen in Fig. 8. These and searches using Rydberg atoms (Yamamoto et al. 2001) have the best chances of excluding ragions of parameter space where the axion might constitute cold dark matter.

### 2.4  Supersymmetric Dark Matter

The appearance of supersymmetry at the TeV scale was originally motivated by the hierarchy problem (Maiani 1979; ’t Hooft 1979; Witten 1981) - why is $`m_Wm_P10^{19}`$ GeV, the only candidate we have for a fundamental mass scale in physics, or alternatively why is the Coulomb potential in an atom so much larger than the Newton potential? The former is $`e^2=𝒪(1)`$, whereas the latter is $`G_Nm^2m^2/m_P^2`$, where $`m`$ is a typical particle mass scale. It is not sufficient simply to set particle masses such as $`m_Wm_P`$, since quantum corrections will increase them again. Many one-loop quantum corrections each yield

$$\delta m_W^2𝒪\left(\frac{\alpha }{\pi }\right)\mathrm{\Lambda }^2,$$
(28)
where $`\mathrm{\Lambda }`$ is an effective cut-off, representing the scale at which the Standard Model ceases to be valid, and new physics appears. If $`\mathrm{\Lambda }m_P`$ or the GUT scale $`10^{16}`$ GeV, the ‘small’ correction (28) will be much larger than the physical value of $`m_W^2`$. It would require ‘unnatural’ fine-tuning to choose the bare value of $`m_W^2`$ to be almost equal and opposite to the ‘small’ correction (28), so that their combination happens to have the right magnitude. Alternatively, supersymmetry introduces an effective cut-off $`\mathrm{\Lambda }`$ by postulating equal numbers of bosons and fermions with identical couplings, in which case the corrections (28) cancel among themselves, leaving

$$\delta m_W^2𝒪\left(\frac{\alpha }{\pi }\right)(m_B^2m_F^2),$$
(29)
which is $`<m_W^2`$ if

$$|m_B^2m_F^2|<𝒪(1)\mathrm{TeV}^2,$$
(30)
i.e., if supersymmetry appears at relatively low energy.

This naturalness argument for low-energy supersymmetry is supported by several pieces of indirect empirical evidence. One is that the strengths of the different gauge interactions measured at low energies, particularly at the LEP accelerator, are consistent with unification at high energies if there are low-mass supersymmetric particles (Ellis et al. 1991$`b`$; Amaldi et al. 1991; Langacker & Luo 1991; Giunti et al. 1991). A second is that LEP and other precision low-energy data are fitted well by the Standard Model if there is a relatively light Higgs boson with mass $`<200`$ GeV, very consistent with the range $`m_h<130`$ GeV predicted by supersymmetry (Okada et al. 1991; Ellis et al. 1991$`a`$; Haber & Hempfling 1991). A third is that the minimal supersymmetric extension of the Standard Model (MSSM) contains a good candidate $`\chi `$ for cold dark matter, which has a suitable relic density if the supersymmetric mass scale $`<1`$ TeV (Ellis et al. 1984; Goldberg 1983). A fourth may be provided by the anomalous magnetic moment of the muon, $`g_\mu 2`$ (Brown et al. 2001; Bennett et al. 2002), if its experimental value deviates significantly from the Standard Model prediction.

When considering the experimental, cosmological and theoretical constraints on the MSSM, it is common to assume that all the unseen spin-0 supersymmetric particles have some universal mass $`m_0`$ at some GUT input scale, and similarly for the unseen fermion masses $`m_{1/2}`$. These two parameters of this constrained MSSM (CMSSM) are restricted by the absences of supersymmetric particles at LEP: $`m_{\chi ^\pm }>103`$ GeV, $`m_{\stackrel{~}{e}}>99`$ GeV, and at the Fermilab Tevatron collider. They are also restricted indirectly by the absence of a Higgs boson at LEP: $`m_h>114.4`$ GeV, and by the fact that $`bs\gamma `$ decay is consistent with the Standard Model, and potentially by the BNL measurement of $`g_\mu 2`$, as seen in Fig. 9 (Ellis et al. 2003$`a`$; Lahanas & Naopoulos 2003).

As shown there, the CMSSM parameter space is also restricted by cosmological bounds on the amount of cold dark matter, $`\mathrm{\Omega }_{CDM}h^2`$. Since $`\rho _\chi =m_\chi n_\chi `$, and the relic number density $`n_\chi 1/\sigma _{ann}(\chi \chi \mathrm{})`$ where $`\sigma _{ann}(\chi \chi \mathrm{})1/m^2`$, the relic density generically increases with increasing sparticle masses. For some time, the conservative upper limit on $`\mathrm{\Omega }_{CDM}h^2`$ has been $`<0.3`$ (Lahanas et al. 2000, 2001a, 2001b; Barger & Kao 2001; Arnowitt & Dutta 2002), but the recent WMAP data allow this to be reduced to $`<0.129`$ at the 2-$`\sigma `$ level. As seen in Fig. 9, this improved upper limit significantly improves the cosmological upper limit on the sparticle mass scale (Ellis et al. 2003$`a`$, Lahanas & Nanopoulos 2003). If there are other important components of the cold dark matter, the CMSSM parameters could lie below the dark (blue) strips in Fig. 9.

In order to facilitate discussion of the physics reaches of different accelerators and strategies for detecting dark matter, it is convenient to focus on a limited number of benchmark scenarios that illustrate the various different supersymmetric possibilities (Battaglia et al. 2001). In many of these scenarios, supersymmetry would be very easy to observe at the LHC, and many different types of supersymmetric particle might be discovered. However, searches for astrophysical dark matter may be competitive for some scenarios (Ellis et al. 2001).

One strategy is to look for relic annihilations out in the galactic halo, which might produce detectable antiprotons or positrons in the cosmic rays (Silk & Srednicki 1984). As discussed here by Carr (2003), both the PAMELA (Pearce et al. 2002) and AMS (Aguilar et al. 2002) space experiments will be looking for these signals, though the rates are not very promising in the benchmark scenarios we studied (Ellis et al. 2001). Alternatively, one might look for annihilations in the core of our galaxy, which might produce detectable gamma rays. As seen in the left panel of Fig. 10, this may be possible in certain benchmark scenarios (Ellis et al. 2001), although the rate is rather uncertain because of the unknown enhancement of relic particles in our galactic core. A third strategy is to look for annihilations inside the Sun or Earth (Silk et al. 1985), where the local density of relic particles is enhanced in a calculable way by scattering off matter, which causes them to lose energy and become gravitationally bound. The signature would then be energetic neutrinos that might produce detectable muons. As also discussed here by Carr (2003), several underwater and ice experiments are underway or planned to look for this signature, and this strategy looks promising for several benchmark scenarios, as seen in the right panel of Fig. 10 <sup>1</sup><sup>1</sup>1It will be interesting to have such neutrino telescopes in different hemispheres, which will be able to scan different regions of the sky for astrophysical high-energy neutrino sources..

The most satisfactory way to look for supersymmetric relic particles is directly via their scattering on nuclei in a low-background laboratory experiment Goodman & Witten 1985). There are two types of scattering matrix elements, spin-independent - which are normally dominant for heavier nuclei, and spin-dependent - which could be interesting for lighter elements such as fluorine. The best experimental sensitivities so far are for spin-independent scattering, and one experiment has claimed a positive signal (Bernabei et al. 1998). However, this has not been confirmed by a number of other experiments, as discussed here by Kraus (2003) and Smith (2003). In the benchmark scenarios the rates are considerably below the present experimental sensitivities, but there are prospects for improving the sensitivity into the interesting range, as also discussed by Kraus (2003) and Smith (2003), as seen in Fig. 11.

Overall, the searches for astrophysical supersymmetric dark matter have discovery prospects (Ellis et al. 2001) that are comparable with those of the LHC (Battaglia et al. 2001).

### 2.5 Superheavy Dark Matter

As discussed here by Kolb (2003), it has recently been realized that interesting amounts of superheavy particles with masses $`10^{14\pm 5}`$ GeV might have been produced non-thermally in the very early Universe, either via pre- and reheating following inflation, or during bubble collisions, or by gravitational effects in the expanding Universe. If any of these superheavy particles were metastable they might be able, as seen in Fig. 12, to explain the ultra-high-energy cosmic rays that seem (Takeda et al. 2002; Abu-Zayyad et al. 2002) to appear beyond the Greisen-Zatsepin-Kuzmin (GZK) cutoff (Greisen 1966; Zatsepin & Kuzmin 1966). Such models have to face some challenges, notably from upper limits on the fractions of gamma rays at ultra-high energies, but might exhibit distinctive signatures such as a galactic anisotropy (Sarkar 2002). Examples of such superheavy particles are the cryptons found in some models derived from string theory, which naturally have masses $`10^{12\pm 2}`$ GeV and are metastable (like protons), decaying via higher-order multiparticle interactions (Ellis et al. 1990; Benakli et al. 1999). The Pierre Auger experiment (Cronin et al. 2002) will be able to tell us whether ultra-high-energy cosmic rays really exist beyond the GZK cutoff, and, if so, whether they are due to some such exotic top-down mechanism, or whether they have some bottom-up astrophysical origin. Following Auger, there are ideas for space experiments such as EUSO (Petrolini et al. 2002) that could have even better sensitivities to ultra-high-energy cosmic rays.

## 3 Calculate it!

Now that we have a good idea of the matter and energy content of the Universe, and some prospects for detecting it, the next task is to calculate it from first principles on the basis of microphysics and laboratory data.

$``$ $`\mathrm{\Omega }_b`$: As Sakharov taught us (Sakharov 1967), baryogenesis requires the violation of charge conjugation C and its combination CP with parity, interactions that violate baryon number B, and a departure from thermal equilibrium. The first two have been observed for quarks, and are expected within the Standard Model. B violation is also expected in the Standard Model, at the non-perturbative level. One might therefore wonder whether the observed cosmological baryon asymmetry could have been generated by the Standard Model alone, but the answer seems to be no (Gavela et al. 1994). However, it might be possible in the MSSM, if it contains additional sources of CP violation beyond the Standard Model (Carena et al. 2003). An attractive alternative is leptogenesis (Fukugita & Yanagida 1986), according to which first the decays of heavy singlet neutrinos create a CP-violating asymmetry $`\mathrm{\Delta }L0`$, and then this is partially converted into a baryon asymmetry by non-perturbative weak interactions.

At the one-loop level, the asymmetry in the decays of one heavy singlet neutrino $`N_i`$ due to exchanges of another one, $`N_j`$, is

$$ϵ_{ij}=\frac{1}{8\pi }\frac{1}{(Y_\nu Y_\nu ^{})_{ii}}\mathrm{Im}\left((Y_\nu Y_\nu ^{})_{ij}\right)^2f\left(\frac{M_j}{M_i}\right),$$
(31)
where $`Y_\nu `$ is a matrix of Yukawa couplings between heavy singlet and light doublet neutrinos. The expression (31) involves a sum over the light leptons, and hence is independent of the CP-violating MNS phase $`\delta `$ and the Majorana phases $`\varphi _{1,2}`$. Instead, it is controlled by extra phase parameters that are not directly accessible to low-energy experiments.

This leptogenesis scenario produces effortlessly a baryon-to-photon ratio $`Y_B`$ of the right order of magnitude. However, as seen in Fig. 13, the CP-violating decay asymmetry (31) is explicitly independent of $`\delta `$ (Ellis & Raidal 2002). On the other hand, other observables such as the charged-lepton-flavour-violating decays $`\mu e\gamma `$ and $`\tau \mu \gamma `$ may cast some indirect light on the mechanism of leptogenesis (Ellis et al. 2002$`a`$, $`b`$). Predictions for these decays may be refined if one makes extra hypotheses.

One possibility is that the inflaton might be a heavy singlet sneutrino (Murayama et al. 1993, 1994). This would require a mass $`2\times 10^{13}`$ GeV, which is well within the range favoured by seesaw models. The sneutrino inflaton model predicts (Ellis et al. 2003$`b`$) values of the spectral index of scalar perturbations, the fraction of tensor perturbations and other CMB observables that are consistent with the WMAP data (Peiris et al. 2003) . Moreover, this model predicts a branching ratio for $`\mu e\gamma `$ within a couple of orders of magnitude of the present experimental upper limit.

$``$ $`\mathrm{\Omega }_{CDM}`$: The relic density of supersymmetric dark matter is calculable in terms of supersymmetric particle masses and Standard Model parameters. The sensitivity to these parameters is quite small in generic regions, but may be larger in some exceptional regions corresponding to ‘tails’ of the MSSM parameter space (Ellis & Olive 2001). At least away from these regions, data from the LHC on supersymmetric parameters should enable the cold dark matter density to be calculated quite reliably.

$``$ $`\mathrm{\Omega }_\mathrm{\Lambda }`$: The biggest challenge may be the cosmological vacuum energy. For a long time, theorists tried to find reasons why the cosmological constant should vanish, but no convincing symmetry to guarantee this was ever found. Now cosmologists tell us that the vacuum energy actually does not vanish. Perhaps theorists’ previous failure should be reinterpreted as a success? If the vacuum energy is indeed a constant, the hope is that it could be calculated from first principles in string or M theory. Alternatively, as argued here by Steinhardt 2003, perhaps the vacuum energy is presently relaxing towards zero, as in quintessence models (Maor et al. 2002). Such models are getting to be quite strongly constrained by the cosmological data, in particular those from high-redshift supernovae and WMAP (Spergel et al. 2003) as seen in Fig. 14, and it seems that the quintessence equation of state must be quite similar to that of a true cosmological constant (Spergel et al. 2003). Either way, the vacuum energy is a fascinating discovery that provides an exciting new opportunity for theoretical physics.

## 4 Test it!

Laboratory experiments have explored the energy range up to about 100 GeV, and quantum gravity must become important around the Planck energy $`10^{19}`$ GeV: where will new physics appear in this vast energy range? There are reasons to think that the origin of particle masses will be found at some energy $`<10^3`$ GeV. If they are indeed due to some elementary scalar Higgs field, this will provide a prototype for the inflaton. If the Higgs is accompanied by supersymmetry, this may provide the cold dark matter that fills the Universe. Some circumstantial evidence for supersymmetry may be provided by the anomalous magnetic moment of the muon and by measurements of the electroweak mixing angle $`\theta _W`$, in the framework of grand unified theories. These would operate at energy scales $`10^{16}`$ GeV, far beyond the direct reach of accelerators. The first hints in favour of such theories have already been provided by experiments on astrophysical (solar and atmospheric) neutrinos, and cosmology may provide the best probes of grand unified theories, e.g., via inflation and/or super-heavy relic particles, which might be responsible for the ultra-high-energy cosmic rays. Cosmology may also be providing our first information about quantum gravity, in the form of the vacuum energy.

The LHC will extend the direct exploration of the energy frontier up to $`10^3`$ GeV. However, as these examples indicate, the unparallelled energies attained in the early Universe and in some astrophysical sources may provide the most direct tests of many ideas for physics beyond the Standard Model. The continuing dialogue between the two may tell us the origins of dark matter and dark energy.