text:
# A Cepheid is No More: Hubble’s Variable 19 in M33

## 1 Introduction

In his seminal work, “A spiral nebula as a stellar system, Messier 33,” Hubble (1926) presented the discovery of 35 Cepheid variables in that galaxy, which were used to determine its distance. No other searches for variables in M33 were undertaken for the next fifty years, until the surveys of van den Bergh et al. (1975); Sandage & Carlson (1983); Kinman et al. (1987) that discovered a total of 77 new Cepheids and recovered many of Hubble’s objects.

The DIRECT project (Kaluzny et al., 1998; Stanek et al., 1998) started in 1996 the first modern CCD-based search for variables in M33, among other targets. Our first two seasons surveyed three $`10\mathrm{}\times 10\mathrm{}`$ fields located in the central area of M33; the variable star content of these fields has been presented in Macri et al. (2001a) and Stanek et al. (2001). Since these fields had been covered by all previous surveys, we cross-correlated our variable star catalog with published catalogs and successfully identified 57 of the 60 known Cepheids. Two of the three missing variables were short period Cepheids: B08 from Sandage & Carlson (1983) (P=3.2 d) and Q19113 from Kinman et al. (1987) (P=13.4 d). These variables are faint and close to our detection threshold; therefore, it was not surprising that they eluded our detection.

The third missing variable turned out to be a surprising object, and is the topic of this paper. It is V19 from Hubble (1926), classified by him as a 54.7-day Cepheid. It is a bright, isolated, and easy to identify star located some $`6\mathrm{}`$ north of the nucleus. Hubble’s light curve exhibited a variation of $`1.1`$ mag, which is not present in our data. Additionally, the mean $`B`$ magnitude of the object increased by $`0.5`$ mag since Hubble’s original observations.

§2 summarizes the existing photometry for this object; §3 addresses the current variability, periodicity, and location of this object in the H-R and PL diagrams; and §4 discusses the possible nature of this object and its known Galactic counterparts.

## 2 Photometry

### 2.1 Data from Hubble (1926)

The discovery of V19 by Hubble (1926) was based on observations obtained at the Mount Wilson 100-inch telescope in 56 epochs between August 1919 and September 1925, with two additional plates from September 1909 and August 1910. Magnitudes were measured using the Argelander method, and were based on local reference stars calibrated through a photometric sequence from the Selected Area 45. Hubble’s original magnitudes suffered from zeropoint and scale errors, as pointed out by Sandage (1983) and Christian & Schommer (1987).

We calculated a new transformation between Hubble’s magnitude scale and the standard system, using our photometry (Macri et al., 2001b) for Hubble’s original comparison stars. We identified 121 of these stars in our CCD images, using the finding charts published by Sandage (1983) and performed a least-squares fit between Hubble’s photographic magnitudes and our CCD-based $`B`$ magnitudes. We found

$$m_B=1.62(\pm 0.04)\left[m_{pg}18.5\right]+19.68(\pm 0.02),$$
(1)
where $`m_B`$ is the $`B`$ magnitude in the standard system and $`m_{pg}`$ is the photographic magnitude as measured by Hubble. The relation has a scatter of 0.18 mag after the rejection of 11 stars which deviated by more than 0.5 mag. Our coefficients agree with those derived by Christian & Schommer (1987). We corrected Hubble’s original measurements using the newly derived transformation, and analyzed the data using a suite of programs developed by the DIRECT project (Kaluzny et al., 1998). We derived a best-fit period of 54.71 days, identical to Hubble’s original value, and an intensity-weighted mean $`B`$ magnitude of $`19.59\pm 0.23`$ mag. Figure 1 shows a phased light curve of V19 based on Hubble’s recalibrated data. The scatter in the light curve is consistent with the r.m.s. scatter of Equation 1.

### 2.2 DIRECT Data

We based our identification of V19 on the finding charts of Hubble (1926), van den Bergh et al. (1975) and Humphreys & Sandage (1980), all of whom mark the same star as the variable. V19 is the DIRECT object D33J013357.1+304512.4, located at R.A.=1h33m57.1s, Dec.=+30d45m12.4s (J2000.0). We ruled out a mis-identification of V19 with a nearby star in Hubble’s original finding chart by identifying all variables in the general vicinity ($`<3\mathrm{}`$). This is a generous search area, given the fact that all other variables discovered by Hubble were recovered within a few arcseconds of previously published positions. The nearest variable of any kind is D33J013357.1+304455.2, a faint ($`V=21.5`$) periodic variable located $`17\mathrm{}`$ south of V19. The nearest Cepheid is D33J013354.9+304532.5, with a period of 36 days and located $`39\mathrm{}`$ north-west of V19. The nearest Cepheid with a period close to Hubble’s original period for V19 is D33J013347.5+304423.2, with a period of 56 days and located $`2.5\mathrm{}`$ to the south-west.

V19 was observed by the DIRECT project a total of 11, 116 and 60 times in $`BVI`$, respectively, during the 1996 and 1997 observing seasons. The description of the reduction and analysis of this data set was presented in Macri et al. (2001b). Additional observations of another DIRECT field which contains V19 were obtained during the 1999 season; for the purposes of this paper, we performed a preliminary reduction using the standard DIRECT pipeline (Kaluzny et al., 1998) and calibrated the resulting photometry using bright, non-variable stars from the database of Macri et al. (2001b). This third year of observations yielded an additional 19, 121 and 35 points in $`BVI`$, respectively. The mean $`BVI`$ magnitudes of V19 obtained from the three seasons spanned by our data are $`19.08\pm 0.05`$, $`18.21\pm 0.05`$ and $`17.30\pm 0.05`$ mag, respectively. The light curves exhibit very little variation, of order 0.05 mag peak-to-peak. §3 contains the analysis of our light curve data333The light curve of V19 and its finding chart can be found at the DIRECT Web page: http://cfa-www.harvard.edu/~kstanek/DIRECT.

### 2.3 Other photometric data

We searched the literature for additional observations of V19 obtained after 1926 and before 1996. The only time series photometry that we were able to locate was that of van den Bergh et al. (1975), who obtained 67 plates of M33 between January 1966 and November 1974 at the Palomar 48-inch Schmidt, using a variety of plate types and filters. The data for V19 consists of 19 points in $`B`$ and 41 points in $`V`$. Unfortunately, the uncertainties quoted by the authors for the individual measurements is 0.7 mag, and there are no data available to check the absolute calibration of their photometric systems.

V19 was observed by Humphreys & Sandage (1980) in their survey for blue and red stars in M33, and was designated as R142 in their catalog. We corrected their tabulated photometry according to Wilson et al. (1990) and obtained $`B=19.1`$ mag and $`BV=0.9`$ mag, consistent with our own determinations. Unfortunately, the authors did not obtain a spectrum for this star. Additionally, V19 is present in the fields surveyed by Kinman et al. (1987), but these authors do not list it in their compilation of variables.

Madore et al. (1985) reported the results of $`JHK`$ observations of V19 and other M33 Cepheids conducted between 1980 and 1983 at a variety of telescopes using a photometer with a $`5\mathrm{}`$ diameter aperture. However, the authors expressed concerns about possible contamination by field stars inside their aperture. Therefore, we analyzed $`K_s`$-band images of V19, obtained by Don McCarthy with the PISCES camera (McCarthy et al., 1998) at the refurbished MMT, under $`0.6\mathrm{}`$ seeing. These images show the presence of another star of similar infrared flux some $`3\mathrm{}`$ away. A preliminary calibration of the MMT data, based on isolated 2MASS stars present in the field of view, yields a mean magnitude for V19 of $`K_s=16.1\pm 0.05`$ mag. This, in turn, implies a $`VK`$ color of $`2.1\pm 0.1`$ mag, which is somewhat redder than, but consistent with, the one predicted for a 55-day Cepheid, $`1.9\pm 0.1`$ mag (c.f. Equations 5 and 9 of Macri et al., 2001c).

