text:
# Indirect signals from light neutralinos in supersymmetric models without gaugino mass unification

## I Introduction

In Supersymmetric models without gaugino-mass unification at the grand unification scale, neutralinos can be lighter than the current lower bound of 50 GeV, which instead occurs in the case of gaugino–universal models. In Refs. lowneu ; lowmass  we discussed the properties of these light neutralinos as relic particles ($`R`$–parity conservation is assumed) and showed that an absolute lower limit of 6 GeV on the neutralino mass $`m_\chi `$ can be placed by applying the most recent determinations of the upper bound on the Cold Dark Matter (CDM) content in the Universe, in combination with constraints imposed on the Higgs and supersymmetric parameters by measurements at colliders and other precision experiments, like the muon anomalous magnetic moment and the rare decay $`bs+\gamma `$. In Refs. lowneu ; lowmass  we also showed that direct detection rates for light relic neutralinos make these particles detectable with WIMP direct search experiments with current technologies. A comparison of our predictions with intervening experimental results was presented in Ref. lowdir  .

In the present paper we examine light neutralinos in connection with the indirect signals which can be produced by neutralino self–annihilations in the galactic halo or inside celestial bodies. We compare our theoretical predictions with current experimental data on measurements of gamma–rays and antiprotons in space and of up-going muons at neutrino telescopes. Results are presented for a wide range of the neutralino mass, from the established lower bound of 6 GeV up to 500 GeV. However, our discussions are focused on light neutralinos (i.e. neutralinos with $`m_\chi <`$ 50 GeV), since these are not usually considered in current literature otherlow .

The structure of the paper is the following. In Sec. II we briefly summarize the gaugino non–universal supersymmetric model and the properties of light neutralinos which arise in this framework. In Sec. III we discuss the dark matter density distribution in the galactic halo, which is relevant to indirect detection signals, especially to the gamma–ray flux. In Sec. IV we present the calculation and comparison with data of the gamma–ray flux: we consider both the signal coming from the galactic center and from high galactic latitudes. In Sec. V we discuss the antiproton signal, whereas in Sec. VI we show our results for the indirect signals at neutrino telescopes. Our conclusions are drawn in Sec. VII.

## II Supersymmetric models without gaugino–mass unification

A typical assumption of supersymmetric models is the unification condition for the three gaugino masses $`M_{1,2,3}`$ at the GUT scale: $`M_1=M_2=M_3`$. This hypothesis implies that at the electroweak scale $`M_10.5M_2`$. Under this unification condition the bound on the neutralino mass is determined to be $`m_\chi >50`$ GeV. This is derived from the lower bound on the chargino mass (which depends on $`M_2`$ but not on $`M_1`$) determined at LEP2: $`m_{\chi ^\pm }>`$ 100 GeV. By allowing a deviation from gaugino–universality, the neutralino can be lighter than in the gaugino–universal models when $`M_1RM_2`$, with $`R<0.5`$. In this case current data from accelerators do not set an absolute lower bound on $`m_\chi `$.

We consider here an extension of the MSSM which allows for a deviation from gaugino–mass universality by the introduction of the parameter $`R`$, varied here in the interval: $`(0.01÷0.5)`$. This range for $`R`$ implies that the accelerator lower bound on the neutralino mass can be moved down to few GeV for $`R0.01`$. The ensuing light neutralinos have a dominant bino component; a deviation from a pure bino composition is mainly due to a mixture of $`\stackrel{~}{B}`$ with $`\stackrel{~}{H}_1^{}`$ lowneu ; lowmass ; lowdir . Notice that our range of $`R`$ includes also the usual model with gaugino–mass universality.

We therefore employ an effective MSSM scheme at the electroweak scale, defined in terms of a minimal number of parameters, only those necessary to shape the essentials of the theoretical structure of MSSM and of its particle content, supplemented by the gaugino non–universality parameter $`R`$. The assumptions that we impose at the electroweak scale are: a) all squark soft–mass parameters are degenerate: $`m_{\stackrel{~}{q}_i}m_{\stackrel{~}{q}}`$; b) all slepton soft–mass parameters are degenerate: $`m_{\stackrel{~}{l}_i}m_{\stackrel{~}{l}}`$; c) all trilinear parameters vanish except those of the third family, which are defined in terms of a common dimensionless parameter $`A`$: $`A_{\stackrel{~}{b}}=A_{\stackrel{~}{t}}Am_{\stackrel{~}{q}}`$ and $`A_{\stackrel{~}{\tau }}Am_{\stackrel{~}{l}}`$. As a consequence, the supersymmetric parameter space consists of the following independent parameters: $`M_2,\mu ,\mathrm{tan}\beta ,m_A,m_{\stackrel{~}{q}},m_{\stackrel{~}{l}},A`$ and $`R`$. In the previous list of parameters we have denoted by $`\mu `$ the Higgs mixing mass parameter, by $`\mathrm{tan}\beta `$ the ratio of the two Higgs v.e.v.’s and by $`m_A`$ the mass of the CP-odd neutral Higgs boson.

In the numerical random scanning of the supersymmetric parameter space we have used the following ranges: $`1\mathrm{tan}\beta 50`$, $`100\mathrm{GeV}|\mu |`$, $`M_21000\mathrm{GeV}`$, $`100\mathrm{GeV}m_{\stackrel{~}{q}},m_{\stackrel{~}{l}}1000\mathrm{GeV}`$, $`\mathrm{sign}(\mu )=1,1`$, $`90\mathrm{GeV}m_A1000\mathrm{GeV}`$, $`3A3`$, in addition to the above mentioned range $`0.01R0.5`$. We impose the experimental constraints: accelerators data on supersymmetric and Higgs boson searches and on the invisible width of the Z boson, measurements of the branching ratio of the $`bs+\gamma `$ decay and of the upper bound on the branching ratio of $`B_s\mu ^++\mu ^{}`$, and measurements of the muon anomalous magnetic moment $`a_\mu (g_\mu 2)/2`$. The range used here for the $`bs+\gamma `$ branching ratio is $`2.1810^4BR(bs+\gamma )4.2810^4`$ bsgamma . For the branching ratio of $`B_s\mu ^++\mu ^{}`$ we employ the upper limit $`BR(B_s\mu ^++\mu ^{})<7.510^7`$ (95% C.L.) bsmumu\_exp ; for the theoretical evaluation we have used the results of bsmumu\_th  with inclusion of the QCD radiative corrections to the $`b`$ Yukawa coupling carena . For the deviation of the current experimental world average of $`a_\mu `$ from the theoretical evaluation within the Standard Model we use the 2$`\sigma `$ range: $`142\mathrm{\Delta }a_\mu 10^{11}474`$; this interval takes into account the recent evaluations of Refs. davier ; hagiwara2 . We notice that gluinos do not enter directly into our loop contributions to $`BR(bs+\gamma )`$ and $`BR(B_s\mu ^++\mu ^{})`$, since we assume flavor-diagonal sfermion mass matrices. Thus, gluinos appear only in the QCD radiative corrections to the $`b`$ Yukawa coupling; in the evaluation of these effects the value of the relevant mass parameter $`M_3`$ is taken at the standard unification value $`M_3=M_2\alpha _3(M_Z)/\alpha _2(M_Z)`$, where $`\alpha _3(M_Z)`$ and $`\alpha _2(M_Z)`$ are the SU(3) and SU(2) coupling constants evaluated at the scale $`M_Z`$.

The new data on the cosmic microwave background from WMAPwmap , used in combination with other cosmological observations, mainly galaxy surveys and Lyman–$`\alpha `$ forest data, are sharpening our knowledge of the cosmological parameters, and in particular of the amount of dark matter in the Universe. From the analysis of Ref. wmap , we obtain a restricted range for the relic density of a cold species like the neutralinos. The density parameter of cold dark matter is bounded at $`2\sigma `$ level by the values: $`(\mathrm{\Omega }_{CDM}h^2)_{\mathrm{min}}=0.095`$ and $`(\mathrm{\Omega }_{CDM}h^2)_{\mathrm{max}}=0.131`$. This is the range for CDM that we consider in the present paper. An independent determination for the content of cold dark matter in the Universe is provided by the Sloan Digital Sky Survey Collaboration sloan ; this new data agrees with the results of Ref. wmap .

We recall that the relic abundance $`\mathrm{\Omega }_\chi h^2`$ is essentially given by $`\mathrm{\Omega }_\chi h^2<\sigma _{\mathrm{ann}}v>_{\mathrm{int}}^1`$, where $`<\sigma _{\mathrm{ann}}v>_{\mathrm{int}}`$ is the thermally–averaged neutralino annihilation cross section times average velocity, integrated from the freeze–out temperature $`T_f`$ to the present one $`T_0`$. The quantity $`\sigma _{\mathrm{ann}}`$ enters also in the calculation of the indirect signals that will be discussed in the following sections. In the evaluation of $`\sigma _{\mathrm{ann}}`$ we have considered the full set of available final states: fermion-antifermion pairs, gluon pairs, pairs of charged Higgs bosons, one Higgs boson and one gauge boson, pairs of gauge bosons anni . We have not included coannihilation in our evaluation of the neutralino relic abundance, since, at variance from a constrained SUGRA scheme, in our effective supersymmetric model a matching of the neutralino mass with other masses is accidental, i.e. not induced by some intrinsic relationship among different parameters. Introducing coannihilation would only produce an insignificant reshuffle in the representative points of the scatter plots displayed in the present paper, without a modification of their borders, which are the only feature of physical significance.

The relic abundance $`\mathrm{\Omega }_\chi h^2`$ of neutralinos lighter than 50 GeV which arise in our class of gaugino non–universal models has a relatively simple structure in terms of dominant diagrams in the annihilation cross section lowneu ; lowmass . Here we just remind that combining our calculation of the relic abundance of light neutralinos with the value of $`(\mathrm{\Omega }_{CDM}h^2)_{\mathrm{max}}`$, an absolute lower bound on the neutralino mass of 6 GeV can be set lowneu ; lowmass . We note that within our present scanning of the supersymmetric parameter space, the lower limit on the neutralino mass shifts to about 7 GeV, when the upper bound on $`BR(B_s\mu ^++\mu ^{})<7.510^7`$ is implemented (this constraint was not included in lowneu ; lowmass ). It is remarkable that a lower limit on $`m_\chi `$ is set not by searches at accelerators, but instead by cosmological arguments.

## III Dark matter in the Galaxy

Signals due to neutralino self–annihilation in the halo depend quadratically on the dark matter density distribution $`\rho (\stackrel{}{r})`$, and are therefore very sensitive to the features of this physical quantity. Two properties are of special relevance: 1) the behavior of the density distribution in the galactic center (GC), 2) the extent of the density contrast (clumpiness), which represents the deviation of the actual density distribution from a smooth distribution.

