Towards an accurate calculation of the neutralino relic density

We re-analyze the neutralino relic density in the minimal supersymmetric standard model. We use exact expressions for the neutralino annihilation cross section into all tree-level final states, including all contributions and interference terms. We demonstrate that several final states may give comparable contributions to the relic density, including the gauge-Higgs and Higgs-Higgs boson channels. This illustrates the importance of performing a complete calculation. We consistently compare the exact results with those of the usual expansion method. We confirm that the expansion gives highly inaccurate results very near resonances and new thresholds and show that very far from such cases the discrepancy is typically only of the order of a few per cent. However, we demonstrate that, near the pseudoscalar Higgs resonance, the expansion can give a sizeable error (of more than 10%) over a significant range of the neutralino mass which can be as large as several tens of GeV. We also perform several related checks and comparisons. In particular, we show that the often employed approximate iterative procedure of computing the neutralino freeze-out temperature gives generally very accurate results, except in the case when the expansion method is used in the vicinity of resonances and thresholds.