Lastly, we retrieved archival HST/WFPC2 frames in which V19 is present (dataset u2c604, proposal # 0538) from the Space Telescope Archive. The images were obtained in the $`F439W`$ and $`F547M`$ filters and show the object (at x=92,y=283 on WF4) to be well isolated, with no detectable companions inside our ground-based seeing disk.

## 3 Analysis

We computed the $`J_S`$ index (Stetson, 1996) using all 237 $`V`$ measurements from the DIRECT data, and found a value of $`J_S=0.807`$, which is above the standard DIRECT threshold of 0.75. However, the overall r.m.s. scatter is only 0.025 mag, below our usual threshold of 0.04 mag. Therefore, the star does not meet the standard criteria for variability adopted by the DIRECT project.

We attempted to find a period from our data, using the standard DIRECT technique of a modified Lafler-Kinman algorithm (Stetson, 1996) and the CLEAN algorithm of Roberts et al. (1987). No convincing periods were found, which is not surprising given that our frame-to-frame random photometric error (arising from the limitations of fixed-position PSF photometry) is of the order of $`0.020.03`$ mag. We conclude that if the 54.7-day periodicity seen by Hubble is still present, its amplitude must be less than 0.05 mag. Long-term follow-up observations using larger telescopes would be most useful to detect any small-scale periodicity and/or variability.

Figure 2 shows the $`B`$-band light curve of V19 from 1910 to 2000, based on the data from Hubble (1926), vhk75, and this work. The intensity-weighted mean magnitude shows a dramatic increase, from $`19.59\pm 0.23`$ mag in 1919-1925 to $`B=19.08\pm 0.05`$ in 1996-1999, while the amplitude of pulsation decreases from 1.1 mag to undetectable levels. Figure 3 shows the location of V19 in the H-R diagrams of Macri et al. (2001b). We also plot the location of 130 Cepheids with $`P>10`$ days present in the inner region of M33 (Macri et al., 2001a). Figure 4 shows the location of V19 in the $`BVI`$ Cepheid P-L relations defined by the same set of variables, using Hubble’s original period of 54.7 days. These figures indicate that the modern optical magnitudes of V19 are somewhat brighter than, but still consistent with, those of 55-day Cepheids in M33.

## 4 Discussion

The cessation of pulsations, or any significant decrease in amplitude (i.e., by a factor of ten), has never been observed in a classical Cepheid. V19 may be the first case, or it may belong to another class of variable stars. We present here four classes of variable stars whose pulsational behaviour may explain the nature of V19.

Population II Cepheids: RU Cam, a 22-day Cepheid which in 1965 abruptly decreased in amplitude from 1 mag to about 0.1-0.2 mag (Demers & Fernie, 1966), was later shown to be a W Virginis variable (Population II Cepheid) that exhibited a highly unstable and modulated light curve thereafter (Kollath & Szeidl, 1993). Even before the 1965 event, Payne-Gaposchkin (1941) had suggested that RU Cam may be a W Virginis variable. However, RU Cam is of little relevance to our case, because in order for V19 to be a Population II Cepheid, it should be about 2 magnitudes fainter (see Alcock et al., 1998) than Population I Cepheids of similar period. Furthermore, W Virginis variables are not known to change their mean magnitude, although a few RV Tauri variables (longer-period Population II Cepheids) do that cyclically.

Peculiar Population I Cepheids: These objects are extremely rare and entirely confined to much shorter periods. A case in point is the classical Cepheid V473 Lyr of $`P=1.49`$ days (e.g., Burki et al., 1986; Andrievsky et al., 1998). It has an amplitude modulation of a factor of 15 which resembles the beating of two closely spaced pulsation modes. This phenomenon is not understood, but it is most likely associated with the interaction of high-order modes. Furthermore, the interval corresponding to low amplitudes in V473 Lyr has a short duration; in general, the phenomenon seems to have very little in common with V19. A more relevant object in this category could be Polaris, now established to be a Population I Cepheid pulsating in the first overtone mode (Feast & Catchpole, 1997). It has a period is 3.97 days, and its amplitude has decreased by a factor of 3 over the past 50 years. Despite earlier reports, this amplitude decrease has stabilized over the past 15 years; the pulsation amplitude is practically constant since 1986 (Hatzes & Cochran, 2000). As in the case of V473 Lyr, the amplitude decrease could be due to high-order mode interactions (Evans et al., 2001). However, the pulsation amplitude of Polaris has always been very low ($`0.1`$ mag in V), which makes the comparison with V19 difficult. Finally, we should note that there are stars inside the Cepheid instability strip that are not variable, for reasons that are still unknown.

UU Herculis stars: These are supergiants with Cepheid-like pulsations that exhibit occasional standstills (Sasselov, 1984). Pulsations cease abruptly and at mid-cycle (corresponding to the mean magnitude of the star), last a couple of months, and then abruptly start again (Sasselov, 1983). The UU Herculis stars may alternate between Cepheid-like and RV Tauri-like pulsations (over timescales of several years) with periods of 30 to 90 days, which are always of relatively low amplitude ($``$0.3 mag in V). UU Her is the only exception, with amplitudes as high as 0.6 mag in V. Most stars in this class have infrared excesses, low metallicities, and show variable emission in the H$`\alpha `$ line (Klochkova et al.,, 1997).

Luminous Blue Variables: Given its luminosity (M$`{}_{\mathrm{V}}{}^{}6.4`$ mag), V19 appears to be more closely related to the Luminous Blue Variables (LBVs), which are not always blue. One example is 164 G Sco (M$`{}_{\mathrm{V}}{}^{}=8.4`$, T=10,000K, $`R175R_{}`$, after Lamers et al., 1998). It shows multiperiodic variability ($`4555`$ days) of very low amplitude ($`0.1`$ mag in V). Occasionally, LBVs experience dramatic outbursts and mass loss episodes which could considerably alter their effective radii and temperatures for extended periods of time. There also exist other luminous supergiants which are not LBVs but exhibit Cepheid-like pulsations, such as V810 Cen (Kienzle et al., 1998). However, their pulsations have always been observed to be low ($`0.1`$ mag) and multiperiodic.

Each of the types described above belongs to a different evolutionary state. The W Virginis (and RV Tauri) variables are most likely low-mass AGB (and tip-of- or post-AGB) stars (hydrogen- and helium-shell burning). The same appears to be true for the majority of UU Hers stars (a related object which ”wanders” across the HR diagram is FG Sge). Hence their lower luminosity and common circumstellar infrared excesses. They often show emission in the absorption profile of the H$`\alpha `$ and other strong lines. There is no evidence that V19 possesses any of the above characteristics and its luminosity is incompatibly higher. The remaining types of variables (LBVs, V810 Cen) are very massive supergiants. At a given pulsation period they are brighter than Population I Cepheids and they are unlikely to be mistaken for Cepheids, due to their low amplitudes and multi-periodic pulsation.

In summary, V19 does not belong to any of the above types: it is too luminous to be a low-mass Population II star, and its pulsation was too regular and high-amplitude 75 years ago to be a massive microvariable supergiant. Therefore, V19 appears to be a high-mass Population I star that behaves more like a low-mass Population II RV Tauri variable.

Given the bright and isolated nature of V19, and its location in the host galaxy, it is likely that most photographic plates of M33 obtained at 2-m class (or larger) telescopes during the 20th century will yield precise photometry for this object. We encourage the astronomical community to contact us if they are in possession of plate material that might help constrain the past behavior of V19. Future observations, including spectroscopy and long-term, multi-band photometric follow-up, will help resolve the nature of this object. Such resolution is desirable in order to avoid misclassifying objects like V19 as classical Cepheids in distant galaxies.

## 5 Acknowledgments

We would like to thank Don McCarthy for the PISCES images of V19 taken at the new MMT. Thanks are due to D. Bersier, N.R. Evans, J. Kaluzny & B. Paczynski for helpful comments. LMM would also like to thank John Huchra for his support and comments. DDS acknowledges support from the Alfred P. Sloan Foundation and from NSF grant No. AST-9970812. KZS was supported by Hubble Fellowship grant HF-01124.01 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.


format: text

text:
# 1 Introduction

## 1 Introduction

Solitons are special solutions of non-linear wave equations. The most relevant characteristic of solitons is that they are localized static solutions. The simplest example consist of a single scalar field $`\varphi `$ in one spatial and one temporal dimensions. Perhaps the most famous one is the sine-Gordon soliton ,.

At first sight, it can be thought that a wave equation with a single scalar field in more than three spatial dimensions with solitonic solutions can be found. However there is a Virial theorem which restricts this possibility. Here we are going to transcript that theorem and it’s proof for the convenience of the reader :

###### Theorem 1

There are no non-trivial static solitary waves of systems with scalar fields when the space dimensionality is three or more and when the Lagrangian has the form:

$$(𝐱,\mathrm{t})=\frac{1}{2}(_\mu \varphi )(^\mu \varphi )\mathrm{U}(\varphi (𝐱,\mathrm{t}))$$
(1)
with $`\varphi =[\varphi _i(𝐱,\mathrm{t});\mathrm{i}=1,\mathrm{},\mathrm{N}]`$ a set of N coupled scalar fields in D space plus one time dimensions, and $`U(\varphi (𝐱,\mathrm{t}))`$ a positive definite potential.

###### Proof 1

A static solution $`\varphi (𝐱)`$ obeys

$$^2\varphi =\frac{U}{\varphi }(𝐱)$$
(2)
where $`^2`$ is the Laplacian in D dimensions. This equation clearly correspond to the extremum condition $`\delta W=0`$ for the static energy functional

$$W[\varphi ]d^Dx\left[\frac{1}{2}_i\varphi _i\varphi +U(\varphi (𝐱))\right]V_1[\varphi ]+V_2[\varphi ]$$
(3)
where the functionals $`V_1`$ and $`V_2`$ stand for the two terms on the right-hand side. Note that not only W but also $`V_1`$ and $`V_2`$ are non-negative. Now, let $`\varphi _1(𝐱)`$ be a static solution. Consider the one-parameter family of configurations

$$\varphi _\lambda =\varphi _1(\lambda 𝐱).$$
(4)
It is easy to check that

$$W[\varphi _\lambda ]=V_1[\varphi _\lambda ]+V_2[\varphi _\lambda ]=\lambda ^{2D}V_1[\varphi _1]+\lambda ^DV_2[\varphi _1].$$
(5)
Since $`\varphi _1(𝐱)`$ is an extremum of $`W[\varphi ]`$, it must in particular make $`W[\varphi _\lambda ]`$ satationary with respect to variations in $`\lambda `$; that is,

$$\frac{d}{d\lambda }W[\varphi _\lambda ]=0at\lambda =1.$$
(6)
Differentiating (5) using (6) gives us

$$(2D)V_1[\varphi _1]=DV_2[\varphi _1].$$
(7)
Since $`V_1`$ and $`V_2`$ are non-negative (7) cannot be satisfied for $`D3`$ unless $`V_1[\varphi _1]=V_2[\varphi _1]=0`$. This means that $`\varphi _1(𝐱)`$ has to be space-independent and equal to one of the zeros of $`U[\varphi ]`$. This is just a trivial solution and the theorem precludes non-trivial space-dependent solutions. q.e.d.

For the case of $`D=2`$, equation (7) tell us that $`V_2[\varphi _1]=0`$. The simplest example of a solution of this kind is the non-linear O(3) model ,, relevant in the description of the statical mechanics of an isotropic ferromagnet.

We are interested in the possible existence of a static solution of a system in a curved space-time in three spatial and one temporal dimensions. We will construct a simple model and look for solitonic solutions numerically.

## 2 The model

Consider the simplest case in a curved space-time. One scalar field whose equation of motion is:

$$^\mu _\mu \varphi \frac{V}{\varphi }=0$$
(8)
Without potencial (8) is the Laplace equation,whose solutions take a maximum or minimum value only at the spatial boundaries. If we solve the Laplace equation in a space with a connected boundary or in a compact space, the solution in every point of the space will necessarily have the same value as in the boundary (i.e. a trivial solution).

So we are going to introduce a potential $`V=(\varphi ^21)^2`$ in analogy with the $`1+1`$ dimensional case. Moreover, we consider the simplest kind of spatialy compact universe, the static Einstein Universe:

$$ds^2=dt^2+d\chi ^2+sin^2\chi (d\theta ^2+sin^2\theta d\phi ^2).$$
(9)
Equation (8) in the space-time corresponding to (9) takes the form:

$$\frac{^2\varphi }{\chi ^2}\frac{^2\varphi }{t^2}+2cot\chi \frac{\varphi }{\chi }+sin^2\chi cot\theta \frac{\varphi }{\theta }+sin^2\chi \frac{^2\varphi }{\theta ^2}+sin^2\chi sin^2\theta \frac{^2\varphi }{\phi ^2}=2(\varphi ^21)2\varphi .$$
(10)
We are interested in static solutions and with spherical simetry so we have

$$\frac{d^2\varphi }{d\chi ^2}+2cot\chi \frac{d\varphi }{d\chi }4\varphi (\varphi ^21)=0$$
(11)
this is an ordinary differential equation of second order. Equation (11) can be separated in 2 ordinary equations of first order. Let:

$$x_1=\varphi x_2=\frac{d\varphi }{d\chi }$$
(12)
so, (11) can be written as the following system:

$$x_2=\frac{dx_1}{d\chi }$$
(13)
$$\frac{dx_2}{d\chi }+2x_2cot\chi 4x_1(x_1^21)=0.$$
(14)
These equations can be integrated by the Runge-Kutta method .

We note that equation (11) is singular at $`\chi =0`$ and $`\chi =\pi `$, so we use the following initial conditions:

$$\varphi (0)=\varphi _a^0\frac{d\varphi (0)}{d\chi }=\dot{\varphi }_a^0=0$$
(15)
$$\varphi (\pi )=\varphi _b^0\frac{d\varphi (\pi )}{d\chi }=\dot{\varphi }_b^0=0$$
(16)
So we have a differential equation (11) with boundary conditions at the two end points. It can be solved using a “shooting to a middle point” method, with $`\varphi _a^0`$ and $`\varphi _b^0`$ as shooting parameters. The main idea of this method is the following:

We perform an integration from $`\chi =0`$ to $`\chi =\frac{\pi }{2}`$ obtaining $`\varphi _a(\frac{\pi }{2})`$. Also we integrate from $`\chi =\pi `$ to $`\chi =\frac{\pi }{2}`$ obtaining $`\varphi _b(\frac{\pi }{2})`$. Now, let construct the following function:

$$F(\varphi _a,\varphi _b)=[\varphi _a\left(\frac{\pi }{2}\right)\varphi _b\left(\frac{\pi }{2}\right),\dot{\varphi }_a\left(\frac{\pi }{2}\right)\dot{\varphi }_b\left(\frac{\pi }{2}\right)]$$
(17)
The first integration allows us to evaluate F obtaining in general $`F(\varphi _a^0,\varphi _b^0)0`$ (figure 1). We are interested in the zeros of this function, because at those points the solutions of the integration match $`\varphi `$ and $`\dot{\varphi }`$ in a smooth way. We used the Newton-Raphson method to find the zeros of a function , and the result is showed in figure 2.

## 3 Conclusion

We have shown that the theorem that precludes the existence of solitonic solutions to systems based in saclar field in $`3+1`$ or more flat space-time dimensions would be false if extended to the curve space-time case. We have done so by explicitly constructing one such solitonic solution of a simple model based on a single scalar field in the static Einstein Universe.

This is another indication333The most widely knows example of this phenomenon arose in the consideration of Einstein-Yang-Mills theory where solitonic solutions have been found while it know that these are no such solutions in Yang-Mills theory in Minkowskian space-time . that the interplay of curved space-time physics and soliton physics allow us for richer phenomena than each either of the two fields on its own.

## 4 Aknowledgments

This work was partially supported by DGAPA-UNAM grant #IN121298 and CONACyT grant #32272E.


text:
# Conductivity of the classical two-dimensional electron gas

## Abstract

We discuss the applicability of the Boltzmann equation to the classical two-dimensional electron gas. We show that in the presence of both the electron-impurity and electron-electron scattering the Boltzmann equation can be inapplicable and the correct result for conductivity can be different from the one obtained from the kinetic equation by a logarithmically large factor.

The description of kinetics of the classical electron gas in semiconductors is usually done with the help of the Boltzmann kinetic equation

$`{\displaystyle \frac{n_𝒑(𝒓,t)}{t}}+𝒗{\displaystyle \frac{n_𝒑(𝒓,t)}{𝒓}}+e𝑬{\displaystyle \frac{n_𝒑(𝒓,t)}{𝒑}}`$ (1)
$`=I_{imp}(n_𝒑(𝒓,t))+I_{ee}(n_𝒑(𝒓,t))+I_{ph}(n_𝒑(𝒓,t))`$ (2)
where $`n_𝒑(𝒓,t)`$ is the electron distribution function, $`𝑬`$ is the external electric field, $`𝒗`$ and $`𝒑`$ are the electron velocity and momentum, and $`I_{imp},I_{ee}`$ and $`I_{ph}`$ are the electron-impurity, electron-electron and electron-phonon scattering integrals which are characterized by the relaxation times $`\tau _i`$, $`\tau _e`$ and $`\tau _{ph}`$ respectively.

At low enough temperatures $`T`$ when $`\tau _i\tau _{ph},\tau _e`$, the conductivity of the sample is determined by the electron-impurity scatterings and is given by the Drude formula

$$\sigma _0=e^2D\nu $$
(3)
Here $`D=v_F^2\tau _i/2`$ (where $`v_F`$ is the Fermi velocity, $`\tau _i=1/v_FN_ia`$, $`N_i`$ is the impurity concentration and $`a`$ is the amplitude of the electron-impurity scattering), and $`\nu `$ is the electron density of states at the Fermi level. To be concrete, we consider the case of a degenerate electron gas, though our results are of a general character. Eq. (3) also holds in the case of intermediate temperatures when $`\tau _{ph}\tau _i\tau _e`$. Indeed, since the electron-electron scatterings conserve the total momentum in the electron system, they do not contribute directly to the resistance of the system.

The form of $`I_{imp}`$ in Eq. (2) is a result of a procedure in which averaging over realizations of the random impurity potential is done before solving the equation of motion. The criteria of applicability of the conventional Boltzmann kinetic equation (2) correspond to the absence of correlations between subsequent scattering events. This leads us to a requirement that the electron gas is weakly interacting

$$\tau _i,\tau _et_{coll}$$
(4)
where $`t_{coll}`$ is the the collision time, which is estimated as $`a/v`$ in classical mechanics, and as $`\frac{\mathrm{}}{ϵ}`$ in quantum mechanics ($`ϵ`$ is the characteristic electron energy).

It has been shown  that in the two-dimensional case even under conditions (4) the quantum interference corrections to the Drude formula logarithmically diverge as temperature $`T`$ goes to zero. This means that in this case the kinetic equation (2) is not applicable at small $`T`$.

In this article we study the case of high enough temperatures when the quantum interference corrections can be neglected. The classical corrections to Eq. (3) associated with correlations between subsequent scatterings were considered in  in the framework of the Lorentz gas model. It was shown that in the two-dimensional case correlational corrections to the diffusion coefficient and to the conductivity are small, but nonanalytical functions of the concentration of scatterers $`N_i`$.

We consider in this article the case of a two-dimensional electron gas with the mean free path $`l_e=v_F\tau _e`$ much smaller than the average distance between the impurities: $`l_eN_i^{1/2}`$. We show that under conditions (4) a combination of electron-electron and electron-impurity scatterings can lead to big correlational corrections to the Drude formula and that the correct expression for the conductivity of the system is

$$\sigma =\sigma _0\left(1+\frac{av_F}{2\pi \eta }\mathrm{ln}\frac{1}{N_i^{1/2}l_e}\right)$$
(5)
where $`\eta `$ is the kinematic viscosity of the electron liquid, which is estimated for a Fermi liquid as $`\eta v_Fl_e`$. Eq. (5) differs from Eq. (3) by a logarithmically large factor, provided

$`{\displaystyle \frac{av_F}{2\pi \eta }}\mathrm{ln}{\displaystyle \frac{1}{N_i^{1/2}l_e}}1.`$

On distances from impurities larger than $`ł_e`$, the motion of the electron system can be described by the Navier-Stokes hydrodynamic equation which can be derived from Eq. (2) with $`I_{imp}=0`$ . In this case the elastic scatterers should be described by a random potential and the averaging over the random realizations of the potential should be done after solving the Navier-Stokes equation. This should in principle give a result for conductivity that is different from the one implied by Eq. (2). The two approaches are equivalent in the approximation which assumes the variation of the hydrodynamic velocity between impurities to be small.

The idea of this paper is that, in fact, in the two-dimensional case the hydrodynamic velocity of electrons is a significantly nonuniform function of coordinates: the velocity near impurities is much smaller than inbetween them. This is in contrast to the three-dimensional case where the spatial variation of the velocity is small and Eqs. (2), (3) are applicable. The reason for the significant nonuniformity of the velocity field is connected to the Stokes paradox : in the presence of an obstacle bypassed by the fluid, the solution of the linearized Navier-Stokes equation in the two-dimensional case grows logarithmically with the distance from the obstacle. This means the absence of a bounded solution. In our case the finite concentration of scatterers provides a bound on the magnitude of the velocity field. However, the magnitude of the average velocity of electrons induced by the electric field turns out to be logarithmically bigger than that which follows from Eq. (2).

In the approximation linear in $`𝑬`$, the Navier-Stokes equation has the form

$$p+\mu \mathrm{\Delta }𝒖+ne𝑬\underset{i}{}𝒇\delta (𝒓𝒓_i)=0$$
(6)
Here $`n`$ is the electron concentration, $`\mu =mn\eta `$ is the dynamic viscosity of the electron fluid (where $`m`$ is the electron mass), and $`𝒓_i`$ are coordinates of impurities labeled by the index ”i”. The last term in Eq. (6) is associated with the force on the fluid by the i-th impurity, whose radius is much smaller than $`l_e`$. Therefore the value of the force $`𝒇`$ can be estimated assuming that to zeroth order in $`a/l_e`$ the distribution function of the electrons near the $`i`$-th impurity ($`|rr_i|l_e`$) is an equilibrium one with the average velocity $`𝒖(𝒓_i)`$. (In the following we will assume that all $`𝒖(𝒓_i)𝒖_0`$ are of the same order.) Consequently we have

$$𝒇=np_Fa𝒖_0.$$
(7)
In the stationary case the conservation of the total momentum of the system gives

$$e𝑬=𝒇N_i/n.$$
(8)
Far from the position of impurities ($`|𝒓𝒓_i|l_e`$) the electron liquid can be considered as an incompressible one. The continuity equation for an incompressible liquid

$`𝒖=0`$

and the Navier-Stokes equation (6) have a solution for the pressure and the velocity

$`𝒖(𝒓)`$ $`=`$ $`C\widehat{x}+\widehat{x}{\displaystyle \frac{C_0}{2}}\mathrm{ln}r{\displaystyle \frac{C_0}{4}}(\widehat{r}\mathrm{cos}\theta +\widehat{\theta }\mathrm{sin}\theta )`$ (9)
$`p(r)`$ $`=`$ $`p_0\mu C_0{\displaystyle \frac{\mathrm{cos}\theta }{r}}`$ (10)
where $`C_0=\frac{f}{2\pi \mu }`$, $`\widehat{r}`$ and $`\widehat{\theta }`$ are unit basis vectors of the cylindrical coordinate system and the x-axis (with the unit vector $`\widehat{x}`$) is oriented parallel to the electric field.

The value of the coefficient $`C`$ is determined using $`𝒖(r=l_e,\theta =\pi )=𝒖_0`$ and Eqs. (7), (8) and (10):

$`C={\displaystyle \frac{\sigma _0E}{ne}}(1{\displaystyle \frac{av_F}{4\pi \eta }}\mathrm{ln}l_e).`$

As a result, the spatially averaged velocity $`\overline{𝒖}`$ of the electronic flow is

$`\overline{𝒖}={\displaystyle \frac{\sigma _0𝑬}{ne}}(1+{\displaystyle \frac{av_F}{2\pi \eta }}\mathrm{ln}{\displaystyle \frac{1}{N_i^{1/2}l_e}}).`$

Expressing the current density as $`𝒋=ne\overline{𝒖}`$, the electric conductivity Eq. (5) is obtained.

The calculations presented above assumed that the viscosity of the electron gas is a well-defined quantity. On the other hand it is well known that the fluctuational corrections to the viscosity of a 2D neutral liquid diverge , which means that hydrodynamic equations are nonlocal in two dimensions. This fact was pointed out in  and discussed in the most complete form in . Below we present a trivial generalization of the procedure  to estimate fluctuational corrections to the viscosity of the two-dimensional charged electron liquid. In this case the low-lying collective excitations are plasmons and the fluctuations of the charge and of the hydrodynamic velocity can be described as a combination of plasmon excitations with different wave vectors $`𝒒`$. The real part of the plasmon frequency is  $`\omega _𝒒=\sqrt{2\pi ne^2|𝒒|/m}`$ , where $`𝒒`$ is the plasmon momentum. Since the imaginary part of the frequency $`\omega _𝒒^{}=\eta |𝒒|^2/2`$ is proportional to $`|𝒒|^2`$ ($`\eta =\mu /mn`$ is the kinematic viscosity of the liquid), fluctuations with sufficiently low $`|𝒒|`$ have arbitrarily large mean free path and therefore plasmons are well-defined excitations. This is the physical reason why the fluctuational mechanism of viscosity is always the basic one at sufficiently small gradients and in samples with no elastically scattering potential. Since the phases of the thermal fluctuations are random it is convenient to introduce a distribution function of plasmon fluctuations $`n(𝒒)`$, normalized in such a way that the quantity $`\frac{d𝒒}{(2\pi )^2}\omega (𝒒)n(𝒒)`$ is the average energy density of fluctuations. The flux of momentum associated with these excitations is

$$\mathrm{\Pi }_{ij}=\frac{d𝒒}{(2\pi )^2}q_j\frac{\omega _𝒒}{q_i}n(𝒒)$$
(11)
The part of the distribution function of plasmons which is proportional to the gradient of the hydrodynamic velocity can be retrieved from the kinetic equation

$$\frac{\omega _𝒒}{𝒒}\frac{n_𝒖(𝒒)}{𝒓}=(n(𝒓,𝒒)n_u(𝒒))/\tau $$
(12)
where $`\tau =\omega ^1`$ is the plasmon mean free time and $`n_𝒖(𝒒)`$ is the equilibrium Bose distribution function corresponding to a local hydrodynamic velocity $`𝒖(𝒓)`$. At small $`\omega `$ and $`|𝒖|`$ one can approximate $`n_𝒖(𝒒)\frac{T}{\omega _𝒒𝒒𝒖(r)}`$. Using Eq. (12) and expanding $`n_𝒖(𝒒)`$ in $`𝒖`$, one can estimate from Eq. (11) the fluctuational contribution to the viscosity

$$\delta \eta \frac{1}{mn}(q\frac{\omega _𝒒}{q})^2\tau \frac{n_0(𝒒)}{\omega _𝒒}q𝑑q$$
(13)
(where $`n_0(𝒒)=n_{𝒖=0}(𝒒))`$. Thus the fluctuational correction to the viscosity can be described only in terms of hydrodynamic quantities: viscosity and plasmon dispersion relation. In two-dimensional pure systems the expression Eq. (13) logarithmically diverges. However, in the presence of the elastically scattering potential the dispersion relation of plasmons which we used holds only at those frequencies for which the plasmon mean free time is shorter than the electron-impurity scattering time $`\tau <\tau _i`$. This determines a cut-off in Eq. (13). As a result, the fluctuational correction to the viscosity can be neglected provided

$`\delta \eta {\displaystyle \frac{T}{\mu }}\mathrm{ln}({\displaystyle \frac{T}{E_F}}{\displaystyle \frac{\sigma _0}{e^2/\mathrm{}}})\eta .`$

In conclusion, we would like to stress that the existence of the large logarithmic factor in the Eq. (5) and the consequent inapplicability of the Boltzmann kinetic equation is a property of the two-dimensional problem. In the three-dimensional case the Stokes paradox does not exist, the difference between $`𝒖_0`$ and $`\overline{𝒖}`$ is small and the Boltzmann kinetic equation can be applied.

This work was partially supported by the Division of Material Sciences, U.S.National Science Foundation under Contract No. DMR-9205144.

We would like to thank A. Zyuzin, A.Andreev, S.Chakrovarty, M.Feigelman, A.I.Larkin, D.Khmelnitskii, S.Kivelson, F.Zhou and P.Kovtun for valuable discussions.


text:
# 1 Introduction

## 1 Introduction

One of the most exciting unresolved issues in particle physics and cosmology is the nature of dark matter (DM) in the Universe. Cosmological and astronomical observations have put, over the years, significant restrictions on its expected properties. In particular, it is widely expected that a hypothetical candidate should be in the form of cold DM (CDM), while a combination of several measurements has narrowed the range of its relic abundance to

$$0.1\text{ }<\mathrm{\Omega }_{\mathrm{CDM}}h^2\text{ }<0.3,$$
(1)
where $`\mathrm{\Omega }_{\mathrm{CDM}}`$ is the ratio of the CDM relic density to the critical density and $`h0.7`$ is a parameter in the Hubble constant $`H_0=100h`$ km/sec/Mpc . Constraints from the Big Bang Nucleosynthesis  and, recently, from observations of the cosmic microwave background radiation  have strongly restricted the abundance of baryonic DM; for example, Ref.  quotes the range $`\mathrm{\Omega }_bh^2`$ $`=`$ $`0.02\pm 0.002`$. Therefore most of the DM in the Universe is expected to be non-baryonic.

From a particle physics point of view, CDM is most naturally made of some weakly interacting (stable) massive particles (WIMPs). Supersymmetry (SUSY) , which over the years has emerged as a leading candidate for new physics beyond the Standard Model (SM), predicts that, in the presence of R-parity, the lightest supersymmetric particle (LSP) is stable. This alone makes the LSP an interesting candidate for a WIMP and DM in the Universe.

A SUSY candidate for the LSP WIMP is by no means unique. In models coupled to gravity there exists the gravitino, the SUSY partner of the graviton. The gravitino has long been known to be a potential candidate for DM, although it generically suffers from a well-known “gravitino problem” . In SUSY models that incorporate the Peccei–Quinn solution to the strong CP problem, there exists the axino, the fermionic partner of the axion. The axino has recently been shown to be an attractive and well-motivated alternative candidate for the LSP and CDM . A number of more speculative possibilities, stemming from, for example, string theories, have also been proposed. However, in some sense the simplest choice is to consider one of the superpartners of the SM particles as a potential LSP. Among these, the lightest neutralino $`\chi `$ stands out as probably the most natural and attractive candidate for the LSP and DM . In this paper we will focus on this case.

In order to be able to perform a reliable comparison between theoretical predictions and improving measurements of the relic abundance and limits from underground DM searches, a precise calculation of the relic density becomes necessary. The constraint (1) is known to often provide a strong restriction on allowed combinations of SUSY masses and couplings. For example, in the minimal supergravity model (also dubbed CMSSM) the condition $`\mathrm{\Omega }_\chi h^2<0.3`$ provides a strong upper bound of a few hundred GeV on the masses of gauginos and scalars over an overwhelming fraction of the parameter space . When combined with recent LEP and other data, the remaining CMSSM parameter space becomes highly restricted . In the future, the range (1) is expected to be significantly narrowed by MAP and Planck, which in turn will lead to even stronger constraints on allowed combinations of SUSY masses and couplings.

One of the key steps in computing the relic density of the neutralino involves calculating the thermal average of its annihilation cross section times the relative velocity. This calculation is technically rather involved and, therefore, in most of the literature, approximate methods have been used. These methods are based on expanding the thermal average in powers of $`xT/m_\chi `$, where $`T`$ is the temperature of the thermal bath and $`m_\chi `$ is the neutralino mass. One usually computes only the first two terms of the expansion. It has been well known that the expansion fails badly near $`s`$-channel resonances  and near new final-state thresholds . In particular, it was emphasized in Ref.  that, owing to the very narrow width of the lightest SUSY Higgs boson $`h`$, the error can be as large as a few orders of magnitude very close to the $`h`$-resonance. However, even though the usual expansion is generally expected to be accurate enough further away from resonances and thresholds, no detailed study of this point has yet been done.

On the other hand, while a general formalism for computing the thermal average precisely was derived a long time ago , a full set of exact, general and explicit analytic expressions for the neutralino pair-annihilation cross section is still not available in the literature. (The thermal average can be obtained by performing a single integral over the cross section, as we will see below.) The cross sections for the neutralino pair-annihilation into the SM fermion-pair ($`f\overline{f}`$) and the lightest Higgs–boson pair ($`hh`$) final states were given in Ref. , and for the $`WW`$ and $`ZZ`$ final states in Ref. . The annihilation into $`f\overline{f}`$ is often most important. However, other final states can also play a considerable role, depending on the model.

We have derived a full set of exact analytic expressions for the cross section of the neutralino pair-annihilation in the minimal supersymmetric standard model (MSSM). We have made no simplifying assumptions about the degeneracy of the left- and right–sfermion masses, and we included all tree-level final states and all intermediate states. We have kept finite widths in $`s`$-channel resonances. While interference terms and usually sub-dominant channels are sometimes neglected in the literature, in some cases they may give significant contributions. In particular, the channels involving gauge–Higgs boson may in some cases be sizeable or even dominant. One example where this is the case is the $`W^\pm H^{}`$ final state, and in this paper we give expressions for the cross section for the chargino exchange contribution. A full set of all the expressions for all the final states will be presented elsewhere .

As a check, we have performed a numerical comparison of our cross section with the results obtained by using the recently released code DarkSusy . We have found, for the same values of input parameters, an impressive agreement, at the level of a few per cent, for all the annihilation channels, which we find reassuring.

We have also made a comparison of the relic abundance computed using our exact formulae with the one obtained using our expansion formulae. So far, no such detailed analysis has been presented in the literature in the case of the MSSM. In Ref.  an analogous comparison was made in the context of the highly constrained minimal supergravity scheme and only relatively close to resonances. We confirm that the expansion gives highly inaccurate results near resonances and new thresholds. We also show that very far from such cases the error is typically rather small, of the order of a few per cent. However, we find that, because of the existence of several resonances ($`Z`$ and the Higgs bosons), the expansion produces large errors, compared to an exact treatment, over a sizeable range of $`m_\chi `$, even of a several tens of GeV. We therefore conclude that the widely used method of expansion may lead to significant errors in a sizeable fraction of the neutralino mass which is typically expected not to significantly exceed $`200\text{GeV}`$ by various naturalness criteria .

In the current version of our code we have used a popular approximate iterative method of determining the freeze-out point. We have examined the accuracy of the procedure and compared our results with the ones obtained using DarkSusy. We have found that the method in general works very well when both exact and expansion methods are used, except in the latter case near resonances and thresholds where it may even break down. The iterative procedure can be safely applied in such special cases, provided exact expressions for the cross section are used and much care is paid to properly integrating the thermal average numerically.

## 2 Calculation of the Relic Density

The relic abundance of some stable species $`\chi `$ is given by $`\mathrm{\Omega }_\chi \rho _\chi /\rho _{crit}`$, where $`\rho _\chi `$ is the relic’s (mass) density and $`\rho _{crit}3H_0^2/8\pi G=1.9\times 10^{29}(h^2)g/cm^3`$ is the critical density. It can be computed by solving the Boltzmann equation, which describes the time evolution of the co-moving number density $`n_\chi `$ in the expanding Universe,111 For a review of calculations of the relic density, see, e.g., Refs. .

$`{\displaystyle \frac{dn_\chi }{dt}}`$ $`=`$ $`3Hn_\chi \sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}\left(n_\chi ^2(n_\chi ^{\mathrm{eq}})^2\right),`$ (2)
where $`n_\chi ^{\mathrm{eq}}`$ is the number density that the species would have in thermal equilibrium, $`H(T)`$ is the Hubble expansion rate, $`\sigma (\chi \chi \mathrm{all})`$ denotes the cross section of the species annihilation into ordinary particles, $`v_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ is a so-called Møller velocity  which is the relative velocity of the annihilating particles, and $`\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ represents the thermal average of $`\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ which will be given below. In the early Universe, the species $`\chi `$ were initially in thermal equilibrium, $`n_\chi `$ $`=n_\chi ^{\mathrm{eq}}`$. When their typical interaction rate $`\mathrm{\Gamma }_\chi `$ became less than the Hubble parameter, $`\mathrm{\Gamma }_\chi \text{ }<H`$, the annihilation process froze out. Since then their number in a co-moving volume has remained basically constant.

In order to calculate the relic density precisely enough, one has to consider the following ingredients: (i) compute all the contributions to the annihilation cross section, (ii) use an exact formula for the thermal average, (iii) rigorously solve the Boltzmann equation, and (iv) include co-annihilation. (This becomes important when the mass of the next LSP is nearly degenerate with the LSP mass .) In this work, we will concentrate on points (i) and (ii) but will also examine to some extent point (iii). In particular, we will compare the usually used expansion formulae with the exact ones. We will not consider here the effect of co-annihilation.

An often used, approximate, although in general quite accurate (see later), solution to the Boltzmann equation is based on solving iteratively the equation

$`x_f^1`$ $`=`$ $`\mathrm{ln}\left({\displaystyle \frac{m_\chi }{2\pi ^3}}\sqrt{{\displaystyle \frac{45}{2g_{}G_N}}}\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}(x_f)x_f^{1/2}\right),`$ (3)
where $`g_{}`$ represents the effective number of degrees of freedom at freeze-out ($`\sqrt{g_{}}9`$) and $`G_N`$ is the Newton constant. Typically one finds that the freeze-out point $`x_fT_f/m_\chi `$ is roughly given by $`1/20`$. The relic density $`\rho _\chi =m_\chi n_\chi `$ at present is given by

$`\rho _\chi `$ $`=`$ $`{\displaystyle \frac{1.66}{M_{\mathrm{Pl}}}}\left({\displaystyle \frac{T_\chi }{T_\gamma }}\right)^3T_\gamma ^3\sqrt{g_{}}{\displaystyle \frac{1}{J(x_f)}},`$ (4)
where $`M_{\mathrm{Pl}}`$ denotes the Planck mass, $`T_\chi `$ and $`T_\gamma `$ are the present temperatures of the neutralino and the photon, respectively. The suppression factor $`(T_\chi /T_\gamma )^3`$ $``$ $`1/20`$ follows from entropy conservation in a comoving volume . Finally, $`J(x_f)`$ is given by

$$J(x_f)_0^{x_f}𝑑x\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}(x),$$
(5)
where $`x=T/m_\chi `$, as defined earlier.

Below we will concentrate on computing $`J(x_f)`$ and in particular on the thermally averaged annihilation cross section times velocity $`\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$. A proper definition which includes separate thermal baths for both annihilating particles is given by

$$\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}(T)=\frac{d^3p_1d^3p_2\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}e^{E_1/T}e^{E_2/T}}{d^3p_1d^3p_2e^{E_1/T}e^{E_2/T}}$$
(6)
where $`p_1=(E_1,𝐩_1)`$ and $`p_2=(E_2,𝐩_2)`$ are the 4-momenta of the two colliding particles, and $`T`$ is the temperature of the bath.222 Note that in many cases one uses another definition of $`<\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}>`$ which involves a single thermal bath for both neutralinos. Compare, e.g., Refs. .  The above expression can be reduced to a one-dimensional integral which can be written in a Lorentz-invariant form as