The most commonly used density distributions can be parametrized by the following spherically–averaged density profile hernquist :

$$\rho (r)=\rho _l\left(\frac{r_{}}{r}\right)^\gamma \left[\frac{1+(r_{}/a)^\alpha }{1+(r/a)^\alpha }\right]^{(\beta \gamma )/\alpha },$$
(1)
where $`r=|\stackrel{}{r}|`$, $`r_{}=8`$ kpc S2  is the distance of the Sun from the galactic center, $`a`$ is a scale length and $`\rho _l`$ is the total local (solar neighborhood) dark matter density. In particular, the isothermal density profile corresponds to $`(\alpha ,\beta ,\gamma )`$ = (2,2,0), the Navarro, Frenk and White (NFW) profile nfw  corresponds to $`(\alpha ,\beta ,\gamma )`$ = (1,3,1) and the Moore et al. profile moore  to $`(\alpha ,\beta ,\gamma )`$ = (1.5,3,1.5). The two latter profiles, both derived from numerical simulations of structure formation, differ noticeably in their behavior at small distances from the GC: $`r^1`$ for the NFW and $`r^{1.5}`$ for the Moore et al. profile, with ensuing large differences in the size of the expected signals for WIMP annihilation from the central region of the Galaxy.

Recent results of extensive numerical simulations, aimed at an analysis of the inner structure of $`\mathrm{\Lambda }CDM`$ halos, strongly disfavor a behavior as singular as $`r^{1.5}`$, but also indicate that a NFW profile is likely not to be adequate at small distances from the GC navarro . It turns out that the density profile is not described by a singular power–law at small distances; rather, in this asymptotic regime, the numerical results are well fitted by a profile whose logarithmic slope $`\delta (r)d(\mathrm{ln}\rho (r))/d(\mathrm{ln}r)`$ is given by $`\delta (r)r^\alpha `$ with $`\alpha 0.17`$. This leads to a non–singular dark matter density distribution function of the form navarro :

$$\rho (r)=\rho _2\mathrm{exp}\left\{\frac{2}{\alpha }\left[\left(\frac{r}{r_2}\right)^\alpha 1\right]\right\},$$
(2)
where $`r_2`$ is the radius where the logarithmic slope is $`\delta =2`$, and $`\rho _2\rho (r_2)`$. These various distributions mainly differ in their behavior at small $`r`$. We wish to stress here that anyway current cosmological simulations are not reliable for radii smaller than an $`r_{\mathrm{min}}0.1`$ – 1 kpc. We also notice that singular profiles are subject of debate in current literature, with analyses pointing to inconsistencies with observational data on rotational curves salucci .

Furthermore we recall that other density profiles are able to describe the dark matter halo, including for instance different classes of logarithmic and power–law potentials, axisymmetric distributions or even triaxial distributions bcfs . In the following we will concentrate on the standard isotropic density profiles of Eq. (1) and on the new profile of Eq. (2) deduced from numerical simulations. For definiteness, we use as a reference model the NFW density distribution, and discuss the deviation from this reference case when the other density profiles are considered. Most signals do not depend or are only mildly dependent on the critical behavior of the dark matter density around the GC: this occurs for the neutrino fluxes from the Earth and the Sun, for antiprotons and for gamma–rays coming from large galactic latitudes. On the contrary, gamma-rays coming from the GC are very sensitive to the inner parts of the Galaxy and the differences between the different halo profiles will be explicitly discussed.

A key parameter for all the density distributions is the value for the total local dark matter density $`\rho _l`$. This parameter can be determined for each density profile assuming compatibility with the measurements of rotational curves and the total mass of the Galaxy bcfs . For instance, a simple modelling of the visible and dark components of the Galaxy showed that $`\rho _l`$ can range from 0.18 GeV cm-3 to 0.71 GeV cm-3 for an isothermal sphere profile, from 0.20 GeV cm-3 to 1.11 GeV cm-3 for a NFW distribution and 0.22 GeV cm-3 to 0.98 GeV cm-3 for a Moore et al. shape bcfs . For definiteness, our results will be presented for $`\rho _l=0.3`$ GeV cm-3 for all the density profiles employed in the present analysis. The parameter $`\rho _l`$ enters as a mere scaling factor in the signal fluxes: the effect of varying $`\rho _l`$ is therefore easily taken into account.

Once the density profile which describes the total dark matter density in the galactic halo is chosen, the actual neutralino density distribution is taken to be:

$$\rho _\chi (r)=\xi \rho (r),$$
(3)
where $`\xi `$ accounts for the fact that neutralinos could be only a fraction of the total cold dark matter ($`\xi 1`$). This characteristic is linked to the actual relic abundance of neutralinos and is accounted for by using the standard rescaling prescription: $`\xi =\mathrm{min}[1,\mathrm{\Omega }_\chi h^2/(\mathrm{\Omega }_{CDM}h^2)_{\mathrm{min}}]`$.

As was noticed in Refs. bss ; ss ; bbm  an effect of density contrast in the dark matter distribution could produce a strong enhancement effect in signals due to $`\chi `$-$`\chi `$ annihilations in the halo. This property was subsequently considered in connection with various signals (gamma rays, positrons, antiprotons) clump ; clump2 , sometimes under the assumption of a strong clumpiness effect, at the level of a few orders of magnitude. However, according to a recent analytical investigation on the production of small-scale dark matter clumps bde , the clumpiness effect would not be large, with the result that the ensuing enhancement on the annihilation signals is limited to a factor of a few. Similar conclusions are also reached in Ref. setal  by using results of high-resolution numerical simulations.

## IV Gamma rays

The flux of gamma rays $`\mathrm{\Phi }_\gamma (E_\gamma ,\psi )`$ originated from neutralino pair annihilation in the galactic halo bss ; clump ; gamma\_history  and coming from the angular direction $`\psi `$ is given by:

$$\mathrm{\Phi }_\gamma (E_\gamma ,\psi )=\frac{1}{4\pi }\frac{\sigma _{\mathrm{ann}}v}{m_\chi ^2}\frac{dN_\gamma }{dE_\gamma }\frac{1}{2}I(\psi ),$$
(4)
where $`\sigma _{\mathrm{ann}}v`$ is the annihilaton cross section times the relative velocity mediated over the galactic velocity distribution function and $`dN_\gamma /dE_\gamma `$ is the energy spectrum of $`\gamma `$-rays originated from a single neutralino pair annihilation. The quantity $`I(\psi )`$ is the integral of the squared dark matter density distribution performed along the line of sight:

$$I(\psi )=_{\mathrm{l}.\mathrm{o}.\mathrm{s}}\rho ^2(r(\lambda ,\psi ))𝑑\lambda (\psi ),$$
(5)
and $`\psi `$ is the angle between the line of sight and the line pointing toward the GC. The angle $`\psi `$ is related to the galactic longitude $`l`$ and latitude $`b`$ by the expression $`\mathrm{cos}\psi =\mathrm{cos}l\mathrm{cos}b`$. A point at a distance $`\lambda `$ from us and observed under an angle $`\psi `$ is therefore located at the galactocentric distance $`r=\sqrt{\lambda ^2+r_{}^22\lambda r_{}\mathrm{cos}\psi }`$. The factor of $`1/2`$ in Eq.(4) is due to the fact that the gamma–ray flux depends on the number of neutralino pairs present in the galactic halo, as pointed out in Ref. clump1 . This factor of $`1/2`$ applies as well to any other indirect detection signal which depends on the annihilation of a pair of Majorana fermions in the galactic halo, like positrons, antiprotons, antideuterium. In the case of dark matter composed of Dirac fermions, the statistical factor would instead be $`1/4`$, if $`\rho (r)`$ describes the total dark matter distribution ascribed to the given Dirac species.

### IV.1 The source spectrum

As far as the annihilation of light neutralinos is concerned (namely for neutralino masses below the thresholds for gauge–bosons, higgs-bosons and $`t`$ quark production), the production of $`\gamma `$–rays in the continuum takes contribution mainly from the hadronization of quarks and gluon pairs produced in the neutralino annihilation process. The subsequent $`\pi ^0`$ production and decay $`\pi ^02\gamma `$ give usually the dominant contribution. In this case the $`\gamma `$–ray energy spectrum is given by:

$$\frac{dN_\gamma }{dE_\gamma }=_{E_\pi ^{\mathrm{min}}}^{E_\pi ^{\mathrm{max}}}P(E_\pi ,E_\gamma )\frac{dN_\pi }{dE_\pi }𝑑E_\pi ,$$
(6)
where $`P(E_\pi ,E_\gamma )=2(E_\pi ^2m_\pi ^2)^{1/2}`$ is the probability per unit energy to produce a $`\gamma `$–ray with energy $`E_\gamma `$ out of a pion with energy $`E_\pi `$, while $`dN_\pi /dE_\pi `$ is the pion yield per annihilation event.

We have evaluated the quantity $`dN_\pi /dE_\pi `$ by means of a Monte Carlo simulation with the PYTHIA package pythia . The Monte Carlo code has been run by injecting $`q\overline{q}`$ and gluon pairs back–to–back at fixed center–of–mass energy $`E_{\mathrm{cm}}=2m_\chi `$. Since quarks and gluons are confined, they contribute to a complex final–state pattern of out–going hadronic strings decaying to physical hadrons through fragmentation. In the Lund string scheme, fragmentation is an intrinsically scale invariant process. This implies that, in the rest frame of the decaying string, the final–state spectrum is invariant in the variable $`xE_F/m_{\mathrm{string}}`$, where $`E_F`$ is the energy of the given final state and $`m_{\mathrm{string}}`$ is the total string mass. Were it not for showered gluons, a $`q\overline{q}`$ pair from neutralino annihilation would produce, in the reference frame of the two annihilating neutralinos, a single hadronic string at rest with $`m_{\mathrm{string}}=2m_\chi `$, subsequently fragmenting to produce (among other particles) a pion spectrum which would be scale–invariant in the variable $`yE_\pi /m_\chi `$. However, due to showering, this scale invariance is significantly broken, since the pion energy–spectrum is given by the superposition of different decaying strings boosted at different energies. Therefore the pion spectrum at a given neutralino mass cannot be obtained from that calculated at a different $`m_\chi `$.