$`\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}(T)`$ $`=`$ $`{\displaystyle \frac{1}{8m_\chi ^4TK_2^2(m_\chi /T)}}{\displaystyle _{4m_\chi ^2}^{\mathrm{}}}𝑑s\sigma (s)(s4m_\chi ^2)\sqrt{s}K_1\left({\displaystyle \frac{\sqrt{s}}{T}}\right),`$ (7)
where $`s=(p_1+p_2)^2`$ is a usual Mandelstam variable and $`K_i`$ denotes the modified Bessel function of order $`i`$.

In computing the cross section $`\sigma (s)`$, it is convenient to introduce a Lorentz-invariant function $`w(s)`$

$$w(s)=\frac{1}{4}𝑑\mathrm{LIPS}|𝒜(\chi \chi \mathrm{all})|^2$$
(8)
where $`|𝒜(\chi \chi \mathrm{all})|^2`$ denotes the absolute square of the reduced matrix element for the annihilation of two $`\chi `$ particles, averaged over initial spins and summed over final spins. The function $`w(s)`$ is related to the annihilation cross section via

$$w(s)=\frac{1}{2}\sqrt{s(s4m_\chi ^2)}\sigma (s).$$
(9)
For a general annihilation process into a two-body final state $`\chi \chi f_1f_2`$, the function $`w(s)`$ can be written as

$`w(s)`$ $`=`$ $`{\displaystyle \frac{1}{32\pi }}{\displaystyle \left[\theta \left(s(m_{f_1}+m_{f_2})^2\right)\beta _f(s,m_{f_1},m_{f_2})\stackrel{~}{w}_{f_1f_2}(s)\right]},`$ (10)
where the summation goes over all possible channels, and

$`\stackrel{~}{w}_{f_1f_2}(s){\displaystyle \frac{1}{8\pi }}{\displaystyle 𝑑\mathrm{\Omega }|𝒜(\chi \chi f_1f_2)|^2},`$ (11)
where

$`\beta _f(s,m_{f_1},m_{f_2})\left[1{\displaystyle \frac{(m_{f_1}+m_{f_2})^2}{s}}\right]^{1/2}\left[1{\displaystyle \frac{(m_{f_1}m_{f_2})^2}{s}}\right]^{1/2}.`$ (12)
Computation of $`w(s)`$ is in general rather involved and final analytic expressions exhibit considerable complexity, as will be illustrated in the next Section. This is one reason why in most of the literature one uses the well-known approximation in terms of expansion in powers of $`x`$, $`\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}a+bx`$. Using the definition (6), the expansion reads

$$<\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}>=\frac{1}{m_\chi ^2}\left[w\frac{3}{2}\left(2ww^{}\right)x+𝒪(x^2)\right]_{s=4m_\chi ^2}a+bx+𝒪(x^2)$$
(13)
where $`w^{}(s)`$ denotes $`dw(s)/d(s/4m_\chi ^2)`$.