We have therefore evaluated the pion yield per annihilation event, $`dN_\pi ^F(m_\chi ,E_\pi )/dE_\pi `$, for each final state $`F=f\overline{f},gg`$, at different neutralino masses: $`m_\chi =`$ 6, 10, 50, 100, 500 and 1000 GeV and for pion energies ranging from $`m_{\pi ^0}`$ to $`m_\chi `$ in 100 equal bins in logarithmic scale. In order to optimize our numerical calculations, we then obtain the pion spectrum at neutralino masses and pion energies different from the ones sampled through a two dimensional numerical interpolation. We have explicitly checked with the PYTHIA Monte Carlo results that our interpolation is accurate at the percent level both on the reconstructed pion yield and on the final gamma–ray spectrum from $`\pi ^0`$ decay, as given by Eq. (6).

The contribution to the $`\gamma `$–ray spectrum from production and decay of mesons other than pions (mostly $`\eta `$, $`\eta ^{}`$, charmed and bottom mesons) and of baryons is usually subdominant as compared to $`\pi ^0`$ decay. These additional contributions can be safely neglected (they typically contribute only up to 10% of the total flux for $`E_\gamma <`$ 1 GeV). A notable exception is given by the hadronization of $`b\overline{b}`$ pairs at low production energies, i.e. for neutralino masses between the production threshold for a $`b`$–meson and about 10 GeV. In this case jet flavor conservation leads to the production of a bottom meson $`B=B^0,B^\pm ,B_s^0`$ with 100 % probability. In the PYTHIA code, 75% of the times the $`B`$ meson is in the excited state and decays through $`B^{}B+\gamma `$ with $`m_B^{}m_B`$ 46 MeV. Since the $`B^{}`$ mesons are produced almost at rest ($`m_B5.3`$ GeV) for $`m_\chi <10`$ GeV, they generate a (slightly boosted) gamma–ray line that dominates the other contributions below $`E_\gamma `$ 100 MeV. We have then included this peculiar contribution to our interpolating procedure note .

Neutralino annihilation into lepton pairs can also produce gamma–rays from electromagnetic showering of the final state leptons. This process can be dominant for $`E_\gamma <`$ 100 MeV, when the neutralino annihilation process has a sizable branching ratio into lepton pairs. In the case of production of $`\tau `$ leptons, their semihadronic decays also produce neutral pions, which then further contribute to the gamma–ray flux. Also these additional contributions are included in our numerical evaluations, again by modelling the gamma–ray production with the PYTHIA Monte Carlo for the same set of neutralino masses quoted above, and by numerically interpolating for other values of $`m_\chi `$.

When the neutralino masses are sufficiently large, the annihilation channels into Higgs bosons, gauge bosons and $`t\overline{t}`$ pairs become kinematically accessible. We compute analytically the full decay chain down to the production of a quark, gluon or a lepton. The ensuing $`\gamma `$–ray spectrum is then obtained by using the results discussed above for quarks, gluons and leptons, by properly boosting the differential energy distribution to the rest frame of the annihilating neutralinos (see e.g. Appendix I in Ref. pbar\_susy  for details).

### IV.2 The geometrical factor

The integral along the line of sight $`I(\psi )`$ in Eq. (5) is the quantity that takes into account the shape of the dark matter profile. For small values of $`\psi `$, $`I(\psi )`$ is very sensitive to possible enhancement of the density at the GC and therefore large differences for the different density profiles are expected.

When comparing to experimental data, Eq.(5) is averaged over the telescope aperture–angle $`\mathrm{\Delta }\psi `$:

$$I_{\mathrm{\Delta }\psi }=\frac{1}{\mathrm{\Delta }\psi }_{\mathrm{\Delta }\psi }I(\psi )𝑑\psi .$$
(7)
The gamma–ray flux is therefore proportional to $`I_{\mathrm{\Delta }\psi }`$.

In our analysis on the galactic center emission, we will employ data from the EGRET experiment egret\_gc ; egret\_gc\_source , whose angular resolution is given by the longitude–latitude aperture: $`|\mathrm{\Delta }l|5^{}`$, $`|\mathrm{\Delta }b|2^{}`$. Table 1 shows the values of $`I_{\mathrm{\Delta }\psi }`$ for the density profiles discussed in Sec. III. The effect of changing the core radius of an isothermal sphere is not negligible: an increase of a factor of 2.3 is obtained by reducing $`a`$ from 3.5 kpc to 2.5 kpc. In the case of singular distributions a small–distance cut–off of $`r_c=0.01`$ pc is assumed (inside $`r_c`$ the density is assumed to be constant). A NFW profile then gives a flux which is about 10 times larger than an isothermal sphere with $`a=3.5`$ kpc. The very steep Moore et al. profile would produce a flux about 60 times larger than a NFW. The Navarro et al. navarro  profile with $`r`$–dependent log–slope sits between the NFW and Moore et al. cases: though not singular, it nevertheless provides a flux which is about 3 times larger than a singular NFW halo. The variability of $`I_{\mathrm{\Delta }\psi }`$ on the dark matter profile can therefore be as large as a factor of 600, comparing the very steep Moore et al. profile with an isothermal sphere with a large core radius. However, the recent critical analysis on numerical simulation of Ref. navarro  implies that a factor of 30 is likely to be a more plausible interval.

In the case of high galactic latitudes, the dependence of $`I_{\mathrm{\Delta }\psi }`$ on the density profile is much milder. At these high latitudes, EGRET identifies a residual gamma–ray flux which is ascribed to a possible extragalactic component, but could as well be due to dark matter annihilation in the galactic halo. We will consider in our analysis two different angular regions: $`|b|>10^{}`$ with the exclusion of $`|l|40^{}`$ and $`10^{}|b|30^{}`$ around the galactic center egret\_extragal  (region A) and $`|b|86^{}`$ (region B). Region B has been considered in the reanalysis of EGRET data by the authors of Ref. waxman , where more stringent limits on the the gamma-ray residual intensity from the galactic poles have been derived. The value of $`I_{\mathrm{\Delta }\psi }`$ for region A is 1.66 for a NFW profile, and ranges from 1.61 for the isothermal sphere with $`a=3.5`$ kpc to 1.80 for the $`r`$–dependent log–slope profile: when pointing away from the critical behavior of the density profiles at galactic center, the line of sight integral is almost universal. In the case of region B, $`I_{\mathrm{\Delta }\psi }=0.67`$ for the NFW distribution, and ranges from 0.62 for the isothermal sphere with $`a=3.5`$ kpc to 0.69 for the $`r`$–dependent log–slope profile. In both cases, the Moore et al. profile gives a line–of–sight integral slightly smaller than in the case of the $`r`$–dependent log–slope profile.

### IV.3 Signal from the galactic center

Data from low galactic latitudes ($`|b|<10^{}`$), including the galactic center region, have been collected by the EGRET telescope egret\_gc . The diffuse gamma–ray flux of the inner galaxy measured by EGRET shows a possible excess over the estimated background at energies larger than about 1 GeV.

Clearly a firm assessment of an excess requires a good knowledge of the standard production of gamma rays in our Galaxy. At the energies of interest for our analysis - namely, from about 100 MeV to tens of GeV - the main production mechanism of $`\gamma `$-rays is the interaction of cosmic rays (mainly protons and helium nuclei) with the interstellar medium (atomic and molecular hydrogen, and helium). In these strong reactions $`\pi ^0`$’s are produced, and hence $`\gamma `$-rays via pion decay: $`\pi ^02\gamma `$. The ensuing spectrum has a bump around 70 MeV, and drops at high energies with an energy power law which follows the progenitor cosmic ray spectrum ($`E^\alpha `$, with $`\alpha `$ 2.7). Another source of $`\gamma `$-rays comes from inverse-Compton scattering of cosmic ray electrons off the interstellar photons. In particular, energetic electrons may scatter off the cosmic microwave background, and off the infrared, optical and ultraviolet radiation arising from stellar activity and dust. The third radiation component originates from electron bremsstrahlung over the interstellar medium, which may be partially or even totally ionized. Bremsstrahlung $`\gamma `$-rays are mostly important in the low energy tail. For a full calculation of these three main radiation components one needs a good knowledge of the physics of cosmic rays and of the interstellar medium in the region of interest. This is particularly unlikely when dealing with the galactic center area.

In the literature, several different results have been achieved on the subject. First of all, the EGRET Collaboration developed a detailed calculation of the $`\gamma `$–ray background at the energies of interest for the detector egret\_gc ; egret\_calc : this calculation shows a clear deficit of $`\gamma `$-rays towards the GC with respect to measurements. The excess in the data is apparent for $`E_\gamma >1`$ GeV, where the shapes of the spectra of the estimated background and the data differ significantly. At lower energies the spectral agreement is instead rather good. Similar conclusions have been drawn in Refs. mori  with some different procedure in the calculation of the background. In this paper a harder, probably unrealistic, nucleonic spectrum is shown not to be anyway sufficient to explain the GeV excess. Some modifications towards harder electron and nucleon spectra are studied in Ref. strong , but a satisfactory agreement with data is not achieved (notice that it is anyway hard to reconcile these hypothesis with galactic cosmic rays measurements). The results of all these analyzes favor the interpretation of the EGRET data in terms of an excess over the background, mostly at energies above 1 GeV.

However, different assumptions on acceleration busching  and diffusion erlykin  of cosmic ray nucleons, and on the spectral shape of primary nucleons in the interstellar space aharonian , have been proposed and lead to a quite good agreement of the EGRET measured flux: almost all the spectral features are reproduced by these calculations. In this case, the EGRET data would be explained in terms of the standard galactic $`\gamma `$–ray production.

At present, it is very difficult to favor one model against the others, both on theoretical and observational basis. This implies that the uncertainty on the calculation of the galactic $`\gamma `$-ray flux, and in particular at the galactic center, is very difficult to quantify.

Due to these open problems on the determination of the background component, we will develop our analysis along two paths. First of all we will discuss whether, and under which conditions, it is possible to set constraints on low–mass neutralinos from the gamma-ray studies. Then we will comment on the possibility for low–mass neutralinos to explain the putative EGRET excess.

The gamma–ray flux from the galactic center inside the angular region $`|\mathrm{\Delta }l|5^{}`$, $`|\mathrm{\Delta }b|2^{}`$ for a NFW matter density profile is shown in Fig. 1 at three representative photon energies: $`E_\gamma =`$ 0.12, 1.5 and 15 GeV. These energies correspond to three energy bins of the EGRET detector. The scatter plots of the two top panels display a peculiar funnel shape at small masses. This is due to the fact that the neutralino flux is bounded from below by the cosmological limit $`\mathrm{\Omega }_\chi h^2(\mathrm{\Omega }_{CDM}h^2)_{\mathrm{max}}`$. This feature is similar to the one we found in Refs. lowneu ; lowmass ; lowdir , in connection with neutralino direct detection rates. The variation in shape of the scatter plots, when $`E_\gamma `$ is increased, is easily understood in terms of Eq.(4). At $`E_\gamma =0.12`$ GeV the $`m_\chi ^2`$ behavior is clearly visible. Energies of the order of 100 MeV are very crucial in offering the possibility to set limits on the very light neutralino sector. By increasing the photon energy, the lightest neutralinos do not have enough phase space to produce photons at this energy (since they annihilate almost at rest in the Galaxy): therefore the gamma–ray flux at very low masses becomes progressively depressed, as $`E_\gamma `$ increases. At $`E_\gamma =1.5`$ GeV a neutralino with a mass of 6 GeV can produce approximately the same flux as a 10-15 GeV neutralino. At $`E_\gamma =15`$ GeV all the neutralinos lighter than 15 GeV obviously do not produce any photon: at this energy the maximal fluxes are obtained for neutralinos with masses around 30 GeV.

The first conclusion which can be drawn from Fig. 1 is that the supersymmetric model considered here is not constrained, at present, by EGRET data for a NFW density profile. An increase in the flux by a factor of 3.3, as would be in the case of the $`r`$–dependent log–slope profile of Eq.(2), is also not enough to set limits.

Larger enhancements in the geometrical factor $`I_{\mathrm{\Delta }\psi }`$ over the NFW case are necessary to set limits. The comparison of the three panels of Fig. 1 shows that the limits come from different energies $`E_\gamma `$, depending on the neutralino mass. For very light neutralinos, i.e. $`m_\chi <10`$ GeV, the lowest energy bin is the relevant one. In this case an enhancement of a factor of 6 would allow to raise the predicted fluxes for $`E_\gamma =0.12`$ GeV at the level of the EGRET measurement, and therefore to start setting limits. For masses in the range $`10\mathrm{GeV}<m_\chi <20\mathrm{GeV}`$ the $`E_\gamma =1.5`$ GeV bin sets more stringent limits, at least on a fraction of the supersymmetric models, but only for a factor of enhancement of at least 15–20 over the NFW case. These factors are pretty large, even though not as large as the one which refers to a Moore et al. profile, which is about 60, as discussed before. For masses around 30–40 GeV the best limits come from the highest energy bin $`E_\gamma =15`$ GeV, where a factor of 20–25 would allow the fluxes to reach the EGRET data. In the case of the standard MSSM, where the neutralino has masses larger than 50 GeV the lowest energy bins are always less constraining that the higher energy ones, as can be seen by comparing the different panels of Fig. 1. Instead, the lower energy bins are crucial for the study of the low–mass neutralinos. We finally comment that a Moore et al. profile would make all the fluxes for $`m_\chi <10`$ GeV incompatible with the data, but this profile is less likely, as we discussed above.

Now, let us turn to a brief discussion of the possibility to explain the EGRET excess by means of low–mass neutralinos. The analysis made above on the behavior of the gamma–ray fluxes in the three representative energy bins of Fig. 1 shows that neutralinos in the mass range $`25\mathrm{GeV}<m_\chi <40\mathrm{GeV}`$ are the ones which may have the possibility to fill the excess in the energy range above 1 GeV, without spoiling the lower energy behavior of the background which is supposed to have an acceptable agreement with the data. We show that indeed the low–mass neutralinos in this mass range are able to explain the EGRET excess in Fig. 2. In this figure we plot the predicted gamma–ray spectra for two representative supersymmetric configuration when, for definiteness, the gamma–ray background as calculated in Ref. egret\_calc  is assumed. In the left panel of Fig. 2 we show a supersymmetric configuration with $`m_\chi =30`$ GeV and a relic abundance in the proper range to explain dark matter: $`\mathrm{\Omega }_\chi h^2=0.12`$. The dominant annihilation channels of these low–mass neutralinos are $`\chi \chi \overline{\tau }\tau `$ and $`\chi \chi \overline{b}b`$ lowneu ; lowmass . Gamma rays coming from annihilation into $`\tau `$’s give a harder spectrum as compared to the $`b`$ channel. In this representative point the two channels have (approximately) the same branching ratio: this gives sizable contribution to the gamma–ray flux in all the energy range from 1 to 10 GeV, which is where the excess in the EGRET data is more pronounced. By allowing the background to fluctuate down by 10% and by using a geometrical factor $`I_{\mathrm{\Delta }\psi }`$ 30 times larger than in the NFW case, we show that a pretty good agreement between the total flux and the data can be obtained. We are not quoting a statistical significance for this agreement since we are not performing a systematic statistical analysis here: however, we are interested in showing that, in addition to the standard MSSM neutralinos with masses larger than 50 GeV ullio\_morselli , also low–mass neutralinos have the capability of explaining the putative EGRET excess. In both cases, low–mass and standard neutralinos, the values of the line–of–sight integrals $`I_{\mathrm{\Delta }\psi }`$ which are able to explain the EGRET excess are much larger than what is provided by a NFW density profile. However, for neutralinos in the 30–40 GeV mass range these enhancement factors are smaller than in the case of heavier neutralinos, due to the $`m_\chi ^2`$ behavior.

The right panel of Fig. 2 shows a second example, with a neutralino of $`m_\chi =40`$ GeV and $`\mathrm{\Omega }_\chi h^2=0.11`$. The background of Ref. egret\_calc  is again scaled down by 10%. In this case the geometrical factor is 32 times the NFW one. The branching fraction of annihilation into $`b`$ quarks is larger than in the previous case. This enhances the contribution to the gamma–ray flux in the 1–3 GeV range without spoiling the agreement at larger energies.

A detailed analysis of the spectral features of the gamma–ray fluxes produced by neutralino annihilation and their comparison with the EGRET data is beyond the scope of the present paper and will be presented elsewhere.

### IV.4 Signal from high galactic latitude

Data from high galactic latitudes have been collected by the EGRET telescope egret\_extragal . An analysis of the measurements taken over the latitudes $`|b|>10^o`$, and excluding the region $`|l|<40^o`$ and $`10^o<|b|<30^o`$ around the galactic center, has been performed in Ref. egret\_extragal . All the identified sources as well as the components due to the interactions of cosmic rays with the galactic disk gas have been subtracted egret\_extragal . The residual flux, averaged over the considered portion of the sky, has been shown to be isotropic and well fitted by the power law $`\mathrm{\Phi }_{HL}^{\mathrm{EGRET}}(E_\gamma )=k(E_\gamma /E_0)^\alpha `$, where $`k=(7.32\pm 0.34)10^6`$ photons cm-2 sec-1 sr-1 GeV-1, $`\alpha =2.10\pm 0.03`$ and $`E_0=451`$ MeV. This spectrum is often referred to as the extragalactic diffuse emission, since no known source inside the Galaxy seems to be responsible for it. One possibility is that it is due to unresolved gamma-ray–emitting blazars. Relying on the analysis by EGRET, one can use the residual flux as an upper bound to any flux due to exotic sources, including annihilation of relic neutralinos. Recently, a re-analysis of EGRET data has been performed in Ref. waxman , taking particular care of the spatial dependence of the observed photons. Working on the integrated flux and taking into account contributions from several galactic tracers, the authors of Ref. waxman  show that the high-latitude $`\gamma `$-ray sky exhibits strong galactic features and could be well accommodated by simple galactic models. Conservative constraints have been set on the flux integrated above 100 MeV and averaged over different sectors of the sky far from the galactic plane. In this scenario, the room left to an unexplained diffuse flux – often considered as an extragalactic background, but also possibly due to exotic galactic sources – is much smaller (by a factor of three, at least) than the one reported in Ref. egret\_extragal . Here we consider the upper limit of Ref. waxman  on this possible residual, isotropic flux $`I_\gamma `$ in the polar region ($`|b|>86^0`$): $`I_\gamma <0.610^5`$ $`\gamma \mathrm{sec}^1\mathrm{cm}^2\mathrm{sr}^1`$, and compare it to our estimates for the $`\gamma `$-ray flux due to neutralino annihilation averaged in the same spatial region.

The results for both estimates are shown in Fig. 3 for a NFW profile. As in the case of the galactic center emission, also for high latitudes we do not have constraints on the supersymmetric parameter space. Contrary to the case of the galactic center region, for high latitudes the geometrical factor $`I_{\mathrm{\Delta }\psi }`$ is practically independent of the halo profile, as discussed before. We therefore conclude that at present the $`\gamma `$–ray signals from high galactic latitudes do not provide any constraint on the supersymmetric parameter space. The situation could change if further studies will show that a much bigger fraction of the EGRET measured flux at high latitudes is due to galactic foreground or when next–generation experiments will provide further information.

## V Antiprotons in cosmic rays

Like in the case of the gamma–ray flux, the production of antiprotons from neutralino annihilation results from the hadronization of quarks and gluons created in the annihilation process pbar\_history ; pbar0 ; Bottino\_Salati ; bergstrom ; pbar\_susy . The differential rate per single annihilation, unit volume and time is defined as:

$$q_{\overline{p}}^{\mathrm{susy}}(r,T_{\overline{p}})=\sigma _{\mathrm{ann}}v\frac{dN_{\overline{p}}}{dT_{\overline{p}}}\frac{1}{2}\left(\frac{\rho _\chi (r)}{m_\chi }\right)^2,$$
(8)
where $`T_{\overline{p}}`$ denotes the antiproton kinetic energy, $`dN_{\overline{p}}/dT_{\overline{p}}`$ (indicated as $`g(T_{\overline{p}})`$ in Refs.Bottino\_Salati ; pbar\_susy ) is the differential antiproton spectrum per annihilation event, and the factor 1/2 accounts for the number of annihilating neutralino pairs. The spectrum $`dN_{\overline{p}}/dT_{\overline{p}}`$ is evaluated by means of the Monte Carlo simulations we already used in Sec. IV.1.

Once antiprotons are produced in the dark halo, they diffuse and propagate throughout the Galaxy. To describe these processes, we follow the treatment of Ref. pbar\_susy , to which we refer for details. Here we only recall that the propagation of antiprotons has been considered in a two-zone diffusion model PaperI ; PaperI\_bis ; revue , defined in terms of six free parameters whose role is to take into account the main physical processes which affect the propagation of cosmic rays in the Galaxy: acceleration of primary nuclei, diffusion, convective wind, reacceleration process and interaction with the interstellar medium. These free parameters are constrained by comparing the fluxes of various cosmic ray species calculated in our diffusion model with observations. In this regard, the most important observable is the measured Boron/Carbon ratio ($`B/C`$), whose analysis within our diffusion model is presented in Ref. PaperI\_bis . The parameters constrained by the $`B/C`$ measurements have been shown to be compatible with a series of other observed species PaperI\_bis ; pbar\_sec ; beta\_rad , further supporting the employed model. Therefore, in the calculation of the primary antiproton flux we only use those values for the propagation parameters which provide a good statistical agreement to the $`B/C`$ data. One of the main results obtained in Ref. pbar\_susy  is that the supersymmetric antiproton flux, when calculated with the selected propagation parameters, is affected by a large uncertainty. At low energy this uncertainty reaches almost two orders of magnitude, while it diminishes to a factor of thirty at higher energies. Only better and more complete data on cosmic rays (both stable and radioactive) could help in reducing this uncertainty.