Note that in case of the expansion of $`<\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}>`$ defined in terms of a single thermal bath for both neutralinos, the corresponding coefficient $`a^{}`$ is the same as Eq. (13) while $`b^{}=b+\frac{3}{2}a`$ . For more details, see, e.g., Ref. .

The expansion method is widely used not only because it is computationally somewhat less involved but also because it is expected to be relatively accurate. This can seen by examining the behavior of the integrand in Eq. (7). For a massive particle like the neutralino for which $`T\text{ }<m_\chi /20`$ and $`\sqrt{s}2m_\chi `$, the argument of the function $`K_1`$ is much larger than unity. Since $`K_1(y)\sqrt{\pi /2y}e^y`$ at $`y1`$, the thermal average (7) can be written as a convolution of the cross section with a function $`f(s)`$:

$`\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$ $``$ $`{\displaystyle _{4m_\chi ^2}^{\mathrm{}}}𝑑s\sigma (s)f(s),`$ (14)
where $`f(s)`$ has an exponential suppression factor, $`f(s)e^{\sqrt{s}/T}`$. Thus one would expect that the usual expansion in powers of $`x`$ should converge quickly .

One should however use this argument with some care. First, the annihilation cross section may change rapidly with $`s`$. It is well-known that this happens, e.g., near $`s`$-channel resonances and thresholds of new final states. In such cases the expansion method produces large errors . Far away from such singular points the expansion is expected to give accurate results but, to our knowledge, this has never been explicitly verified in the literature. Finally, one would want to examine actually how far from resonances one can start relying on the usual approximation. As we will see later, the range of $`m_\chi `$ where the expansion fails to be reasonably accurate can be actually quite large.

## 3 Analytic Results

As stated in the Introduction, we have derived a full set of exact, analytic expressions for the cross sections for the neutralino pair-annihilation processes into all allowed (tree-level) final states in the general MSSM. We have included all contributing diagrams as well as all interference terms and kept finite widths of all $`s`$-channel resonances. We have made no simplifying assumptions about sfermion masses although we assumed that there are no CP violating phases in SUSY parameters. Here we will give one analytic example of our exact expressions. A full set of results for all the final states will be given elsewhere .

Next, starting from the exact expressions we have computed the coefficients $`a`$ and $`b`$ using Eq. (13) for all the partial annihilation channels. In the literature there exist several analytic formulae for the expansion coefficients, including , but, due to different conventions and complexity of expressions, comparison is not always doable. We have checked our results for the $`a`$-coefficients in appropriate limits against published results and agreed in all cases. In the next Section we will present a numerical comparison of the approximation in terms of the usual expansion with our exact results.

Let us begin by introducing the relevant MSSM parameters.333We follow the convention of Ref. . The lightest neutralino is a mass eigenstate given by a linear combination of two neutral gauginos and two neutral higgsinos

$`\chi `$ $`=`$ $`N_{11}\stackrel{~}{B}+N_{12}\stackrel{~}{W}_3+N_{13}\stackrel{~}{H}_1^0+N_{14}\stackrel{~}{H}_2^0.`$ (15)
The neutralino mass matrix is determined by the $`U(1)_Y`$ and $`SU(2)_L`$ gaugino mass parameters $`M_1`$ and $`M_2`$ (and we impose the usual GUT relation $`M_1=\frac{5}{3}\mathrm{tan}^2\theta _WM_2`$), the Higgs/higgsino mass parameter $`\mu `$ and $`\mathrm{tan}\beta =v_2/v_1`$ which is the ratio of the vacuum expectation values of the two neutral Higgs fields.

There are two neutral scalar Higgs bosons $`h`$ and $`H`$ and one pseudoscalar $`A`$ plus a pair of charged Higgs $`H^\pm `$. (We will typically suppress the Higgs charge assignment except where this may lead to ambiguities.) We use expressions given in Ref.  for computing radiatively corrected Higgs masses.

Other relevant parameters which determine the masses of scalars and various couplings are the squark soft mass parameters $`m_Q`$, $`m_U`$ and $`m_D`$, the slepton soft mass parameters $`m_L`$ and $`m_E`$, and the pseudoscalar mass $`m_A`$. We also include the trilinear terms $`A_i`$ ($`i=t,b,\tau `$) of the third generation which are important in determining the masses and couplings for the stop, sbottom and stau sfermions, respectively.

There are a number of final states into which the neutralino can pair-annihilate. They include: $`\chi \chi f\overline{f}`$, $`WW`$, $`ZZ`$, $`Zh`$, $`ZH`$, $`ZA`$, $`W^\pm H^{}`$, $`AA`$, $`AH`$, $`Ah`$, $`HH`$, $`Hh`$, $`hh`$ and $`H^+H^{}`$. Among these, in the gaugino region, fermion pair $`f\overline{f}`$ final states usually give dominant contributions, unless the exchanged sfermion masses $`m_{\stackrel{~}{f}}`$ are large and in some special cases, as discussed earlier. These final states are also always kinematically allowed (except for $`t\overline{t}`$) but, due to the $`s`$-wave suppression , their cross section is proportional to $`m_f^2/m_W^2`$. For the higgsino-like neutralinos the states $`ZZ`$ and $`WW`$ become very important once kinematically allowed . However, these final states can give significant contributions even in the case of the gaugino-like neutralino. This is because $`WW`$ and $`ZZ`$ are not $`s`$-wave suppressed, unlike the channels $`hh`$, $`Hh`$, $`HH`$, $`AA`$, $`H^+H^{}`$ and $`ZA`$ . What is interesting is that the cross sections for the other $`s`$-wave unsuppressed states $`W^\pm H^{}`$, $`Zh`$, $`ZH`$, $`Ah`$ and $`AH`$, once kinematically allowed, can be even larger than those of $`WW`$ and $`ZZ`$. In fact, they can even dominate over those of $`f\overline{f}`$. These points will be illustrated with numerical examples in Section 4.

Here, we present a set of expressions for a $`t`$\- and $`u`$-channel chargino ($`\chi _{1,2}^\pm `$) exchange contribution to the process $`\chi \chi `$ $``$ $`W^+H^{}`$. This contribution is often dominant. (There are also $`s`$-channel diagrams involving $`h`$, $`H`$ and $`A`$ exchange. Since typically $`m_hm_Hm_Am_{H^\pm }`$, the Higgs resonances will be outside of the kinematically allowed region.) The function $`\stackrel{~}{w}`$, Eq. (11) in this case reads

$`\stackrel{~}{w}_{W^+H^{}}^{\chi ^\pm }`$ $`=`$ $`{\displaystyle \frac{1}{m_W^2}}{\displaystyle \underset{i,j=1}{\overset{2}{}}}\left[m_{\chi _i^\pm }m_{\chi _j^\pm }I_{ij}^{(2)}+m_\chi m_{\chi _i^\pm }I_{ij}^{(1)}+I_{ij}^{(0)}\right],`$ (16)
where

$`I_{ij}^{(2)}`$ $`=`$ $`(C_{+i}^{HW^{}}C_{+j}^{HW}+D_{+i}^{HW^{}}D_{+j}^{HW})`$ (17)
$`\times \left[𝒯_2(sm_W^2m_{H^\pm }^22m_\chi ^2)𝒯_1+G_{WH}^{T(1)}𝒯_0+6m_W^2m_\chi ^2𝒴_0\right]`$
$`+(C_{+i}^{HW^{}}C_{+j}^{HW}D_{+i}^{HW^{}}D_{+j}^{HW})\left[\mathrm{\hspace{0.17em}6}m_W^2m_\chi ^2𝒯_0𝒴_2+G_{WH}^{Y(1)}𝒴_0\right],`$
$`I_{ij}^{(1)}`$ $`=`$ $`\mathrm{Re}(C_{+i}^{HW^{}}C_j^{HW}+D_{+i}^{HW^{}}D_j^{HW})[2𝒯_2+2(2m_\chi ^2m_W^2)𝒯_1`$ (18)
$`2(m_\chi ^2m_W^2)(m_\chi ^2+2m_W^2)𝒯_0+2m_W^2𝒴_1G_{WH}^{Y(2)}𝒴_0]`$
$`+\mathrm{Re}(C_{+i}^{HW^{}}C_j^{HW}D_{+i}^{HW^{}}D_j^{HW})[\mathrm{\hspace{0.17em}6}m_W^2𝒯_1+6m_W^2(m_\chi ^2m_{H^\pm }^2)𝒯_0`$
$`2𝒴_2(sm_{H^\pm }^24m_W^2)𝒴_1G_{WH}^{Y(3)}𝒴_0],`$
$`I_{ij}^{(0)}`$ $`=`$ $`(C_i^{HW^{}}C_j^{HW}+D_i^{HW^{}}D_j^{HW})[(sm_\chi ^22m_W^2)𝒯_2G_{WH}^{T(2)}𝒯_1`$ (19)
$`G_{WH}^{T(3)}𝒯_0+(s2m_W^2)𝒴_2+G_{WH}^{Y(4)}𝒴_0]`$
$`+(C_i^{HW^{}}C_j^{HW}D_i^{HW^{}}D_j^{HW})\left[\mathrm{\hspace{0.17em}6}m_W^2m_\chi ^2𝒯_1m_\chi ^2𝒴_2+G_{WH}^{Y(5)}𝒴_0\right].`$
Obviously the total contribution to $`w(s)`$ from $`W^\pm H^{}`$ is twice that from $`W^+H^{}`$. The symbols $`C_{\pm i}^{HW}`$ and $`D_{\pm i}^{HW}`$ in the above equations are combinations of several coupling constants:

$`C_{\pm i}^{HW}`$ $``$ $`C_S^{\chi _i^+\chi H^{}}C_V^{\chi _i^+\chi W^{^{}}}\pm C_P^{\chi _i^+\chi H^{}}C_A^{\chi _i^+\chi W^{^{}}},`$ (20)
$`D_{\pm i}^{HW}`$ $``$ $`C_S^{\chi _i^+\chi H^{}}C_A^{\chi _i^+\chi W^{^{}}}\pm C_P^{\chi _i^+\chi H^{}}C_V^{\chi _i^+\chi W^{^{}}}.`$ (21)
The relevant couplings describe neutralino-chargino-$`W^{}`$/$`H^{}`$ interaction terms in the Lagrangian as

$`_{\mathrm{int}}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{2}{}}}\overline{\chi _i^{}}\gamma ^\mu \left(C_V^{\chi _i^+\chi W^{}}C_A^{\chi _i^+\chi W^{}}\gamma _5\right)\chi W_\mu ^{}`$ (22)
$`+{\displaystyle \underset{i=1}{\overset{2}{}}}\overline{\chi _i^{}}\left(C_S^{\chi _i^+\chi H^{}}C_P^{\chi _i^+\chi H^{}}\gamma _5\right)\chi H^{}+\mathrm{h}.\mathrm{c}.`$
in the convention of Ref. . The two pairs of charginos are denoted by $`\chi _i^\pm `$ ($`i=1,2`$) and $`m_{\chi _i^\pm }`$ represent their masses. The functions $`𝒯_k`$ $`=`$ $`𝒯_k`$($`s`$, $`m_\chi ^2`$, $`m_{H^\pm }^2`$, $`m_W^2`$, $`m_{\chi _i^\pm }^2`$, $`m_{\chi _j^\pm }^2`$), $`𝒴_k`$ $`=`$ $`𝒴_k`$($`s`$, $`m_\chi ^2`$, $`m_{H^\pm }^2`$, $`m_W^2`$, $`m_{\chi _i^\pm }^2`$, $`m_{\chi _j^\pm }^2`$), $`k`$ $`=`$ $`0,1,2`$, and the functions $`G_{WH}^{T(13)}`$, $`G_{WH}^{Y(15)}`$ are listed in the Appendix.

## 4 Numerical Results

In this Section we will illustrate the points made above by presenting several numerical results. We will concentrate on the gaugino-like neutralino. We remind the reader that a nearly pure gaugino (bino) is the most natural choice for the CDM both in the MSSM  and in the CMSSM where in most cases the LSP comes out to be a nearly pure bino .

We begin by plotting in Figs. 1 and 2 the function $`J(x_f)`$, Eq. (5), computed using exact and expansion methods as a function of $`m_\chi `$. We take $`\mathrm{tan}\beta =10`$, $`m_0m_Q=m_U=m_D=m_L=m_E=500\text{GeV}`$, $`m_A=400\text{GeV}`$, $`A_t=A_b=600\text{GeV}`$ and $`\mu =500\text{GeV}`$. To obtain sparticle and Higgs mass spectra we have used the code DarkSusy. While several analytic expressions and numerical codes (e.g., SUSPECT ) are available in the literature, here we have used DarkSusy to facilitate the comparison of our results for the cross sections with those produced by the package. For the values given above we obtain the following Higgs masses: $`m_h=117.3\text{GeV}`$, $`m_H=400.5\text{GeV}`$ and $`m_{H^\pm }=407.9\text{GeV}`$.

We show both the total and the individual contributions from all the final states. The solid and dotted lines represent the exact and the expansion results, respectively. For the sake of comparison of both methods, for now we do not impose experimental constraints. Note that $`m_\chi `$ increases with increasing $`M_2`$ while we keep $`\mu `$ fixed at a rather large value. This means that in these Figures (and also the following ones below), the neutralino is mostly a gaugino. Even at $`m_\chi =400\text{GeV}`$, $`N_{11}^2+N_{12}^2>0.9`$. (However, between $`400\text{GeV}`$ and $`500\text{GeV}`$ the LSP turns into a nearly pure higgsino.) In this region the effect of co-annihilation with the lightest chargino and the next lightest neutralino, which we do not consider here, is not important and can be safely neglected.

For each final state, we compute $`J(x_f)`$, Eq. (5), numerically via Eq. (7) by using our exact expressions for the cross sections ($`J_{\mathrm{exact}}(x_f)`$). Next, we use the cross sections to compute the first two expansion coefficients $`a`$ and $`b`$, as in Eq. (13), and to compute $`J_{\mathrm{exp}}(x_f)=ax_f+\frac{1}{2}bx_f^2`$. For consistency, in both cases we determine the freeze-out point $`x_f`$ by employing the iterative procedure (3). We will come back to the accuracy of the iterative procedure below.

In the behavior of the total $`J(x_f)`$ one can clearly recognize three peaks which are due to the $`s`$-channel $`Z`$, $`h`$ and $`A/H`$ exchange contributions. In particular, we can see that the $`A`$-resonance is very broad. (Due to a near mass degeneracy $`m_Hm_A`$, the usually narrower $`H`$-pole is ‘burried” underneath the $`A`$-pole.) This is caused by the amplification of the coupling $`Ab\overline{b}\mathrm{tan}\beta `$ at larger $`\mathrm{tan}\beta `$ and the fact that, in the resonance region, the amplitude-square for the $`A`$-exchange is proportional to $`s`$, unlike the case of $`h`$ and $`H`$ where it goes like $`s4m_\chi ^20`$ at $`\sqrt{s}2m_\chi `$. The $`Z`$-pole is also clearly visible in all the partial $`f\overline{f}`$ final states while the Higgs poles do not contribute to the $`\nu \overline{\nu }`$ final state since neutrinos do not couple to Higgs bosons. The $`H`$-pole also gives a visible enhancement to the $`WW`$, $`ZZ`$ and $`hh`$ final states but not to $`W^\pm H^{}`$, $`ZA`$, $`HH`$, $`Hh`$, $`AA`$ and $`H^+H^{}`$ which do not become kinematically allowed until $`m_\chi >m_H/2`$. Likewise, the $`A`$-pole shows a resonance in the $`Zh`$ final state but not in $`W^\pm H^{}`$, $`ZH`$, $`AH`$ and $`Ah`$.

For the pseudoscalar Higgs exchange in $`\chi \chi f\overline{f}`$, the expansion method gives a discontinuity very close to the $`A`$-pole, $`m_\chi m_A/2=200\text{GeV}`$. For this channel, which is clearly dominant in the resonance region, the $`b`$-coefficient is large and negative and the expansion method and iterarative procedure break down. This can be clearly seen in all the $`f\overline{f}`$ windows (including $`\nu \nu `$) and also in the total value of $`J_{\mathrm{exp}}(x_f)`$. In the close vicinity of the poles, one can see the well-known large difference  (note a log scale) between $`J_{\mathrm{exact}}(x_f)`$ and $`J_{\mathrm{exp}}(x_f)`$ while away from the poles, both contributions are similar.

We can also notice that $`J_{\mathrm{exact}}(x_f)`$ and $`J_{\mathrm{exp}}(x_f)`$ show a significantly different enhancement in the total values and in the up-type quark-pair final state in the region $`m_\chi \text{ }>m_t=175\text{GeV}`$ where the $`t\overline{t}`$ final state becomes kinematically allowed. In particular, $`J_{\mathrm{exact}}(x_f)`$ increases gradually starting from $`m_\chi `$ some $`10\text{GeV}`$ below $`m_t`$ while $`J_{\mathrm{exp}}(x_f)`$ exhibits a sharp increase only at $`m_\chi =m_t`$. This is because the exact treatment takes into account the thermal distribution of momenta of the annihilating neutralinos, while the expansion method neglects it . The same effect is visible near the thresholds of all the final states which only become open at some non-zero value of $`m_\chi `$, especially of those involving one or two Higgs bosons.

As we stated in Section 3, once kinematically allowed, the contributions from the $`W^\pm H^{}`$ final state can be comparable, or even larger (if sfermion masses are increased) than the $`f\overline{f}`$ final state. In fact, the $`ZH`$, $`Ah`$ and, at much larger $`m_\chi `$ also $`AH`$, final state contributions are not much smaller. We can also see that the channels $`ZZ`$ and $`WW`$, as well as $`Zh`$, despite not being $`s`$-wave suppressed, are usually subdominant, except near $`m_\chi m_A/2`$ because of the pseudoscalar Higgs exchange contributing to those channels, and when the higgsino admixture of the neutralino increases which is the case with increasing $`m_\chi `$ while keeping $`\mu `$ fixed.

As a way of verifying our results, we have performed a numerical comparison of our exact results for the cross section (more precisely of the function $`w(s)`$) with the same quantity computed using DarkSusy. We have found an excellent agreement at the level of a few per cent for all the tree-level final states when we fixed the same values of input parameters. This gives us some confidence that the (highly complex) expressions that we have derived (and also those in DarkSusy), are correct.

As mentioned above, in the current version of our numerical code we have used the iterative method (3) which is widely used in the literature. We applied it to both the exact formula for $`\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}`$, Eq. (7), and to the usual approximation (13). The results are presented in Fig. 3 as a solid and dashed line, respectively. For comparison, we also display $`x_f`$ using DarkSusy (dotted line) where a different numerical procedure is used. It is clear that there is a general agreement among the three different procedures and also that $`x_f`$ can be only roughly approximated by the usually quoted value $`0.05`$. Away from the resonances and thresholds the expansion method can be safely used in determining $`x_f`$ while near such special points this is not the case. In particular, near the $`A`$-pole the expansion gives $`\sigma v_{\mathrm{M}\mathrm{ø}\mathrm{l}}<0`$ and the iterative procedure breaks down. In contrast, the exact numerical integration (6) provides reliable results both near and away from special points. However, we have carefully verified that this is the case only when the integration in Eq. (6) is performed with high enough numerical precision.

In Fig. 4a, we plot the ratio $`\mathrm{\Omega }_{\mathrm{exp}}/\mathrm{\Omega }_{\mathrm{exact}}`$ versus $`m_\chi `$ for the same set of parameters as before. We denote by $`\mathrm{\Omega }_{\mathrm{exact}}`$ and $`\mathrm{\Omega }_{\mathrm{exp}}`$ the neutralino relic density calculated with the exact (7) and approximate (13) formulae for the thermal average, respectively. The relic abundance in both cases is computed by solving Eq. (3) iteratively and using Eq. (4). From now on we impose all applicable experimental constraints, including the lightest Higgs mass constraint $`m_h>113.5\text{GeV}`$  and a lower bound $`104\text{GeV}`$ for chargino masses. In the case of $`bs\gamma `$ we have followed DarkSusy in using a somewhat old range $`1.0\times 10^4`$ $`<`$ $`B(bs\gamma )`$ $`<`$ $`4.0\times 10^4`$  and have not yet included new SUSY contributions at large $`\mathrm{tan}\beta `$ that have been recently pointed out . We feel that this is sufficient for the purpose of this analysis. The solid line represents the experimentally allowed region while the dotted line corresponds to the low $`m_\chi `$ range excluded by collider experiments.

It is evident that the ratio $`\mathrm{\Omega }_{\mathrm{exp}}/\mathrm{\Omega }_{\mathrm{exact}}`$ varies considerably with $`m_\chi `$. Near resonances it is very large as discussed above. If we conservatively demand a $`10\%`$ accuracy in the relic density computations then the approximate method fails to satisfy this criterion over a broad range $`160\text{GeV}\text{ }<m_\chi \text{ }<220\text{GeV}`$, which is a significant fraction of the neutralino mass range. (In fact, various naturalness constraints, even if only indicative, give a rough constraint $`m_\chi \text{ }<200\text{GeV}`$ .) For larger values of $`\mathrm{tan}\beta `$ the range of $`m_\chi `$ where $`0.9<\mathrm{\Omega }_{\mathrm{exp}}/\mathrm{\Omega }_{\mathrm{exact}}<1.1`$ is not satisfied becomes even larger.

It is also clear that even further away from resonances the ratio $`\mathrm{\Omega }_{\mathrm{exp}}/\mathrm{\Omega }_{\mathrm{exact}}`$ exhibits a rather complex behavior. For $`220\text{GeV}\text{ }<m_\chi \text{ }<260\text{GeV}`$ the new channels $`W^\pm H^{}`$, $`ZA`$ and $`ZH^0`$ successively kick in, and become actually even somewhat more important than $`f\overline{f}`$. Since, as we mentioned above, the expansion method does not work properly near new channels thresholds , the ratio varies quite rapidly over the mentioned range of $`m_\chi `$, although not as much as near resonances. Finally, far away from resonances and thresholds, the expansion method seems to work remarkably well.

In Fig. 4b we show the relic abundance $`\mathrm{\Omega }_\chi h^2`$ as a function of $`m_\chi `$. Again, the solid (dotted) line represents the range of $`m_\chi `$ allowed (excluded) by the experimental constraints. The two horizontal dashed lines delimit the cosmologically expected range (1) of $`\mathrm{\Omega }_\chi h^2`$. Two features should be mentioned. Over the range of smaller $`m_\chi `$, in particular near resonances, the relic abundance varies rapidly. On the other hand, at larger $`m_\chi `$, where new channels involving the “heavy” Higgs bosons $`H^\pm `$, $`A`$ and $`H`$ open up and can become important, $`\mathrm{\Omega }_\chi h^2`$ exhibits a wide plateau. In both regions it is clearly important to compute $`\mathrm{\Omega }_\chi h^2`$ accurately if one wants to achieve an accuracy expected from future determinations of $`\mathrm{\Omega }_{\mathrm{CDM}}h^2`$. In particular, it is clear from Fig. 4b that at larger $`m_\chi `$ a relatively small difference of a few per cent in determining $`\mathrm{\Omega }_\chi h^2`$ could exclude on cosmological grounds the range of $`m_\chi \text{ }>280\text{GeV}`$ which in Fig. 4b is (barely) allowed. This again illustrates the uncertainties involved in computing $`\mathrm{\Omega }_\chi h^2`$ and in determining the cosmologically allowed regions of SUSY parameters.

## 5 Conclusions

While we are certainly still some time away from a “precision era” for measuring cosmological parameters, much effort has been focussed on improving theoretical calculations of the relic abundance of the neutralino. We have derived a full set of analytic expressions for the neutralino annihilation cross sections into all tree-level final states. In this paper, we have presented one example of such an expression for the neutralino annihilation into $`W^+H^{}`$ which can be even dominant. We have also computed the first two terms in the expansion of the exact formulae in powers of $`x`$ and we have compared the integral of the thermally averaged product of the cross section and velocity and the relic abundance computed both ways. We have confirmed the well-known inaccuracy of the approximate method near reasonances and threshold and showed that, far away from such cases the expansion method works well. However, due to the presence of several resonances and thresholds in the MSSM, the expansion method may lead to sizeable errors over a relatively large range of the neutralino mass of up to several tens of GeV. We have also demonstrated that the iterative way of computing the freeze-out point works rather well for both the exact and expansion methods of computing the thermal average but not in the latter case near resonances and thresholds.

## Acknowledgements

T.N. would like to thank the CERN Theory Division for its kind hospitality.

## Appendix

Here we give expressions for the auxiliary functions used in the text. First, we define

$`D(s,x,y_1,y_2)`$ $``$ $`x+{\displaystyle \frac{y_1+y_2}{2}}{\displaystyle \frac{s}{2}},`$ (23)
$`F(s,x,y_1,y_2)`$ $``$ $`{\displaystyle \frac{1}{2}}\sqrt{s4x}\sqrt{s(\sqrt{y_1}+\sqrt{y_2})^2}\sqrt{1{\displaystyle \frac{(\sqrt{y_1}\sqrt{y_2})^2}{s}}}.`$ (24)
If we define $`DD(s,x,y_1,y_2)`$, $`FF(s,x,y_1,y_2)`$, $`t_\pm (s,x,y_1,y_2)D\pm F`$ and $`(𝒯_i,𝒴_i)`$ $``$ $`(𝒯_i,𝒴_i)(s,x,y_1,y_2,z_1,z_2)`$ ($`i`$ $`=`$ $`0,1,2`$), then

$`(s,x,y_1,y_2,z)`$ $`=`$ $`{\displaystyle \frac{1}{2F}}\mathrm{ln}\left|{\displaystyle \frac{t_+(s,x,y_1,y_2)z}{t_{}(s,x,y_1,y_2)z}}\right|,`$ (25)
and

$`𝒯_0`$ $`=`$ $`{\displaystyle \frac{1}{z_1z_2}}[(s,x,y_1,y_2,z_1)(s,x,y_1,y_2,z_2)],`$
$`𝒯_1`$ $`=`$ $`{\displaystyle \frac{1}{z_1z_2}}[z_1(s,x,y_1,y_2,z_1)z_2(s,x,y_1,y_2,z_2)],`$
$`𝒯_2`$ $`=`$ $`1+{\displaystyle \frac{1}{z_1z_2}}[z_1^2(s,x,y_1,y_2,z_1)z_2^2(s,x,y_1,y_2,z_2)],`$
$`𝒴_0`$ $`=`$ $`{\displaystyle \frac{1}{z_1+z_22D}}[(s,x,y_1,y_2,z_1)+(s,x,y_1,y_2,z_2)],`$
$`𝒴_1`$ $`=`$ $`{\displaystyle \frac{1}{z_1+z_22D}}[2(z_1D)(s,x,y_1,y_2,z_1)2(z_2D)(s,x,y_1,y_2,z_2)],`$
$`𝒴_2`$ $`=`$ $`1+{\displaystyle \frac{1}{z_1+z_22D}}[z_1(z_12D)(s,x,y_1,y_2,z_1)`$
$`+z_2(z_22D)(s,x,y_1,y_2,z_2)].`$
The expressions for $`G_{WH}^{T(i)}`$ and $`G_{WH}^{Y(j)}`$, where $`i=1,2,3`$ and $`j=1,\mathrm{},5`$, in Eqs. (17-19) are given by

$`G_{WH}^{T(1)}`$ $``$ $`s(m_\chi ^2+2m_W^2)m_\chi ^4m_\chi ^2(m_{H^\pm }^2+3m_W^2)m_W^2m_{H^\pm }^2,`$
$`G_{WH}^{T(2)}`$ $``$ $`s(m_\chi ^2+2m_W^2)2m_\chi ^4+m_\chi ^2(m_{H^\pm }^2m_W^2)2m_W^2(m_{H^\pm }^2+m_W^2),`$
$`G_{WH}^{T(3)}`$ $``$ $`(m_\chi ^2m_{H^\pm }^2)(m_\chi ^2m_W^2)(m_\chi ^2+2m_W^2),`$
$`G_{WH}^{Y(1)}`$ $``$ $`G_{WH}^{T(1)},`$
$`G_{WH}^{Y(2)}`$ $``$ $`3sm_W^212m_W^2m_\chi ^2+3m_W^2m_{H^\pm }^23m_W^4,`$
$`G_{WH}^{Y(3)}`$ $``$ $`s^2s(2m_\chi ^2+2m_{H^\pm }^2+3m_W^2)+2m_\chi ^4+(2m_\chi ^2+m_{H^\pm }^2)(m_{H^\pm }^2+3m_W^2)2m_W^4,`$
$`G_{WH}^{Y(4)}`$ $``$ $`s(m_\chi ^22m_W^2)(m_\chi ^2m_{H^\pm }^2)+4m_W^2m_\chi ^4`$
$`+m_\chi ^2(m_{H^\pm }^45m_W^2m_{H^\pm }^2+2m_W^4)2m_W^4m_{H^\pm }^2,`$
$`G_{WH}^{Y(5)}`$ $``$ $`sm_\chi ^2(m_\chi ^2m_W^2)m_\chi ^6m_\chi ^4(m_{H^\pm }^2+3m_W^2)+m_\chi ^2m_W^2(2m_{H^\pm }^2+3m_W^2).`$