At variance with the gamma ray signal, the supersymmetric antiproton flux measured at Earth is almost insensitive to the specific form of the dark matter distribution function, among those discussed in Sec. III. Indeed, these various distributions differ mainly at the galactic center, while at solar neighborhood differences are very mild. Since charged particles, such as antiprotons, suffer enormous energy redistributions, gains and losses, it is almost impossible for an antiproton produced around the galactic center to reach the Earth. This property was shown in Ref. origin , and quantified in Ref. pbar\_susy  for the case of a NFW distribution and an isothermal one.

In the present work, we use directly the results reached in Ref. pbar\_susy , where the function

$$C_{\mathrm{susy}}^{\mathrm{prop}}(T_{\overline{p}})=\frac{\mathrm{\Phi }_{\overline{p}}(r_{},T_{\overline{p}})}{\mathrm{{\rm Y}}g(T_{\overline{p}})}$$
(9)
was calculated. In this equation, $`\mathrm{\Phi }_{\overline{p}}(r_{},T_{\overline{p}})`$ is the interstellar antiproton flux after propagation and $`\mathrm{{\rm Y}}`$ is the supersymmetric flux factor:

$$\mathrm{{\rm Y}}=\frac{1}{2}\xi ^2\frac{\sigma _{\mathrm{ann}}v}{m_\chi ^2}.$$
(10)
The quantity $`C_{\mathrm{susy}}^{\mathrm{prop}}(T_{\overline{p}})`$ may be considered as a measure of how the source flux $`q_{\overline{p}}^{\mathrm{susy}}(r,T_{\overline{p}})`$ is modified by the propagation of antiprotons in the Galaxy before reaching the heliosphere. In the results presented in the following, we have calculated the antiproton flux $`\mathrm{\Phi }_{\overline{p}}(r_{},T_{\overline{p}})`$ according to Eq. (9), where the $`C_{\mathrm{susy}}^{\mathrm{prop}}(T_{\overline{p}})`$ function has been taken directly from Ref. pbar\_susy  for a few representative combinations of the propagation parameters and source spectra $`g(T_{\overline{p}})`$. We have calculated the quantities entering the factor $`\mathrm{{\rm Y}}`$ as described in Sec. II.

### V.1 Secondary antiprotons

Antiprotons in the Galaxy are also produced via standard interactions. Proton and helium cosmic rays interact with the interstellar hydrogen and helium nuclei, producing quarks and gluons that subsequently can hadronize into antiprotons. A calculation of this secondary antiproton flux has been done in Ref. pbar\_sec , to which we refer for details. Here we only emphasize the main results obtained in that work: i) The antiproton flux has been evaluated consistently by employing the propagation parameters as derived from a full and systematic analysis on stable nuclei PaperI ; ii) This secondary antiproton flux is in very good agreement with the data taken from the experiments bess bess95-97 ; bess98 , ams ams98 , caprice caprice  (see Fig. 14 in Ref. pbar\_susy ); iii) The uncertainty on the final flux due to propagation is about 20%, with a slight dependence on the energy. Another important source of uncertainty of order 20%-25% resides in the nuclear production cross sections, in particular when considering the interactions over the interstellar helium.

### V.2 Constraints on a primary antiproton source?

As discussed above, the secondary antiproton flux already provides a satisfactory agreement with current experimental data, and then no much room is left to primary contributions. This situation suggests that antiproton data could be used to place significant constraints on supersymmetric parameters. However, one has to notice that, as shown in pbar\_susy , the supersymmetric primary flux is affected by uncertainties much larger than those related to the secondary flux. This is due to the fact that the sources of the latter are located in the galactic disk. On the contrary, the relic neutralinos are expected to be distributed in the whole galactic halo and then produce an antiproton flux much more sensitive to the astrophysical parameters.

To show quantitatively how the experimental data could constrain the supersymmetric parameters, in Fig. 4 we display the antiproton flux evaluated at $`T_{\overline{p}}=0.23`$ GeV for a full scan of our supersymmetric model described in Sec. II. As expected, the scatter plot is prominent at small masses. Furthermore, it is remarkable that for $`m_\chi <`$ 25 GeV the scatter plot is funnel-shaped. The reason is the same as the one given above in connection with Fig. 1. The two panels of Fig. 4 correspond to two different sets of the propagation parameters. One, used in the left panel, is the set giving the best fit to the $`B/C`$ ratio, while the other, hereby denoted as the conservative set and used in the right panel provides the lowest (secondary and primary) antiproton fluxes. A spherical cored isothermal distribution for dark matter has been used. However, as mentioned before, a different choice does not significantly modifies the scatter plots. The shaded region denotes the amount of primary antiprotons which can be accommodated at $`T_{\overline{p}}=0.23`$ GeV without entering in conflict with the BESS experimental data bess95-97 ; bess98  and secondary antiproton calculations pbar\_sec .

From the right panel of Fig. 4 we conclude that, within the current astrophysical uncertainties, one cannot derive any constraint on the supersymmetric parameters, if one assumes a very conservative attitude in the selection of the propagation parameters. It is worth noticing that even within this choice, some supersymmetric configurations at very small masses are close to the level of detectability. As a further comment on the left panel of Fig. 4, we wish to stress that any further breakthrough in the knowledge of the astrophysical parameters would allow a significant exploration of small mass configurations, in case the conservative set of parameters is excluded. Should the effect of antiproton propagation turn out to be equivalent to the one obtained with the best fit set, the analysis of cosmic antiprotons would prove quite important for exploring very light neutralinos. This is particularly true for neutralino masses below $`<`$ 15 GeV, in view of the typical funnel shape displayed in the scatter plots.

## VI Upgoing muons at neutrino telescopes

Indirect evidence for WIMPs in our halo may be obtained at neutrino telescopes by measurements of the upgoing muons, which would be generated by neutrinos produced by pair annihilation of neutralinos captured and accumulated inside the Earth and the Sun others\_mu ; noi\_mu . The process goes through the following steps: capture by the celestial body of the relic neutralinos through a slow–down process due essentially to neutralino elastic scattering off the nuclei of the macroscopic body; accumulation of the captured neutralinos in the central part of the celestial body; neutralino–neutralino annihilation with emission of neutrinos; propagation of neutrinos (we have included the $`\nu _\mu `$$`\nu _\tau `$ vacuum oscillation effect with parameters: $`\mathrm{\Delta }m^2=3\times 10^3`$ eV2, $`\mathrm{sin}\theta =1`$) and conversion of their $`\nu _\mu `$ component into muons in the rock surrounding the detector; propagation and detection of the ensuing upgoing muons in the detector.

The various quantities relevant for the previous steps are calculated here according to the method described in Refs. noi\_mu , to which we refer for further details. We just recall that the neutrino flux due to the annihilation processes taking place in a distant source like the Sun, as a function of the neutrino energy $`E_\nu `$, is given by

$$\frac{dN_\nu }{dE_\nu }=\frac{\mathrm{\Gamma }_A}{4\pi d^{\mathrm{\hspace{0.17em}2}}}\underset{F,f}{}B_{\chi f}^{(F)}\frac{dN_{f\nu }}{dE_\nu },$$
(11)
where $`\mathrm{\Gamma }_A`$ is the annihilation rate inside the macroscopic body griest , $`d`$ is the distance from the source, $`F`$ denotes the $`\chi `$$`\chi `$ annihilation final states, $`B_{\chi f}^{(F)}`$ denotes the branching ratios into heavy quarks, $`\tau `$ leptons and gluons in the channel $`F`$; $`dN_{f\nu }/dE_\nu `$ is the differential distribution of the neutrinos generated by the hadronization of quarks and gluons and the subsequent hadronic semileptonic decays. The annihilation rate is given by $`\mathrm{\Gamma }_A=C/2\mathrm{tanh}^2(t_0/\tau _A)`$ griest , where $`t_0`$ is the age of the macroscopic body, $`\tau _A=(CC_A)^{1/2}`$, $`C_A`$ is the annihilation rate proportional to the neutralino–neutralino annihilation cross–section and $`C`$ denotes the capture rate per effective volume of the body. The capture rate $`C`$ is calculated here as in Refs. gould ; helio , the WIMP velocity distribution in the galactic inertial frame being described by a Maxwellian function with a dispersion velocity of 270 km s-1. We recall that the capture rate by the Earth is favored when the WIMP mass is close to the nuclear mass of one of the main chemical components of the Earth: Oxygen, Magnesium, Silicon in the mantle and Iron in the core gould . We have neglected the contributions of the light quarks directly produced in the annihilation process or in the hadronization of heavy quarks and gluons, because these light particles stop inside the medium (Sun or Earth) before their decay. For the case of the Sun we have also considered the energy loss of the heavy hadrons in the solar medium.

### VI.1 Neutrinos from the center of the Earth

In the left panel of Fig. 5 we show the expected upgoing muon flux integrated for $`E_\mu >`$1 GeV, as a function of $`m_\chi `$, and compared to the present experimental upper bounds on the same quantity from the experiments Superkamiokande sk , MACRO macro , and AMANDA amanda . For $`m_\chi <`$ 40 GeV the signal from the Earth presents several peaks due to resonant capture on Oxygen, Silicon and Magnesium (we recall that the process of capture on Earth is driven by the coherent neutralino–nucleus cross section). These elements are almost as abundant in Earth as Iron, which is the most important target nucleus for neutralino capture at higher masses. The dip at $`m_\chi M_Z/2`$ is due to the rescaling prescription of Eq. (3), since the resonance in the $`Z`$–exchange annihilation cross section reduces the relic abundance $`\mathrm{\Omega }_\chi h^2`$. Moreover, for $`m_\chi <`$ 25 GeV, the branching ratio $`B_{\chi \tau }^{(\stackrel{~}{\tau })}`$ to the $`\tau \overline{\tau }`$ final state, which is the one with the highest neutrino yield per annihilation, is suppressed. This last property is due to the fact that, in this range of $`m_\chi `$, the final state to $`b\overline{b}`$ in the annihilation cross section is required to be the dominant one in order to keep the relic abundance $`\mathrm{\Omega }_\chi h^2`$ below its cosmological upper bound lowmass . This, together with the fact that lower $`\chi `$ masses imply softer $`\nu `$ spectra that produce less $`\mu `$’s above threshold, explains why the up-going muon signal expected for light neutralinos ($`m_\chi <`$ 50 GeV) turns out to be below the level reached at higher masses.

It is worth noting that substantial modifications to the standard Maxwellian velocity distribution in the solar neighborhood have been considered in recent years, with two conflicting models. Damour and Krauss dk  proposed the existence of a solar-bound population, generated by WIMPs which scattered off the Sun surface and were then set, by perturbations from other planets, into orbits crossing the Earth, but not the Sun. The existence of this low-speed population would make the capture of WIMPs by the Earth very efficient with an ensuing dramatic enhancement of the up-going muon flux as compared to the standard case bdeku . On the contrary, Gould and Alam ga  used arguments based on calculations of asteroid trajectories to conclude that solar-bound WIMPs could evolve quite differently, inducing a significant suppression in the up-going muon flux from the Earth, as compared to the standard Maxwellian case, for WIMP masses larger than 150 GeV. A very recent re–analysis of this problem le  supports the conclusions of Gould and Alam, though with a less dramatic suppression effect. Our results have been derived using the standard Maxwellian distribution. The results of Ref. le  would not significantly alter our conclusion on the detectability of light neutralinos. The maximal fluxes are obtained for resonant capture on $`O`$, $`Si`$ and $`Mg`$ nuclei in the mantle: in this situation no suppression occurs. For neutralino masses away from the resonant condition, a reduction factor up to 0.8 for $`m_O<m_\chi <m_{Fe}`$ and up to 2–3 for $`m_\chi <10`$ GeV is possible le : however, in these cases, the upgoing muon flux is already very depressed, as it is shown in Fig. 5, left panel. In conclusion, using the standard Maxwellian distribution, the present measurements of up-going muons from the Earth put some constraints on neutralino configurations for masses above 50 GeV. For lighter neutralinos, explorations by neutrino telescopes would require a substantial increase in sensitivity while keeping a low energy threshold (close to 1 GeV). This in turn would imply a sizable extension of the telescope and a dense array of photomultipliers, which is certainly feasible, but very expensive.

### VI.2 Neutrinos from the Sun

In the right panel of Fig. 5 we show the up-going muon flux expected from the Sun, integrated for $`E_\mu >`$1 GeV, as a function of $`m_\chi `$. The signal is compared to the present experimental upper bounds on this flux coming from the experiments Superkamiokandesk , MACROmacro , AMANDAamanda  and Baksanbaksan . Also in this case the signal level turns out to be suppressed for $`m_\chi <`$ 50 GeV as compared to what is obtained at higher masses, the reasons of this behavior being the same as in the case of the Earth. On the other hand, the enhancement of the signal at $`m_\chi m_W`$ is due to a peculiar behavior of the neutralino–nucleon spin–dependent cross section, which drives the neutralino capture in the Sun (mainly on hydrogen). This cross section reaches its maximum whenever the $`Z`$–exchange channel dominates, and this requirement is verified when the neutralino–$`Z`$ coupling, proportional to the combination $`a_3^2a_4^2`$, is maximal spin . By numerical inspection we have verified that this last quantity is significantly peaked for $`m_\chi M_W`$. In this range of masses the annihilation channel to $`W^+W^{}`$ opens up and dominates the annihilation cross section. We conclude here that investigations of light neutralinos by up-going muons from the Sun do not provide favourable prospects.

## VII Conclusions and Perspectives

We have considered the most relevant indirect strategies for detecting the presence of relic neutralinos in our Galaxy through the products of their self-annihilation. This includes annihilations taking place directly in the galactic halo or inside celestial bodies (Earth and Sun). Our investigation has been performed in the frame of an effective supersymmetric model at the electroweak scale with no assumption on gaugino mass universality at the GUT scale. The range of the neutralino mass taken into consideration brackets a wide interval, from 6 GeV up to 500 GeV. While the low extreme is decided by the lower bound of 6 GeV, established in Ref. lowneu ; lowmass , the upper extreme of 500 GeV is chosen only for convenience. Actually, no model-independent upper limit for the neutralino mass is available, apart for a generic value of order of 1 TeV, beyond which the raison d’être of supersymmetry fades away. Though our calculations span over the wide range of the neutralino masses recalled above, our discussions were focused on light neutralinos, i.e. neutralinos with masses below 50 GeV: this value corresponds to the lower bound of $`m_\chi `$ when gaugino-mass unification is assumed. Indeed, light neutralinos are rarely considered in the literature, though their properties are quite interesting, as we already proved in connection with their cosmological properties and their detectability by current experiments of direct WIMP search lowneu ; lowmass ; lowdir . Thus the present paper is the natural continuation of our previous investigations on neutralinos of small mass. Different galactic dark matter distributions have been considered, from the cored isothermal one to the profiles obtained by numerical cosmological simulations in Refs. nfw ; moore , including also the most recent ones of Ref. navarro .

Now we summarize the main results of the present paper:

* For the $`\gamma `$ rays we have considered separately fluxes from the galactic center and high latitude regions, and compared our predictions with the EGRET data. Our numerical results have been provided employing as a reference DM distribution the NFW profile. We have shown that in this case the EGRET data at all angles do not put any constraints on the supersymmetric flux. The minimum gap between the theoretical predictions and the data occurs at light masses and is of almost one order of magnitude. We have discussed by how much this gap changes in terms of the dark matter distribution. This variation is relevant only for signals coming from the galactic center. For this sector, we have shown that using the $`r`$–dependent log–slope distribution of Ref. navarro  the gap between data and supersymmetric fluxes is reduced by a factor of 3 with respect to the NFW profile. Only profiles as steep as the Moore et al. one, disfavored by recent simulations, could exclude some light neutralino configurations. We have also shown that neutralinos of masses around 30–40 GeV could explain the EGRET excess in case of a significant enhancement effect as compared to the NFW distribution. A general word of caution concerns the fact that the background due to conventional cosmic rays production mechanisms still suffers from sizeable uncertainties.
* We have shown that in the case of cosmic antiprotons, no constraint on the supersymmetric parameters can be derived, if one assumes a very conservative attitude in the selection of the propagation parameters. However, it is remarkable that indeed the signal at very small masses is close to the level of detectability. Some breakthrough in the knowledge of the astrophysical parameters could allow a significant exploration of small mass configurations. This is particularly true for neutralino masses below about 15 GeV, in view of the typical funnel shape displayed in the scatter plots.
* The present measurements of up-going muons from the center of the Earth put some constraints on neutralino configurations for masses above 50 GeV. For lighter neutralinos, explorations by neutrino telescopes require a substantial increase in sensitivity with an energy threshold close to 1 GeV. Investigations of light neutralinos by up-going muons from the Sun are very disfavored.

We wish to recall that, according to the measurements of the HEAT Collaboration heat , the spectrum of the positron component of cosmic rays shows some enhancement between 7 and 20 GeV. This is only a mild effect which, as shown in Ref. heat , could be explained in terms of conventional secondary production mechanisms. Alternatively, some authors have interpreted this effect as a deviation from a pure secondary flux, which could be generated by neutralino self-annihilation kane ; be . This hypothesis, to be a viable one, requires that the neutralino flux is enhanced by a sizable factor. Since the measured positrons must be created in a region around the Earth of a radius of a few kpc origin , an enhancement would imply a significant dark matter overdensity in that region, with implications for the $`\overline{p}`$ signal. This scenario appears strongly model dependent, and as such not suitable for setting constraints to supersymmetric parameters.

Antideuterons in space as a signal of neutralino self-annihilation, which were shown in Ref.dbar  to be a very promising investigation means, will be considered in a forthcoming paper.

###### Acknowledgements.

We acknowledge Research Grants funded jointly by the Italian Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR), by the University of Torino and by the Istituto Nazionale di Fisica Nucleare (INFN) within the Astroparticle Physics Project. We thank the Referee for suggesting the inclusion of the upper bound on the $`BR(B_s\mu ^++\mu ^{})`$ decay in our analysis.


format: text

text:
# Statistical Challenges in Modern Astronomy

## I The glorious history of astronomy and statistics

Astronomy is perhaps the oldest observational science111The historical relationship between astronomy and statistics is described in references ,  and elsewhere. Our Astrostatistics monograph gives more detail and contemporary examples of astrostatistical problems .. The effort to understand the mysterious luminous objects in the sky has been an important element of human culture for at least $`10^4`$ years. Quantitative measurements of celestial phenomena were carried out by many ancient civilizations. The classical Greeks were not active observers but were unusually creative in the applications of mathematical principles to astronomy. The geometric models of the Platonists with crystalline spheres spinning around the static Earth were elaborated in detail, and this model endured in Europe for 15 centuries. But it was another Greek natural philosopher, Hipparchus, who made one of the first applications of mathematical principles that we now consider to be in the realm of statistics. Finding scatter in Bablylonian measurements of the length of a year, defined as the time between solstices, he took the middle of the range – rather than the mean or median – for the best value.

This is but one of many discussions of statistical issues in the history of astronomy. Ptolemy estimated parameters of a non-linear cosmological model using a minimax goodness-of-fit method. Al-Biruni discussed the dangers of propagating errors from inaccurate instruments and inattentive observers. While some Medieval scholars advised against the acquisition of repeated measurements, fearing that errors would compound rather than compensate for each other, the usefulnes of the mean to increase precision was demonstrated with great success by Tycho Brahe.

During the 19th century, several elements of modern mathematical statistics were developed in the context of celestial mechanics, where the application of Newtonian theory to solar system phenomena gave astonishingly precise and self-consistent quantitative inferences. Legendre developed $`L_2`$ least squares parameter estimation to model cometary orbits. The least-squares method became an instant success in European astronomy and geodesy. Other astronomers and physicists contributed to statistics: Huygens wrote a book on probability in games of chance; Newton developed an interpolation procedure; Halley laid foundations of actuarial science; Quetelet worked on statistical approaches to social sciences; Bessel first used the concept of ”probable error”; and Airy wrote a volume on the theory of errors.

But the two fields diverged in the late-19th and 20th centuries. Astronomy leaped onto the advances of physics – electromagnetism, thermodynamics, quantum mechanics and general relativity – to understand the physical nature of stars, galaxies and the Universe as a whole. A subfield called “statistical astronomy” was still present but concentrated on rather narrow issues involving star counts and Galactic structure Trumpler53 . Statistics concentrated on analytical approaches. It found its principle applications in social sciences, biometrical sciences and in practical industries (e.g., Sir R. A. Fisher’s employment by the British agricultural service).

## II Statistical needs of astronomy today

Contemporary astronomy abounds in questions of a statistical nature. In addition to exploratory data analysis and simple heuristic (usually linear) modeling common in other fields, astronomers also often interpret data in terms of complicated non-linear models based on deterministic astrophysical processes. The phenomena studied must obey known behaviors of atomic and nuclear physics, gravitation and mechanics, thermodynamics and radiative processes, and so forth. ‘Modeling’ data may thus involves both the selection of a model family based on an astrophysical understanding of the conditions under study, and a statistical effort to find parameters for the specified model. A wide variety of issues thus arise:

* Does an observed group of stars (or galaxies or molecular clouds or $`\gamma `$-ray sources) constitute a typical and unbiased sample of the vast underlying population of similar objects?
* When and how should we divide/classify these objects into 2, 3 or more subclasses?
* What is the intrinsic physical relationship between two or more properties of a class of objects, especially when confounding variables or observational selection effects are present?
* How do we answer such questions in the presence of observations with measurements errors and flux limits?
* When is a blip in a spectrum (or image or time series) a real signal rather than a random event from Gaussian (or often Poissonian) noise or confounding variables?
* How do we interpret the vast range of temporally variable objects: periodic signals from rotating stars or orbiting extrasolar planets, stochastic signals from accreting neutron stars or black holes, explosive signals from magnetic reconnection flares or $`\gamma `$-ray bursts?
* How do we model the points in 2, 3, …, 6-dimensional points representing photons in an image, galaxies in the Universe, Galactic stars in phase space?
* How do we quantify continuous structures seen in the sky such as the cosmic microwave background, the interstellar and intergalactic gaseous media?
* How do we fit astronomical spectra to highly non-linear astrophysical models based on atomic physics and radiative processes, including confidence limits on the best-fit parameters?

From a superficial examination of the astronomical literature222Such bibliometric measures are easily accomplished as the entire astronomical research literature is on-line (in full text at subscribing institutions) through the NASA-supported Astrophysics Data System, http://adsabs.harvard.edu/abstract\_service.html., we can show that such questions are very common today. Of $`15,000`$ refereed papers published annually, 1% have “statistics” or “statistical” in their title, 5% have “statistics’ in their abstract, 10% treat time-variable objects, $`510`$% (est.) present or analyze multivariate datasets, and $`510`$% (est.) fit parametric models. Accounting for overlaps, we roughly estimate that around $`3,000`$ distinct studies each year require non-trivial statistical methodologies. Roughly 10% of these are principally involved with statistical methods; indeed, some of these purport to develop new methods or improve on established ones.

## III Astrostatistics today

We thus find that astronomy and astrophysics today requires a vast range of statistical capabilities. In statistical jargon, it helps for astronomers to know something about: sampling theory, survival analysis with censoring and truncation, measurement error models, multivariate classification and analysis, harmonic and autoregressive time series analysis, wavelet analysis, spatial point processes and continuous surfaces, density estimation, linear and non-linear regression, model selection, and bootstrap resampling. In some cases, astronomers need combinations of methodologies that have not yet been fully developed (§7 below).

Faced with such a complex of challenges, mechanical exposure to a wider variety of techniques is a necessary but not sufficient prerequisite for high-quality statistical analyses. Astronomers also need to be imbued with established principles of statistical inference; e.g., hypothesis testing and parameter estimation, nonparametric and parametric inference, Bayesian and frequentist approaches, and the assumptions underlying and applicability conditions for any given statistical method.

Unfortunately, we find that the majority of the thousands of astronomical studies requiring statistical analyses use a very limited set of classical methods. The most common tools used by astronomers are: Fourier transforms for temporal analysis (developed by Fourier in 1807), least squares regression and $`\chi ^2`$ goodness-of-fit (Legendre in 1805, Pearson in 1900, Fisher in 1924), the nonparametric Kolmogorov-Smirnov 1- and 2-sample nonparametric tests (Kolmogorov in 1933), and principal components analysis for multivariate tables (Hotelling in 1936).

Even traditional methods are often misused. Feigelson & Babu  found that astronomers use interchangeably up to 6 different fits for bivariate linear least squares regression: ordinary least squares (OLS), inverse regression, orthogonal regression, major axis regression, the OLS mean, and the OLS bisector. Not only did this lead to confusion in comparing studies (e.g., in measuring the expansion of the Universe via Hubble’s constant, $`H_o`$), but astronomers did not realize that the confidence intervals on the fitted parameters can not be correctly estimated with standard analytical formulae. Similarly, Protassov et al.  found that the majority of astronomical applications of the $`F`$ test, or more generally the likelihood ratio test, are inconsistent with asymptotic statistical theory.

But, while the average astronomical study is limited to often-improper usage of a limited repertoire of statistical methods, a significant tail of outliers are much more sophisticated. The maximization of likelihoods, often developed specially for the problem at hand, is perhaps the most common of these improvements. Bayesian approaches are also becoming increasingly in vogue.

In a number of cases, sometimes buried in technical appendices of observational papers, astronomers independently develop statistical methods. Some of these are rediscoveries of known procedures; for example, Avni et al.  and others recovered elements of survival analysis for treatments of left-censored data arising from nondetections of known objects. Some are quite possibly mathematically incorrect; such as various revisions to $`\chi ^2`$ for Poissonian data that assume the resulting statistic still follows the $`\chi ^2`$ distribution. On rare occasions, truly new and correct methods have emerged; for example, astrophysicist Lynden-Bell  discovered the maximum-likelihood estimator for a randomly truncated dataset, for which the theoretical validity was later established by statistician Woodroofe .

A growing group of astronomers, recognizing the potential for new liaisons with the accomplishments of modern statistics, have promoted astrostatistical innovation through cross-disciplinary meetings and collaborations. Fionn Murtagh, an applied mathematician at Queen’s University (Belfast) with long experience in astronomy, and his colleagues have run conferences and authored many useful monographs (e.g., , ,  and ). We at Penn State have run a series of Statistical Challenges in Modern Astronomy meetings with both communities in attendance (e.g.,  and ). Alanna Connors has organized brief statistics sessions at large astronomy meetings, and we have organized brief astronomy sessions at large Joint Statistical meetings. We wrote a short volume called Astrostatistics  intended to familiarize scholars in one discipline with relevant issues in the other discipline. Other series conferences are devoted to technical issues in astronomical data analysis but typically have limited participation by statisticians. These include the dozen Astronomical Data Analysis Software and Systems (e.g., ), several Erice workshops on Data Analysis in Astronomy (e.g., ), and the new SPIE Astronomical Data Analysis conferences (e.g., ).

Most importantly, several powerful astrostatistical research collaborations have emerged. At Harvard University and the Smithsonian Astrophysical Observatory, David van Dyk worked with scientists at the $`Chandra`$333The $`Chandra`$ X-ray Observatory is one of NASA’s Great Observatories. It was launched in 1999 with a total budget around $2 billion. X-ray Center on several issues, particularly Bayesian approaches to parametric modeling of spectra in light of complicated instrumental effects. At Carnegie Mellon University and the University of Pittsburgh, the Pittsburgh Computational Astrophysics group addressed several issues, such as developing powerful techniques for multivariate classification of extremely large datasets and applying nonparametric regression methods to cosmology. Both of these groups involved academics, researchers and graduate students from both fields working closely for several years to achieve a critical mass of cross-disciplinary capabilities.

Other astrostatistical collaborations must be mentioned. David Donoho (Statistics at Stanford University) works with Jeffrey Scargle (NASA Ames Research Center) and others on applying advanced wavelet methods to astronomical problems. James Berger (Statistics at Duke University) has worked with astronomers William Jefferys (University of Texas), Thomas Loredo (Cornell University), and Alanna Connors (Eureka Inc.) on Bayesian methodologies for astronomy. Bradley Efron (Statistics at Stanford University) has worked with astrophysicist Vehé Petrosian (also at Stanford) on survival methods for interpreting $`\gamma `$-ray bursts. Philip Stark (Statistics at University of California, Berkeley) has collaborated with solar physicists in the $`GONG`$ program to improve analysis of oscillations of the Sun (helioseismology). More such collaborations exist in the U.S., Europe and elsewhere.

## IV The Virtual Observatory: A new imperative for astrostatistics

A major new trend is emerging in observational astronomy with the production of huge, uniform, multivariate databases from specialized survey projects and telescopes444An enormous collection of catalogs, and some of the underlying imaging and spectral databases, are already available on-line. Access to many catalogs is provided by Vizier (http://vizier.u-strasbg.fr). The NASA Extragalactic Database (NED, http://ned.ipac.caltech.edu), SIMBAD stellar database (http://simbad.u-strasbg.fr), and ADS (footnote 2) give integrated access to many catalogs and bibliographic information. Raw data are available from all U.S. space-based observatories; see, for example, the Multi-mission Archive at Space Telescope (MAST, http://archive.stsci.edu) and High Energy Astrophysics Science Archive Research Center (HEASARC, http://heasarc.gsfc.nasa.gov).. But they are heterogeneous in character, reside at widely dispersed locations, and accessed through different database systems. Examples include:

1. $`10^810^9`$-object catalogs of stars and stellar extragalactic objects (i.e., quasars). These include the all-sky photographic optical USNO-B1 catalog, the all-sky near-infrared 2MASS catalog, and the wide-field Sloan Digital Sky Survey (SDSS). Five to ten photometric values, each with measured heteroscedastic measurement errors (i.e., different for each data point), are available for each object.
2. $`10^510^6`$-galaxy redshift catalogs from the 2-degree Field (2dF) and SDSS spectroscopic surveys. The main goal is characterization of the hierarchical, nonlinear and anisotropic clustering of galaxies in a 3-dimensional space. But the datasets also include spectra for each galaxy each with $`10^3`$ independent measurements.
3. $`10^510^6`$-source catalogs from various multiwavelength wide-field surveys such as the NRAO Very Large Array Sky Survey in one radio band, the InfraRed Astronomical Satellite Faint Source catalog in four infrared bands, the Hipparcos and Tycho catalogs of star distances and motions, and the X-ray Multimirror Mission Serendipitous Source Catalogue in several X-ray bands now in progress. These catalogs are typically accompanied by large image libraries.
4. $`10^210^4`$-object samples of well-characterized pre-main sequence stars, binary stars, variable stars, pulsars, interstellar clouds and nebulae, nearby galaxies, active galactic nuclei, gamma-ray bursts and so forth. There are dozens of such samples with typically $`1020`$ catalogued properties and often with accompanying 1-, 2- or 3-dimensional images or spectra.
5. Perhaps the most ambitious of such surveys is the planned Large-aperture Synoptic Survey Telescope (LSST) which will survey much of the entire optical sky every few nights. It is expected to generate raw databases in excess of 10 PBy (petabyte) and catalogs with $`10^{10}`$ entries.

An international effort known as the Virtual Observatory (VO) is now underway to coordinate and federate these diverse databases, making them readily accessible to the scientific user (Szalay01, ; Brunner01, ). Considerable progress is being made in the establishment of the necessary data and metadata infrastructure and standards, interoperability issues, data mining, and technology demonstration prototype services555See http://www.ivoa.net and /http://us-vo.org for entry into Virtual Observatory projects.. But scientific discovery requires more than effective recovery and distribution of information. After the astronomer obtains the data of interest, tools are needed to explore the datasets. How do we identify correlations and anomalies within the datasets? How do we classify the sources to isolate subpopulations of astrophysical interest? How do we use the data to constrain astrophysical interpretation, which often involve highly non-linear parametric functions derived from fields such as physical cosmology, stellar structure or atomic physics? These questions lie under the aegis of statistics.

A particular problem relevant to statistical computing is that, while the speed of CPUs and the capacity of inexpensive hard disks rise rapidly, computer memory capacities grow at a slower pace. Combining the largest optical/near-infrared object catalogs today produces a table with $`>1`$ billion objects and up to a dozen columns of photometric data. Such large datasets effectively preclude use of all standard multivariate statistical packages and visualization tools (e.g., R and GGobi) which are generally designed to place the entire database into computer memory. Even sorting the data to produce quantiles may be computational infeasible.

The Virtual Observatory of the 21st century thus presents new challenges to statistical capability in two ways. First, some new methodological developments are needed (§V). Second, efficient access to both new and well-established statistical methods are needed. No single existing software package can provide the vast range of needed methods. We are now involved in developing a prototype system called VOStat to provide statistical capabilities to the VO astronomer. It is based on concepts of Web services and distributed Grid computing. Here, the statistical software and computational resources, as well as the underlying empirical databases, may have heterogeneous structures and can reside at distant locations.

## V Some grand methodological challenges for the coming decade

While it is risky to prognosticate the directions of future research, and judgments will always differ regarding the relative importance of research goals, we can outline a few “grand challenges” for astrostatistical research for the next decade or two.

### V.1 Multivariate analysis with measurement errors and censoring

Traditional multivariate analysis is designed mainly for applications in the social and human sciences where the sources of variance are largely unknowable. Measurement errors are usually ignored, or are considered to be exogenous variables in the parametric models Fuller87 . But astrophysicists often devote as much effort to precise determination of their errors as they devote to the measurements of the quantities of interest. The instruments are carefully calibrated to reduce systematic uncertainties, and background levels and random fluctuations are carefully evaluated to determine random errors. Except in the simple case of bivariate regression (Boggs87, ; Feigelson92, ; Akritas96, ), this information on measurement errors is usually squandered.

While heteroscedastic measurement errors with known variances is common in all physical sciences, only astronomy frequently has nondetections when observations are made at new wavelengths of known objects. These are datapoints where the signal lies below (say) 3 times the noise level. Here again, modern statistics has insufficient tools. Survival analysis for censored data assumes that the value below which the data point must lie is known with infinite precision, rather than being generated from a distribution of noise. Astronomer Herman Marshall  makes an interesting attempt to synthesize measurement errors and nondetections, but statistician Leon Gleser  argues that he has only recovered Fisher’s failed theory of fiducial distributions. Addressing this issue in a self-consistent statistical theory is a profound challenges that lies at the heart of interpreting the data astronomers obtain at the telescope.

### V.2 Statistical inference and visualization with very-large-N datasets

The need for computational software for extremely large databases – multi-terabyte image and spectrum libraries and multi-billion object catalogs – is discussed in section 4. A suite of approximate methods based on flowing data streams or adaptive sampling of large datasets resident on hard disks should be sought. Visualization methods involving smoothing, multidimensional shading and variable transparency, should be brought into the astronomer’s toolbox. Here, considerable work is being conducted by computer scientists and applied mathematicians in other applied fields so that independent development by astrostatisticians might not be necessary to achieve certain goals.

### V.3 A cookbook for construction of likelihoods and Bayesian computation

While the concepts of likelihoods and their applications in maximum likelihood estimation, Bayes Theorem and Bayes factors are becoming increasingly well-known in astronomical research, the applications to real-life problems is still an art for the expert rather than a tool for the masses. Part of the problem is conceptual; astronomers need training in how to construct likelihoods for familiar parametric situations (e.g., power law distributions or a Poisson process). Part of the problem is computational; astronomers need methods and software for the oft-complex computations. Many such methods, such as Markov chain Monte Carlo, are already well-established and can be directly adopted for astronomy . For example, astronomers are often not fully aware of the broad applicability of the EM Algorithm for maximizing likelihoods McLachlan97 666The seminal study of the EM Algorithm is Dempster, Laird & Rubin in 1977 , which is one of the most frequently cited papers in statistics. However, the method was independently derived three years earlier by astronomer Leon Lucy  as an “iterative technique for the rectification of observed distributions” based on Bayes’ Theorem. This study is widely cited in the astronomical literature; its most frequent application is in image deconvolution where it is known as the Lucy-Richardson algorithm..

### V.4 Links between astrophysical theory and wavelets

Wavelet analysis has become a powerful and sophisticated tool for the study of features in data. Originally intended mainly for modelling time series, astronomers also use it increasingly for spatial analysis of images Freeman02 ; Scargle98 . In some ways it can be viewed as a generalization of Fourier analysis in which the basis function need not be sinusoidal in shape and, most importantly, the pattern need not extend over the entire dataset. Wavelets are thus effective in quantitatively describing complicated overlapping structures on many scales, and can also be used for signal denoising and compression. In addition, wavelets have a strong mathematical foundation.

Despite its increasing popularity in astronomical applications, wavelet analysis suffers a profound limitation in comparison with Fourier analysis. A peak in a Fourier spectrum is immediately interpretable as a vibrational, rotational or orbital rotation of solid bodies. A bump or a continuum slope in a wavelet decomposition often has no analogous physically intuitive interpretation. We therefore recommend that astrophysicists seek links between physical theory – often involving continuous media such as turbulent plasmas in the interstellar medium and hierarchical structure formation in the early Universe – and wavelets. One fascinating example is the demonstration that the wavelet spectrum and Lyapunov exponent of the quasi-periodic X-ray emission from Sco X-1, which reflects the processes in an accretion disk around a neutron star, exhibit a transient chaotic behavior similar to that of water condensing and dripping onto an automobile windshield or a dripping handrail (Young96, ).

### V.5 Time series models for astrophysical phenomena

The quasi-periodic oscillation of Sco X-1 is only one of many examples of complex accretional behavior onto neutron stars and black holes seen in X-ray and $`\gamma `$-ray astronomy. The accreting Galactic black hole GRS 1915+105 exhibits a bewildering variety of distinct states of stochastic, quasi-periodic and explosive behaviors. The prompt emission from gamma-ray bursts show a fantastic diversity of temporal behaviors from simple smooth fast-rise-exponential-decays to stochastic spiky profiles. Violent magnetic reconnection flares on the surfaces of the Sun and other magnetically active stars also show complex behaviors. Many of these datasets are multivariate with time series available in several spectral bands often showing lags or hardness ratio variations of astrophysical interest.

There are also important astronomical endeavors which seek astrophysically interesting signals amidst the oft-complex noise characteristics of the detectors. The Arecibo, Parkes and VLA radio telescopes, for example, conduct searches for new radio pulsars or for extraterrestrial intelligences in nearby planetary systems. The Laser Interferometer Gravitational-Wave Observatory (LIGO) and related detectors search for both continuing periodic signals and brief bursts from perturbations in space-time predicted by Einstein’s General Relativity. Here the signals sought are orders of magnitude fainter than instrumental variations.

## VI Infrastructure needed to advance astrostatistics

The current quality of statistical analyses in astronomical research often begs for improvment. There is both inadequate research on important new challenges (§V) and inadequate application of known advanced methods to astronomical problems (§3). Astronomy clearly needs needs a strong and rapid surge of energy in statistical expertise. Three types of activities should be promoted:

In the U.S., the typical curriculum leading to a career in astronomical research requires zero or one course in statistics at the undergraduate level, and zero at the graduate level. Analogously, the curriculum of statisticians includes virtually no coursework in astronomy or other physical science. While statisticians can learn basics from “Astronomy 101” courses given at all universities, the statistical training of astronomers is not as easily accomplished. New curricular products summarizing the applicable statistical subfields, short training workshops for graduate students and young scientists, and effective statistical consulting are all needed.

While several astrostatistical research groups are making exciting progress (§3), the total effort is too small to impact the bulk of astronomical research. Very roughly, astrostatistical funding is currently $1M of the $1B spent annually on astronomical research. This fraction is far below that spent in biomedical or other non-physical-science fields. Though top academic leaders of statistics have expressed great enthusiasm for astronomy and astrostatistics, we can not pull them away from biostatistics and business applications without a major increase in funding. We might seek, for example, $`1020`$ cross-disciplinary research groups active at any one time at the end of a decade’s growth.

For various policy and cultural reasons, astronomers rarely purchase the large commercial statistical software packages, preferring to write their own software as needs arise. This approach has contributed to the narrow methodological scope of astronomical research. Avenues for improving this situation are emerging. R is a large statistical software package with the flexible command-line interface preferred by astronomers that has recently emerged (http://www.r-project.org). A wide variety of specialized packages and codes are also available on-line (http://www.astro.psu.edu/statcodes). The new Web services concept being developed within the context of a Virtual Observatory permits coordinated access to heterogeneous software developed specifically for astronomical applications.

At Penn State, we are in the early stages of developing a Center for Astrostatistics to help attain these goals (http://www.astrostatistics.psu.edu). This is an inter-disciplinary Center to serve the astronomy and statistics communities around the nation and worldwide, seeking to bring advances in statistics into the toolbox of astronomy and astrophysics. The Center’s Web site will maintain the popular StatCodes, build an instructional library of $`R`$ programs, coordinate with the nascent VOStat Web service, and develop an archive of annotated links to selected statistical literature applicable to astronomy (and vice versa). The site is also planned to include tutorial handbooks and curricular products developed specifically for astrostatistical needs.

###### Acknowledgements.

We thank the referee, Peter Shawhan (Caltech), for very helpful improvements to the paper. This work was supported by NSF grant DMS-0101360.