text:
# Criticality and convergence in Newtonian collapse

## I introduction

The Newtonian isothermal gas is a relevant description for important astrophysical systems. It is indeed a good description for cold molecular clouds found in the galactic interstellar medium. In these clouds the typical cooling time is usually much shorter than the dynamical time; as a result they are in thermal equilibrium with whatever steady heating source is present, for example the blackbody radiation of the cosmic microwave background, extragalactic ultraviolet radiation, or local stars. The molecular clouds can be described as isothermal. Since they are the prime location for star formation, their physics determines the whole star formation history in the galaxy. A detailed understanding of the physics of these systems will lead to predictions for the properties of star formation.

There are two possible theoretical approaches for the self-gravitating Newtonian gas: statistical physics and dynamics. The traditional statistical physics approach deals with equilibrium configurations of the system and gave the first description of Jeans gravitational instability which is a fundamental behavior in the physics of self-gravitating systems (Chandrasekhar chandrasekhar1939 , Padmanabhan padmanabhan1990 ). Using the modern tools of statistical field theory and renormalization groups, de Vega, Sanchez and Combes devega1996  and Semelin et al. semelin1999  give a theoretical basis to the observed self-similar properties of the gas. These statistical methods are useful, generally speaking, to describe systems undergoing second order phase transition and exhibiting critical behavior. Investigating the critical behavior of the gas is also the subject of the present work, where we choose to describe the system through dynamics rather than statistical physics and through the renormalization group.

The basic behavior of a clump of isothermal gas is the following. In a sufficiently dense and/or cold configuration, the gas will collapse due to its self-gravity, while in a diffuse and/or hot one it will expand due to its internal pressure forces. A critical behavior appears at the transition between these two regimes.

In fact, Choptuik choptuik1993  showed the existence of a so-called critical phenomenon in the spherical collapse of a massless scalar field by general relativistic numerical simulations. The simulations were performed for a one-parameter family of initial data sets. He found a self-similar solution with critical properties at the threshold of black hole formation. In this first case, the self-similarity holds for a discrete set of transformations only (periodicity in logarithmic time). He showed that the mass of the formed black hole follows a scaling law, and that the phenomenon is universal in the sense that it does not depend on the one-parameter family of initial data sets. Soon after that, similar critical phenomena were observed in the spherical collapse of radiation fluid by Evans and Coleman ec1994 . In this second case, self-similarity holds for a continuous set of transformations. Koike, Hara and Adachi kha1995  gave a clear picture of these phenomena by using the renormalization group approach. Recently, critical phenomena have been found in perfect fluids with more general equation of state maison1996 ; kha1999 ; nc2000  and in a variety of systems. Astrophysical implications of the critical phenomena have been studied in the scenario of primordial black hole formation by Yokoyama yokoyama1998  and Niemeyer and Jedamzik nj1998 . See gundlach1999  for a recent review of critical phenomena in gravitational collapse.

The solution with critical properties is only one of several known self-similar solutions. See Carr and Coley cc1999  for a recent review of self-similar solutions in general relativity. It was conjectured by Carr carr1993  that a generic solution gradually evolves toward a self-similar form, which suggest that one of the self-similar solutions may be an attractor solution. This is called the self-similarity hypothesis. In Newtonian gravity, a self-similar solution, the so-called Larson-Penston solution penston1969 ; larson1969  is known to be an attractor solution in the spherical collapse of isothermal gas hn1997 ; hm2000 ; ti1999 . In general relativity, Harada and Maeda hm2001  showed by numerical simulations that the general relativistic counterpart of the Larson-Penston solution, which was discovered by Ori and Piran op1987 , is an attractor solution in the spherical collapse of a perfect fluid with $`P=kϵ`$ ($`0<k0.03`$), where $`P`$ and $`ϵ`$ are the pressure and energy density, respectively. This convergence is quite important in general relativity because it provides a strong counterexample against the cosmic censorship conjecture in spherical collapse op1987 ; op1990 ; harada1998 ; hm2001 . See hin2002  for a recent review of cosmic censorship and related topics. In the context of the self-similarity hypothesis, the stability criterion for general relativistic self-similar solutions against kink mode was recently obtained by Harada harada2001 , through a generalization of the Newtonian analysis op1988 .

The present paper shows that critical phenomena exist in Newtonian gravity also. We concentrate on a spherically symmetric system of isothermal gas in Newtonian gravity. The existence of a critical phenomenon was predicted by Harada and Maeda hm2001  and Maeda and Harada mh2001 . In these papers, based on the renormalization group analysis, it was predicted that a critical phenomenon is associated with the first member of Hunter’s family of continuous self-similar solutions hunter1977 . The criticality appears in the scaling law for the collapsed mass, which is satisfied for continuous values of the mass. The value of the critical exponent was predicted to be $`0.10567`$. The critical phenomenon observed in this paper agrees very well with the above predictions, and gives a precise meaning to the collapsed mass in Newtonian gravity. We have also observed that the scaling law holds over many orders of magnitude in terms of the parameter.

The organization of this paper is as follows. In Sec. II, we apply the renormalization group analysis to the Newtonian isothermal gas system and draw the expected properties for the critical phenomenon. In Sec. III, we present the results from numerical simulations of Newtonian collapse. In Sec. IV we study the response of the system to temperature fluctuations and discuss astrophysical applications. In Sec. V, we summarize the paper.

## II Renormalization group analysis

In this section, we review the renormalization group analysis of the critical behavior in a spherical system of Newtonian isothermal gas. Part of the following discussion is based on Maeda and Harada mh2001 .

The dynamics of a spherically symmetric isothermal gas system is described in the following set of equations:

$`{\displaystyle \frac{\rho }{t}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}(r^2\rho v)=0,`$ (1)
$`{\displaystyle \frac{}{t}}(\rho v)+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}(r^2\rho v^2)=c_\mathrm{s}^2{\displaystyle \frac{\rho }{r}}\rho {\displaystyle \frac{Gm}{r^2}},`$ (2)
$`{\displaystyle \frac{m}{t}}+v{\displaystyle \frac{m}{r}}=0,`$ (3)
$`{\displaystyle \frac{m}{r}}=4\pi r^2\rho ,`$ (4)
where $`\rho =\rho (t,r)`$, $`v=v(t,r)`$, $`m=m(t,r)`$, $`c_\mathrm{s}`$ and $`G`$ denote the density, radial velocity, total mass inside the radial coordinate, sound speed, and gravitational constant, respectively. The sound speed $`c_\mathrm{s}`$ is related to the temperature $`T`$ of the gas as

$$c_\mathrm{s}=(kT/m)^{1/2},.$$
(5)
where $`k`$ and $`m`$ are the Boltzmann constant and mass of each particle, respectively.

In order to apply the renormalization group approach, it is convenient to redefine the physical quantities as follows:

$$\widehat{r}\frac{r}{c_\mathrm{s}},\widehat{v}\frac{v}{c_\mathrm{s}},\widehat{m}\frac{m}{c_\mathrm{s}^3}.$$
(6)
In terms of these new quantities, the basic equations reduce to the following form:

$`{\displaystyle \frac{\rho }{t}}+{\displaystyle \frac{1}{\widehat{r}^2}}{\displaystyle \frac{}{\widehat{r}}}(\widehat{r}^2\rho \widehat{v})=0,`$ (7)
$`{\displaystyle \frac{}{t}}(\rho \widehat{v})+{\displaystyle \frac{1}{\widehat{r}^2}}{\displaystyle \frac{}{\widehat{r}}}(\widehat{r}^2\rho \widehat{v}^2)={\displaystyle \frac{\rho }{\widehat{r}}}\rho {\displaystyle \frac{G\widehat{m}}{\widehat{r}^2}},`$ (8)
$`{\displaystyle \frac{\widehat{m}}{t}}+\widehat{v}{\displaystyle \frac{\widehat{m}}{\widehat{r}}}=0,`$ (9)
$`{\displaystyle \frac{\widehat{m}}{\widehat{r}}}=4\pi \widehat{r}^2\rho .`$ (10)
We can see that the above hatted system does not contain any dependence on the temperature (or sound speed), which simplifies the analysis.

The initial conditions for the bare system are specified by the density and velocity profiles at the initial time $`t=t_\mathrm{i}<0`$ ($`t=0`$ is the time when the density singularity occurs in the center or, more physically, the time of core formation). Suppose the initial density and velocity profiles are given by functions $`\rho _\mathrm{i}(r)`$ and $`v_\mathrm{i}(r)`$,

$$\rho (t_\mathrm{i},r)=\rho _\mathrm{i}(r),v(t_\mathrm{i},r)=v_\mathrm{i}(r),$$
(11)
respectively. Using the relations of the hatted quantities to the bare quantities, we specify the initial conditions for the hatted system by

$$\widehat{\rho }(t_\mathrm{i},\widehat{r})=\rho _\mathrm{i}(c_\mathrm{s}\widehat{r}),\widehat{v}(t_\mathrm{i},\widehat{r})=\frac{v_\mathrm{i}(c_\mathrm{s}\widehat{r})}{c_\mathrm{s}}.$$
(12)
Therefore, the parameter-dependent system (1)–(4) with unique initial conditions is equivalent to the parameter-independent system (7)–(10) with a one-parameter family of initial conditions. The parameter is $`c_\mathrm{s}`$ or $`T`$.

Although we avoid to fully repeat the analysis by Koike, Hara and Adachi kha1995 , we need to review it briefly. We introduce the variables

$$\tau \mathrm{ln}(t),x\mathrm{ln}\frac{\widehat{r}}{t}.$$
(13)
$`\tau `$ is the scaling variable in the renormalization group transformation. Suppose a self-similar solution $`H_{\mathrm{ss}}(x)`$ exists, with a unique relevant (unstable) mode, which turns out to be the critical solution. This self-similar solution can be considered as a fixed point of the renormalization group transformation. The uniqueness of the unstable mode implies that the fixed point $`H_{\mathrm{ss}}`$ has a stable manifold of codimension one. This manifold is referred to as a critical surface. Any one-parameter family $`\{H_{(p)}(x)|p𝐑\}`$ of initial data sets with parameter $`p`$ generically has an intersection $`H^{}`$ with the critical surface of the fixed point $`H_{\mathrm{ss}}`$. $`H^{}`$ will be driven to $`H_{\mathrm{ss}}`$ under the renormalization group transformation, and hence $`H^{}`$ is the initial data set with the critical value $`p^{}`$. We consider an initial data $`H_\mathrm{i}(x)`$ in the one-parameter family, close to $`H^{}`$, i.e.,

$$h(0,x)=H_\mathrm{i}(x)=H^{}(x)+ϵF(x),$$
(14)
where $`h(\tau ,x)`$ is a solution of the partially differential equations and $`ϵ=pp^{}`$. Then, for large $`\tau _0`$, we have

$$h(\tau _0,x)H_{\mathrm{ss}}(x)+ϵe^{\kappa \tau _0}F_{\mathrm{rel}}(x),$$
(15)
where $`\kappa `$ and $`F_{\mathrm{rel}}(x)`$ are the eigenvalue and eigenfunction of the relevant mode, respectively. Now we choose $`\tau _0`$ so that the first term and the second term are of the same order, i.e.,

$$ϵe^{\kappa \tau _0}=O(1).$$
(16)
At $`\tau _0`$, the solution is deviating from the critical solution. If the collapse bounces to expansion, which will be the case for a subcritical collapse, the above condition obviously holds for the velocity field at the moment of the bounce. On the other hand, if the finite mass near the center shrinks considerably faster than for the critical collapse, which will be the case for the supercritical case, we have a $`\tau _0`$ for which the above condition also holds. Since the deviation from the critical collapse grows at $`x=O(1)`$, the relation $`\widehat{r}=e^{x\tau }`$ implies that the radius of the core collapse or bounce is given by $`\widehat{r}=O(e^{\tau _0})=O(ϵ^\gamma )`$, where $`\gamma =1/(\text{Re}\kappa )`$.

In order to see clearly that self-similar solutions exist in the present hatted system, we introduce the following quantities:

$$U\widehat{v},P4\pi \widehat{r}^2\rho ,M\frac{G\widehat{m}}{t}.$$
(17)
Then, the equations transform to

$`\dot{P}+(1+zU)P^{}+zPU^{}=0,`$ (18)
$`(\dot{U}P+U\dot{P})(U+zU^2+z)P^{}(1+2zU)PU^{}2zP+MPz^2=0,`$ (19)
$`\dot{M}M^{}M+PU=0,`$ (20)
$`zM^{}=P,`$ (21)
where the dot and prime denote the partial derivatives with respect to $`\tau `$ and $`x`$, respectively, and $`ze^x`$. The self-similar solutions are obtained by numerical integration of the ordinary differential equations which follow from the above equations with the assumption $`P=P(x)`$, $`U=U(x)`$ and $`M=M(x)`$. If we assume the existence of the critical solution, using the relation

$$m=\frac{c_\mathrm{s}^3}{G}M(x)e^x\widehat{r},$$
(22)
we obtain the scaling law for the collapsed mass or the “core” mass for the supercritical collapse,

$$m_{\mathrm{core}}|pp^{}|^\gamma .$$
(23)
For the subcritical case, using the relation

$$\rho =\frac{P(x)}{4\pi }\widehat{r}^2,$$
(24)
we obtain the scaling law for the maximum density,

$$\rho _{\mathrm{max}}|pp^{}|^{2\gamma }.$$
(25)
In fact, for Eqs. (18)–(21), there is a discrete set of self-similar solutions with analyticity at the center and at the sonic point, such as the homogeneous free-fall, Larson-Penston, Hunter (a), (b), (c) and (d) solutions, and so on. Larson-Penston solution describes the coherent collapse of the cloud for $`t<0`$. This solution can be extended beyond $`t=0`$ to a late-time solution in which a finite mass collapses to the center and grows with time $`t`$. Hunter (a) solution for $`t<0`$ describes an “exploding collapse” in which the mass of the central collapsing region gets smaller and smaller, converging to zero as $`t`$ reaches $`0`$, being surrounded by an expanding envelope. The central density still blows up to a “singularity” at $`t=0`$. This solution can also be extended beyond $`t=0`$ to a late-time solution in which the whole cloud expands away; the singularity disappears and the central density keeps decreasing.

Maeda and Harada mh2001  found a unique relevant (unstable) mode for the Hunter (a) solution but no relevant mode for the Larson-Penston solution. It was concluded that the former is a critical solution while the latter is an attractor solution. In the former case, the eigenvalue $`\kappa `$ of the unique relevant mode was evaluated to $`\kappa 9.4637`$ by solving the eigenvalue problem. Other members of Hunter’s family have more than one relevant modes, which implies that these solutions are not critical solutions. The homogeneous free-fall solution suffers from kink instability op1988 . The value of the critical exponent $`\gamma `$ is $`0.10567`$. Since we have seen that the parameter for the family of initial data sets is given by the temperature $`T`$, the scaling laws in the present system are given by

$$m_{\mathrm{core}}|TT^{}|^\gamma $$
(26)
for the supercritical case, and

$$\rho _{\mathrm{max}}|TT^{}|^{2\gamma },$$
(27)
for the subcritical case. It should be noted that the parameter is not necessarily the temperature but may be the mass or the radius of the cloud if we fix the temperature.

## III numerical simulations

We have performed numerical simulations of the spherical collapse of an isothermal Newtonian gas. The basic equations for the evolution of the system are given by Eqs. (1)–(4). The code is based on the Lagrangian description of hydrodynamics and a finite difference scheme. We have chosen Lagrangian (comoving) radial coordinates, attached to spherical shells with infinitesimal thickness. The mass inside coordinate $`\overline{r}`$ is consequently constant throughout the simulation. We call such coordinates mass coordinates. Shells are initially positioned at constant intershell distance and their mass is computed from the desired density profile (we do not use equal mass shells). With these choices, we achieve high accuracy for the dynamics in the dense central region where the shells collapse during the simulations. The number of shells used in the simulation is $`10^4`$. One simulation with $`10^5`$ shells has been run to check the influence of numerical resolution. No effect has been detected. This is proof of the accuracy of the code. In this section, we adopt units in which $`G=k/m=1`$ and the total mass $`m_{\mathrm{tot}}`$ and initial radius $`r_{\mathrm{surf},\mathrm{i}}`$ of the cloud are unity.

We have computed the evolution of four types of initial density profiles. Details of the models for the initial density profiles are summarized in Table 1. The initial velocity profile is always set to be identically zero. This assumption seems the most relevant for astrophysical applications, but the same dynamics should develop for different initial velocity profiles. Model 2 is our reference model, in which the initial density profile is given by a quartic polynomial and both the density and the density gradient (i.e., pressure forces) are continuous at the cloud surface. We have computed collapse simulations for various values of the temperature. (Actually, we have computed for various sound speeds and calculated the temperature because the sound speed is rather dynamically meaningful.)

For each model, we find a critical temperature $`T^{}`$, which is of the same order of magnitude as the virial temperature estimated from the initial gravitational energy of the system. For temperatures just below $`T^{}`$, the cloud collapses and the central density goes to infinity in a finite time, producing a singularity. We call this collapse a supercritical collapse. For temperatures just above the critical temperature, the cloud begins to collapse but bounces afterward, and the density finally goes to zero everywhere. We call this incomplete “collapse” a subcritical collapse. The results for the critical temperture for each model are summarized in Table 1. In this table, we quote only the number of digits that are common between the critical values computed using the two different resolutions, $`10^4`$ and $`3\times 10^4`$.

In the critical and supercritical cases, we can estimate the time $`t_{\mathrm{sing}}`$ of singularity formation at the center in different ways. First we can just define it as the instant when the numerical density reaches a very high threshold (e.g., $`\rho 10^{15}`$ while initially of the order of one). Or we can assume the convergence to a self-similar behavior and use the relation $`\rho _\mathrm{c}(t_{\mathrm{sing}}t)^2`$ for the central density $`\rho _\mathrm{c}`$. Both approaches have been tried and give the same results. When the temperature is fine-tuned to the critical one, a case that we call the critical case for convenience, the simulations show that the collapse approaches the Hunter (a) solution for a while before it deviates. The finer the tuning of the temperature, the closer the approach and the later the deviation. Figure 1 shows density profiles in terms of the dimensionless quantity $`4\pi \rho (t_{\mathrm{sing}}t)^2`$ for a tuning of the temperature to the critical temperature $`|\delta T|/T^{}|TT^{}|/T^{}10^{15}`$, along with the theoretical Hunter (a) solution. In this case we estimate $`t_{\mathrm{sing}}`$ using the $`\rho _\mathrm{c}(t_{\mathrm{sing}}t)^2`$ relation. This estimation is local in time and it gives the value which $`t_{\mathrm{sing}}`$ would take if Hunter (a) solution was fully realized. We can see that the density profile converges to Hunter (a) solution. A discrepancy appears at large radii because the numerical configuration has a finite mass while the theoretical Hunter (a) solution has an infinite mass. This discrepancy moves to larger and larger value of $`r/(c_\mathrm{s}t)`$ as we approach $`t_{\mathrm{sing}}`$. Therefore, we associate the Hunter (a) solution to the critical case.

In the near-critical supercritical case, the collapse first approaches the critical solution. Figure 2 shows the evolution of the nondimensional central density parameter

$$Q\mathrm{ln}[4\pi G\rho _\mathrm{c}(t_{\mathrm{sing}}t)^2],$$
(28)
as a function of the central density $`\rho _\mathrm{c}`$ itself. Here again, we estimate $`t_{\mathrm{sing}}`$ locally. We compute $`t_{\mathrm{sing},n}`$ at each time step $`n`$, using the following relation

$$\rho _{\mathrm{c},n}(t_{\mathrm{sing},n}t_n)^2=\rho _{c,n1}(t_{\mathrm{sing},n}t_{n1})^2,$$
(29)
where $`t_n`$ and $`\rho _{\mathrm{c},n}`$ are the time and central density at the $`n`$-th step. We can then deduce $`Q`$. In this figure, the theoretical values for self-similar solutions, such as the homogeneous free-fall, Larson-Penston and Hunter (a) solutions, are also plotted. We find that the curve of the critical collapse approaches the value of the Hunter (a) (critical) solution. However, it afterward deviates and approaches the value of the Larson-Penston solution. This happens because the tuning of the temperature to the critical temperature is not perfect, and because the critical solution is unstable against such fluctuations. During and after the approach to the critical solution, the surrounding region expands away and leaves a finite-mass “core” in the central region. The density field in the central collapsed core at late times is well described by the Larson-Penston solution. This is the case early on for a supercritical collapse also, where the temperature is not fine-tuned to the critical temperature. This type of collapse also produce a finite (more massive) core and an expanding shell. We can see the convergence to Larson-Penston solution in Fig. 3. Therefore we identify the Larson-Penston solution as an attractor solution.

Figure 4 shows the velocity profile around the center for a supercritical collapse. We can see that a finite-size collapsing core is formed. Moreover the physical radius of the collapsing core is almost constant. Figure 5 shows the time evolution of the collapsing mass, for different (supercritical) values of the temperature. We can see that at all supercritical temperature, the collapsing mass tends to a constant and therefore the “core” mass is well defined at the time of the singularity. Figure 6 shows the relation between the core mass and the temperature fluctuation $`|TT^{}|`$. The numerical results reproduce the scaling law derived by the renormalization group analysis in the previous section. The numerical value of the critical exponent is $`0.11`$, which agrees well with the theoretical value $`0.10567`$ determined from the eigenvalue problem of the linear mode analysis. It should be noted that the power-law relation holds approximately for $`|\delta T|/T^{}0.01`$.

For the subcritical collapse, the central density reaches a maximum and then bounces away. Figure 7 shows the relation between the maximum density and the temperature excess $`TT^{}`$. It shows that the maximum density also follows a scaling law, as is expected from the renormalization group analysis. The measured value $`0.22`$ of the critical exponent agrees with the theoretical value $`0.21134`$. Again, it should be noted that the power-law relation holds approximately for $`\delta T/T^{}0.01`$.

Although we have concentrated on model 2, the above features have been checked for models 1 and 3. They are identical in all three models except for the specific value of the critical temperature. Model 4, the top-hat model, is special because the cloud is initially homogeneous. In this case, the central region of the cloud remains homogeneous and freely falls to a singularity if the condition $`m_{\mathrm{tot}}/r_{\mathrm{surf},\mathrm{i}}>\pi ^2c_\mathrm{s}^2/2`$ is satisfied ti1999 . If the condition is not satisfied, the collapse converges to Larson-Penston (attractor) solution.

## IV astrophysical implication

In this section, we study the mass spectrum of collapsed cores resulting from temperature fluctuations around the critical case, according to the treatment of Yokoyama yokoyama1998 . We build a toy model which shows how the critical gravitational collapse may affect the properties of the interstellar medium, and in particular the processes of fragmentation in molecular clouds and star formation. It should be emphasized that the purpose of this section is not to derive the whole shape of the observed initial mass function (IMF) but to show how a link can be established between an initial fluctuation field, and the resulting mass spectrum for the collapsed objects. Numerous additional physical processes would have to be introduced to derive an IMF.

We will use

$$M_{\mathrm{core}}(T)=KM_{\mathrm{tot}}\left(1\frac{T}{T^{}}\right)^\gamma ,$$
(30)
for the masses of cores formed by collapsing isothermal gas clouds with mass $`M_{\mathrm{tot}}`$ and a temperature which is slightly less than $`T^{}`$. $`K`$ is a constant and $`\gamma 0.11`$ is independent of initial density profile. This relation is theoretically valid for $`TT^{}T^{}`$. We can check on Fig. 6 that it is actually valid for $`|TT^{}|/T^{}0.01`$ or $`M_{\mathrm{core}}0.3M_{\mathrm{tot}}`$.

We assume a Gaussian probability distribution for the temperature fluctuations for simplicity,

$$P(T)=\frac{1}{\sqrt{2\pi }\sigma }\mathrm{exp}\left(\frac{(TT_\mathrm{m})^2}{2\sigma ^2}\right),$$
(31)
where $`T_\mathrm{m}`$ is the temperature with the maximum probability and $`\sigma `$ is the dispersion, which has the dimension of a temperature. For simplicity, we consider a collection of clouds of the same mass and size with a distribution of temperature obeying the above probability function. Equation (31) allows us to compute the fraction of clouds collapsing to cores

$$\beta =_{T^{}(1K^{1/\gamma })}^T^{}P(T)𝑑T.$$
(32)
The lower integration limit reflects the fact that a core with a mass heavier than the initial gas cloud cannot be formed.

We can now determine the mass sprectrum of cores formed from a collection of gas clouds with masses $`M_{\mathrm{tot}}`$. We define the mass spectrum as the number $`dN`$ of cores per logarithmic mass bin, normalized by

$$_{\mathrm{}}^{\mathrm{ln}M_{\mathrm{tot}}}\frac{dN}{d(\mathrm{ln}M_{\mathrm{core}})}d(\mathrm{ln}M_{\mathrm{core}})=1.$$
(33)
This mass function is given by

$`{\displaystyle \frac{dN}{d(\mathrm{ln}M_{\mathrm{core}})}}`$ $`=`$ $`\beta ^1P(T(M_{\mathrm{core}})){\displaystyle \frac{dT}{d(\mathrm{ln}M_{\mathrm{core}})}},`$ (34)
$`=`$ $`{\displaystyle \frac{T^{}}{\sqrt{2\pi }\beta \sigma \gamma }}\left({\displaystyle \frac{M_{\mathrm{core}}}{KM_{\mathrm{tot}}}}\right)^{1/\gamma }\mathrm{exp}\left({\displaystyle \frac{1}{2\sigma ^2}}\left[T^{}\left\{1\left({\displaystyle \frac{M_{\mathrm{core}}}{KM_{\mathrm{tot}}}}\right)^{1/\gamma }\right\}T_\mathrm{m}\right]^2\right).`$
From numerical simulations, we have obtained $`\gamma 0.11`$ and $`K1.2`$. There are two parameters $`pT^{}/T_\mathrm{m}`$ and $`q\sigma /T_\mathrm{m}`$ in the mass spectrum. The mass spectra are depicted for $`(p,q)=(1.0,0.1)`$ and $`(p,q)=(1.0,0.7)`$ in Fig. 8(a), while for $`(p,q)=(0.1,0.1)`$ and $`(p,q)=(0.1,0.7)`$ in Fig. 8(b). Core masses heavier than the total mass are cut off. The mass $`M_p`$ associated with the peak in the spectra is smaller than $`M_{\mathrm{core}}=M_{\mathrm{tot}}`$ for $`(p,q)=(1.0,0.1)`$ and $`(0.1,0.1)`$, while equal to $`M_{\mathrm{core}}=M_{\mathrm{tot}}`$ for $`(p,q)=(1.0,0.7)`$ and $`(0.1,0.7)`$ \[these regions are however not fully relevant because of the restrictions to be applied to Eq. (30\]. $`M_\mathrm{p}`$ and the average mass $`M_{\mathrm{av}}`$ of the formed core defined as

$$M_{\mathrm{av}}\frac{_0^{M_{\mathrm{tot}}}M_{\mathrm{core}}\frac{dN}{dM_{\mathrm{core}}}𝑑M_{\mathrm{core}}}{_0^{M_{\mathrm{tot}}}\frac{dN}{dM_{\mathrm{core}}}𝑑M_{\mathrm{core}}},$$
(35)
are summarized in Table 2. Although $`M_\mathrm{p}`$ can be just equal to $`M_{\mathrm{tot}}`$, the average mass $`M_{\mathrm{av}}`$ is always smaller.

The main feature of the mass spectrum is a steep slope in the low-mass tail, due to the small value of the critical exponent $`\gamma `$. This fact has several implications. First, the cores are created with nearly the mass of the initial gas cloud. In the formalism originally introduced by Press and Schechter ps1974 , it is assumed that overdense regions with a single mass scale, if they collapse, result in bound objects with the same mass. However, our result implies that the mass of the collapsed object is not the same as the mass of the initial overdense region, but may be rather smaller if the collapse is near-critical. Second, the critical collapse alone cannot account for the shape of the observed IMF, which shows a negative exponent on most of the mass spectum while we find a large positive exponent. Other phenomenon are responsible for the negative exponent.

## V summary

We have analyzed the critical behavior in the collapse of a Newtonian isothermal gas by the renormalization group approach and in numerical simulations. The renormalization group analysis in the Newtonian case is similar to the general relativity case. We have observed the critical phenomenon in the numerical simulation, which can be compared to those observed during black hole formation in the collapse of radiation fluid in general relativity. The scaling law for the formed black hole mass with critical exponent $`0.36`$ in general relativity is replaced by a scaling law for the collapsed core mass with critical exponent $`0.11`$. The critical solution is a continuous self-similar solution, specifically the first member of Hunter’s family of solutions.

For the supercritical collapse, the collapse near the center converges (early for low temperature, late for almost critical temperature) to another continuous self-similar solution, the Larson-Penston solution. Generic collapse converges to this attractor solution, which suggests the universality of self-similarity in gravitational collapse. This phenomenon is an example for the validity of the self-similarity hypothesis in gravitational collapse.

Some implications of Newtonian critical phenomena to astrophysics have been discussed. The critical behavior alone does not determine the IMF, clumps mass spectrum in molecular clouds and so on. However it suggest some limitations to Press-Schechter’s model. The isothermal gas is a good approximation for a gas system with a short cooling time, while it will not be valid for optically thick clouds, where adiabatic gas model are more relevant. Consequently, critical phenomena should be investigated for more general equations of state. Research along this line is now under investigation maeda2002 .

Newtonian gravity is a simple model for more complicated gravitational theories, such as general relativity. Newtonian critical behaviors in gravitational collapse provide a good laboratory for investigating various effects on the critical phenomena in general relativity, such as deviation from spherical symmetry, rotation, external environment, a variety of matter fields, and so on.

###### Acknowledgements.

We are grateful to T. Koike, B.J. Carr and R. Tavakol for helpful comments. TH and BS were supported by the Grant-in-Aid for Scientific Research (Nos. 05540 and 00273) from the Japanese Ministry of Education, Culture, Sports, Science and Technology.


format: text

text:
# 1 Introduction

## 1 Introduction

The existence of cosmic dark matter is required by a multitude of observations and arguments, such as the excessive peculiar velocities of galaxies within clusters of galaxies or gravitational arcs, indicating much deeper gravitational potentials within clusters than those inferred to be present from only the luminous matter content. Also an epoch of Big Bang nucleosynthesis predicts a fractional contribution of baryons to the critical density, $`\mathrm{\Omega }_b`$, significantly smaller than the total $`\mathrm{\Omega }`$ in form of clumpy matter inferred to exist from observations of large-scale galactic peculiar velocities (e.g., $`\mathrm{\Omega }0.3`$). Establishing the nature of the dark matter is one of the outstanding problems and challenges in cosmology. Virtually all proposed candidates require physics beyond the standard model of particle physics. A candidate particularly well-suited to the formation of the observed large-scale structure, is a weakly interacting particle with comparatively small velocities before the onset of structure formation (i.e. cold dark matter, herafter CDM).

Due to the small thermal velocity of CDM, fluctuations survive from the early universe on all scales. CDM is a bottom up scenario in which structures develop as small clumps collapse and undergo series of merging resulting in the hierarchical formation of massive dark matter halos. These halos are the hosts of baryonic systems, from galaxies to clusters, which cooled and condensed by dissipating energy. Cosmological N-body simulations describe well this non-linear complex series of mergers, accretion events and violent relaxation which are in general agreement with observations but are subject to some difficulties on galactic scales. An overview of these issues was recently presented by . It remains to be seen if the discrepancies on small scales may be solved by the particulars of baryonic physics (e.g. star formation, feedback, galactic black holes), or if indeed, the paradigm of cold dark matter is being challenged. In contrast, the excellent performance of cold dark matter (with inital adiabatic, scale-invariant density perturbations) on large scales seems still undiminished.

On the particle physics side, supersymmetric models are believed to provide the most promising approach to physics beyond the standard model. Supersymmetry can cure in principle several conceptual shortcomings of this model, although definitively compelling scenarios are still to come. Of interest for us here is the prediction of, as yet undetected, bosonic particles for each fermionic standard model particle, and vice versa, expected to lie in the approximate mass range of $`100GeV`$ \- $`1TeV`$ which are actively searched for at present and future colliders. In particular, the lightest supersymmetric particle (LSP), under the assumption of its stability (R-parity conservation) or at least quasi stability on the cosmological time scale, may serve as the dark matter particle as long as it is neutral (i.e. a neutralino). Moreover, its interactions are believed such that it would behave as CDM, and its relic abundance may be naturally of order $`\mathrm{\Omega }1`$ within the minimal supersymmetric extension of the standard model (MSSM), for review see . This is of course a very exciting possibility.

There is a large number of ongoing underground experiments attempting to detect neutralinos within the solar system by ”direct” scattering on nuclei. An alternative approach for the detection of neutralinos is via their occasional annihilation in dark matter halos, and observation of the resultant gamma-rays (and/or neutrinos). However, when it comes to quantitative predictions in relation to experimental detection, one has to tackle with astrophysical uncertainties such as the halo modelling of the astrophysical object under consideration, as well as with particle physics uncertainties related to our ignorance of the physics underlying supersymmetry breaking. The possible detection of such a gamma-ray signal from the spiral galaxy M31 (Andromeda), at a distance of 700kpc, is the subject of the present paper. We will deal mainly with two complementary issues: (i) we give different models of the dark halo component in addition to the disk and bulge components previously studied in the literature; (ii) assuming this dark matter is accounted for by the lightest neutralino of the MSSM, we analyse the gamma-ray fluxes that can originate from the pair annihilation of these particles, in terms of model assumptions such as minimal supergravity.

The paper is organized as follows: section 2 is devoted to M31 halo modelling. In section 3 we first recall briefly the basic ingredients and tools we use for the study of the supersymmetric dark matter signature and then present the main predictions. Section 4 explores the discovery potential of the CELESTE experiment. Section 5 is devoted to a quantitative discussion of possible signal enhancement due to clumpiness and black hole presence in M31.

## 2 Modeling the neutralino halo around M31

### 2.1 Rotation curve of M31 without dark halo

The late-type Sb spiral galaxy M31 lies at a distance of 700 kpc . The visible part mostly consists of a bulge and a disk. The present analysis is based upon the investigation of the neutral hydrogen content of M31, and the model-independent derivation of the velocity field, performed by Braun . The rotation curve, after correction for the ellipticity which the spiral exhibits in the inner 5 kpc, is well fitted with two mass components which are both traced by optical observations: (i) A bulge – with total mass $`7.8\pm 0.5\times 10^{10}`$ M – and assuming a mass-to-light ratio of $`\mathrm{{\rm Y}}_{bulge}=6.5\pm 0.4\mathrm{{\rm Y}}_{B,}`$ where $`\mathrm{{\rm Y}}_{B,}`$ is the mass-to-light ratio for the Sun. This value may be compared to the limit of $`3.7\mathrm{{\rm Y}}_{bulge}5.7`$ that has been derived from synthetic stellar models  assuming the same bulge color and a population age lying in the range from 9.5 Gyr to 15.5 Gyr; (ii) The disk with a mass of $`1.22\pm 0.05\times 10^{11}`$ M in the inner 28 kpc.

These contributions of the disk and the bulge components to the rotation curve are plotted in Fig. 1. Note that this figure reproduces Fig. 8b of . Braun concludes that no dark halo is necessary to account for the velocity field inside M31. This result relies, nevertheless, on the crucial assumption that the mass-to-light ratio of the disk is $`\mathrm{{\rm Y}}_{disk}=6.4\pm 0.4\mathrm{{\rm Y}}_{B,}`$. Such a large value does not agree with estimates based on the blue color of the disk and on synthetic spectra of young stellar populations which it contains . The mass-to-light ratio $`\mathrm{{\rm Y}}_{disk}`$ of a purely stellar component should actually not exceed $``$ $`3.8\mathrm{{\rm Y}}_{B,}`$. In addition, a disk as massive as that proposed by Braun should generally be unstable. We therefore feel that Braun has overestimated the importance of the disk so that a halo around M31 is indeed a viable possibility.

### 2.2 Rotation curve of M31 including a dark halo

Since a disk and a bulge may not be enough to model the rotation velocity of M31, we have assumed the presence of an additional mass component in terms of a spherical halo whose mass density profile is generically given by

$$\rho _\chi (r)=\rho _0\left(\frac{r_0}{r}\right)^\gamma \left\{\frac{r_0^\alpha +a^\alpha }{r^\alpha +a^\alpha }\right\}^ϵ.$$
(1)
A cored isothermal profile with core radius $`a`$ corresponds to $`\gamma =0`$, $`\alpha =2`$ and $`ϵ=1`$. A NFW profile  is obtained with $`\gamma =1`$, $`\alpha =1`$, and $`ϵ=2`$, whereas Moore’s distribution  is recovered if $`\gamma =ϵ=3/2`$ and $`\alpha =1`$.

For a given halo profile our mass models have two free parameters: the mass-to-light ratio of the disk and that of the bulge. While there is ample freedom on the relative importance of the disk and the halo, the observations set stringent constraints on the structure of the neutralino halo and actually favours a NFW  distribution. This is illustrated in Fig. 2 where the central singularity of the neutralino density has been set equal to $`\gamma =1/2`$ (left panel) and $`\gamma =3/2`$ (right panel). The contributions of the bulge and the disk correspond to the dotted and dashed curves, respectively, whereas the dashed-dotted lines denote the halo. The global rotation – solid – curves fail to match the observation points for medium range of R (a) and below 5 kpc (b). We have therefore disregarded these profiles in what follows and have concentrated on the NFW ($`\gamma =1`$) case. This is not surprising insofar as a $`1/r`$ spherical profile leads to the same rotation curve as a disk with constant surface mass density.

The neutralinos potentially concealed around M31 should annihilate and produce high-energy photons. The corresponding flux at Earth, $`n_\gamma `$ – per unit of time, surface, and solid angle – may be expressed as

$$\frac{dn_\gamma }{dtdSd\mathrm{\Omega }}=\frac{1}{4\pi }\frac{\sigma vN_\gamma }{m_\chi ^2}_{\mathrm{los}}\rho _\chi ^2𝑑s=\frac{1}{4\pi }\frac{\sigma vN_\gamma }{m_\chi ^2}𝒥(R)$$
(2)
where $`m_\chi `$ is the neutralino mass and $`\sigma vN_\gamma `$ denotes the thermally averaged annihilation rate yielding $`N_\gamma `$ gamma-rays in the final state. The astrophysical part of the expression consists in the integral $`𝒥`$ along the line of sight of the neutralino density squared $`\rho _\chi ^2`$. Assuming a spherical halo with radial extension $`R_{\mathrm{max}}`$ leads to

$$𝒥(R)=\mathrm{\hspace{0.33em}2}_0^{\sqrt{R_{\mathrm{max}}^2R^2}}\rho ^2\left(\sqrt{s^2+R^2}\right)𝑑s,$$
(3)
for a direction with impact parameter $`R`$ off the centre of M31. Because the density decreases steeply at large distances (see Eq. 1), our results are not sensitive to the actual value of the radial cut-off $`R_{\mathrm{max}}`$. The next step is the sum of the line of sight integral $`𝒥`$ over the solid angle subtended by the source

$$\mathrm{\Sigma }=𝒥(R)𝑑\mathrm{\Omega }.$$
(4)
We are interested in the number $`I_\gamma `$ of high-energy photons – collected per unit of time and surface – that originate from a circular region with angular radius $`\theta _{\mathrm{obs}}`$. The previous expression simplifies into

$$\mathrm{\Sigma }=\mathrm{\hspace{0.33em}2}\pi _0^{\theta _{\mathrm{obs}}}𝒥(R)\mathrm{sin}\theta d\theta ,$$
(5)
where $`R/D=\mathrm{tan}\theta \theta `$ and $`D700`$ kpc is the distance to M31. Integrating relation (2) over the source leads to the gamma-ray signal

$$I_\gamma =\left(3.18\times 10^{13}\text{photons}\text{cm}^2\text{s}^1\right)\left\{\frac{\sigma vN_\gamma }{10^{25}\text{cm}^3\text{s}^1}\right\}\left\{\frac{500\text{GeV}}{m_\chi }\right\}^2\mathrm{\Sigma }_{19},$$
(6)
where $`\mathrm{\Sigma }_{19}`$ denotes $`\mathrm{\Sigma }`$ expressed in units of $`10^{19}\text{GeV}^2\text{cm}^5`$.

In order to estimate the gamma-ray emission of an hypothetical neutralino halo around M31, we have considered an intermediate NFW halo with bulge and disk mass-to-light ratios of $`\mathrm{{\rm Y}}_{bulge}=4.2\mathrm{{\rm Y}}_{B,}`$ and $`\mathrm{{\rm Y}}_{disk}=4.2\mathrm{{\rm Y}}_{B,}`$ (left panel of Fig. 3). The parameter $`\mathrm{\Sigma }_{19}`$ of expression (6) reaches 1 in the inner 3.5 kpc which corresponds to $`5`$ mrad, field of view of CELESTE and 1.2 up to 28 kpc which is the M31 radius. A fiducial model for the neutralino has been assumed with $`m_\chi =500`$ GeV and annihilation cross section multiplied by mean velocity and average number of $`\gamma `$ in final state: $`\sigma vN_\gamma =\mathrm{\hspace{0.17em}10}^{25}\text{cm}^3\text{s}^1`$. The corresponding gamma-ray fluxes are respectively $`3.2\times 10^{13}`$ and $`3.8\times 10^{13}\text{photons}\text{cm}^2\text{s}^1`$ as presented in Table 1. In the right panel of Fig. 3, a more massive halo has been considered. The mass-to-light ratios of the bulge and of the disk are $`\mathrm{{\rm Y}}_{bulge}=3.5\mathrm{{\rm Y}}_{B,}`$ and $`\mathrm{{\rm Y}}_{disk}=2.5\mathrm{{\rm Y}}_{B,}`$, respectively. In this more favourable situation for detection, the gamma-ray emission is increased by a factor of 3 with respect to the previous case. Recent data at radial distance $`15r30`$ kpc, employing different data reduction analysis than that by Braun imply the presence of more mass at large radius than in Braun’s curve and a flatter mass distribution . Therefore one may well be close to the most favourable case.

## 3  Supersymmetric model predictions

In the present section we concentrate on more specific particle physics model predictions of the $`\gamma `$ \- ray fluxes. The aim is to implement, test and use a computational tool which allows to scan over various observables related to supersymmetric dark matter. This is achieved through an interfacing of two public codes, DarkSUSY  and SUSPECT , which we dub hereafter DSS (DarkSUSY-SUSPECT). Significant features on both particle physics and cosmology sides are thus combined within our approach. For instance, we will study the effects of universality versus non-universality of the soft supersymmetry breaking parameters and/or the (necessary) requirement of radiative electroweak symmetry breaking, both on the cosmological relic density of the neutralino LSP as well as on the gamma-fluxes from LSP annihilation in the halo of M31 as modelled in the previous section.

### 3.1 MSSM parameterization

We will focus mainly on the minimal supergravity scenario (mSUGRA) , where the soft breaking of supersymmetry occurs in a hidden sector, which communicates with the visible sector only via gravitational interactions, and translates to the observable particle physics sector in the form of soft masses and trilinear couplings among the scalar fields as well as mass terms in the gaugino sector . We make the usual simplifying assumption of common universal values of these parameters at the grand unified theory (GUT) energy scale $`M_{GUT}`$, i.e. $`m_{scalars}(M_{GUT})m_0,M_{gauginos}(M_{GUT})m_{1/2},A_{trilinear}(M_{GUT})A_0`$. Thus, starting from the common scalar soft supersymmetry breaking (SSB) mass

$$m_{\stackrel{~}{Q}}^2=m_{\stackrel{~}{U}}^2=m_{\stackrel{~}{D}}^2=m_{\stackrel{~}{L}}^2=m_{\stackrel{~}{E}}^2=m_{\stackrel{~}{H_1}}^2=m_{\stackrel{~}{H_2}}^2m_0^2,$$

the common SSB gaugino mass

$$M_1=M_2=M_3m_{1/2},$$

the common SSB trilinear coupling

$$A_t=A_b=A_\tau A_0$$

and unified gauge couplings $`\alpha _1=\alpha _2=\alpha _3\alpha _{GUT}`$, all taken at the GUT scale, the relevant low energy quantities are obtained from the renormalization group evolution of these parameters from $`M_{GUT}`$ down to a scale of the order of the electroweak scale. At this scale electroweak symmetry breaking is required through the minimization equations ,

$`{\displaystyle \frac{1}{2}}M_Z^2`$ $`=`$ $`{\displaystyle \frac{\overline{m}_1^2\mathrm{tan}^2(\beta )\overline{m}_2^2}{\mathrm{tan}^2(\beta )1}}`$ (7)
$`\mathrm{sin}2\beta `$ $`=`$ $`{\displaystyle \frac{2B\mu }{\overline{m}_1^2+\overline{m}_2^2}}`$ (8)
where $`\overline{m}_i^2m_{H_i}^2+\mu ^2+\text{radiative corrections}`$, $`\mu `$ is the supersymmetric mixing parameter of the two higgs doublets superfields and $`B`$ is the corresponding SSB parameter in the higgs potential. For a given $`\mathrm{tan}\beta `$ at the electroweak scale (defined as the ratio of the two Higgs vacuum expectation values $`<H_2^0>/<H_1^0>`$), one determines from the above equations the $`\mu `$ parameter (up to a sign ambiguity) and the B parameter consistent with the physical value of the $`Z`$ boson mass $`M_Z`$. In SUSPECT  the procedure is carried out including radiative corrections to the above EWSB conditions, renormalization group evolution to 1-loop order for the soft parameters, and to 2-loop order for gauge and Yukawa couplings (including threshold corrections from the supersymmetric spectrum). The full MSSM mass spectrum and couplings are computed and fed to DarkSUSY .

Various phenomenological constraints are taken into account (e.g. consistency with top, bottom and $`\tau `$ masses, present experimental limits on superpartner and Higgs masses, limits from bottom decay $`bs\gamma `$, no charged LSP, …), some of which are implemented in SUSPECT and others in DarkSUSY222We have also disactivated the $`g_\mu 2`$ constraint. Clearly a more refined treatment should require consistency with the standard model predictions . Furthermore, absence of charge and/or color breaking minima is not checked within our (exploratory) approach  .

Apart from mSUGRA with independent parameters $`m_0,m_{1/2},A_0,\mathrm{tan}\beta ,sign(\mu )`$, we will also consider an alternative and less constrained case study with 7 free parameters at the electroweak scale, namely $`\mu ,M_2,m_A,\mathrm{tan}\beta ,m_{\stackrel{~}{q}},A_b,A_t`$. Here $`M_2`$ is the soft wino mass (we still assume the GUT relation between $`M_2`$ and soft bino mass $`M_1`$), $`m_A`$ is the physical CP-odd higgs mass, $`m_{\stackrel{~}{q}}`$ the common value of all soft scalar masses at the electroweak scale. Furthermore, we do not require the EWSB constraints Eq.(7,8) in this context. It should be understood, however, that such configurations are considered here only as a test case of the sensitivity of gamma fluxes and neutralino relic density to large departure from more physically motivated parameterizations such as mSUGRA. We will refer to these configurations as low energy universality (LEU).
Finally, we will also consider configurations motivated by the focus point behaviour within mSUGRA . This allows in principle large $`m_0`$ values (in the TeV range) without having to fine-tune the parameters in Eq.(7). This comes about due to a peculiar behaviour in the running of $`m_{H_2}^2`$ from the GUT scale to the electroweak scale, which renders this soft mass fairly insensitive to the supersymmetry parameters provided that $`\mathrm{tan}\beta `$ is moderate or large. The implication for neutralino dark matter can be important as the LSP would acquire a non negligible higgsino component (in contrast to the generic mSUGRA almost purely bino prediction). This would affect the relic density as well as the gamma fluxes from LSP annihilation into W or Z pairs, when kinematically allowed, are no more suppressed as compared to the fermion-anti-fermion channels. Nonetheless, the “naturalness” of the focus point is actually moderated by a high sensitivity to the top quark mass in relic density calculations . We will thus adopt hereafter a qualitative standpoint where $`m_0`$ is allowed to be in the TeV range, disregarding fine-tuning issues. We dub these configurations FP-mSUGRA.

### 3.2 Relic density

In addition to the phenomenological or theoretical constraints mentioned in the previous section, one should impose conservatively  that the LSP relic density be in the cosmologically favoured region 0.1

<
[-0.07cm]
Ωχh2

<
[-0.07cm]
0.30.1

<
[-0.07cm]
subscriptΩ𝜒superscript2

<
[-0.07cm]
0.30.1\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\
[-0.07cm] $\sim$}}~{}\Omega_{\chi}h^{2}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\
[-0.07cm] $\sim$}}~{}0.3. However, we will also consider in the region 0.025

<
[-0.07cm]
Ωχh2

<
[-0.07cm]
0.10.025

<
[-0.07cm]
subscriptΩ𝜒superscript2

<
[-0.07cm]
0.10.025\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\
[-0.07cm] $\sim$}}~{}\Omega_{\chi}h^{2}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\
[-0.07cm] $\sim$}}~{}0.1 where we will apply a renormalization procedure on the expected fluxes.

Over the last decade the neutralino LSP relic density has been extensively studied in the literature (for reviews cf. ). For simplicity we focus mainly on those regions of the parameter space where co-annihilation effects  are unimportant in the estimate of $`\mathrm{\Omega }_\chi h^2`$. For large $`\mathrm{tan}\beta `$ values, Higgs resonances can enhance drastically the LSP annihilation leading to very small $`\mathrm{\Omega }_\chi h^2`$. Such effects are taken into account in our study. The actual bottom and top quark masses become a key issue in this case . We come back to this point later on when discussing benchmarks points.

When far from co-annihilation and resonance regions, the relic density can be estimated through :

$$\mathrm{\Omega }_\chi h^2\frac{1.07\times 10^9x_f}{g_{}^{1/2}M_{pl}(GeV)(a+\frac{b}{2x_f})}$$
(9)
where $`x_fm_\chi /T_f25`$ provides the approximate freeze-out temperature, and $`g_{}`$ counts the massless degrees of freedom at the typical temperature. Here $`a`$ and $`b`$ are defined as usual through the Taylor expansion in relative velocity of the thermally averaged annihilation cross-section times the neutralino velocity where $`<\sigma v>`$  is:

$$<\sigma v>=a+b\frac{T}{m_\chi }+O(\frac{T^2}{m_\chi ^2})$$
(10)
Note that the relation (10) is no more a good approximation in cases where resonance effects on the annihilation rate become important .

### 3.3 mSUGRA versus general MSSM

In this section we compare the results of simulations for MSSM (DarkSUSY) with mSUGRA models (DSS).

#### 3.3.1 Benchmark points comparisons

The SUSY benchmark models have been proposed by  to provide a common way of comparing the SUSY discovery potential of the future accelerators such as LHC or Linear Colliders. The thirteen SUSY scenarios correspond to thirteen configurations of the five mSUGRA parameters with the trilinear coupling parameter $`A_0`$ set to 0. The models fulfill the conditions imposed by LEP measurements, the $`g_\mu 2`$ result (which in our case will be not fulfilled), and the relic density constraint: $`0.1<\mathrm{\Omega }_\chi h^2<0.3`$.

The aim of the present section is to derive the gamma ray fluxes for some of these benchmark models with our current MC simulation programs: DarkSUSY  and SUSPECT , which were described previously. The value of $`\mathrm{\Omega }_\chi h^2`$ is calculated in the DarkSUSY part where only $`\chi ^+`$$`\chi _2^0`$ co-annihilations and annihilations related to the Higgs sector channels are included. As a consequence of these software limitations, the physically accessible SUSY domain in this study corresponds to models B,C,G,I and L of . The simultaneous use of the SUSPECT and DarkSUSY package allows to perform RGE evolution from the GUT scale to EWSB scale. Table 2 presents the initial GUT scale parameter values as well as other input parameters to DarkSUSY, calculated by SUSPECT ($`\mu `$, $`M_1`$, $`M_2`$, $`M_3`$, $`m_{\stackrel{~}{q}}`$, $`m_A`$). The GUT scale and EWSB scale values are also quoted as well as neutralino mass and gaugino fraction, $`R_g`$.

Table 3 presents the resulting relic density values with different theoretical assumptions corresponding either to mSUGRA or to LEU as described in section 3.1 whereas Table 4 gives results for $`\gamma `$-ray flux, $`\mathrm{\Phi }_\gamma `$ respectively. In order to be free from the Galactic Centre halo modelling, a renormalisation factor on the astrophysical part of the $`\gamma `$ flux computation has been applied. This allows a meaningful comparison between our results and those of . In case of comparison between mSUGRA and LEU, one has to choose some common features of the two configurations. We do this by choosing for the low energy universal masses in the LEU case, some of the representative values for the soft masses (and the $`\mu `$ parameter) at the electroweak scale, as given by the mSUGRA analysis. In the case of the soft scalar masses the choice is not unique. The numbers given in Tables 3 and 4 illustrate typical sensitivities.

As far as benchmark points with large $`\mathrm{tan}\beta `$ are concerned, the relic density is very sensitive to the input bottom mass since the latter controls the position of the s-channel CP-odd Higgs exchange pole .

The effect of applying radiative corrections has been found to be below 25% on flux predictions. The results obtained here are also compared to those in , where a different mSUGRA MC has been used and a more complete set of co-annihilation channels has been included. The observed differences on flux predictions of at most $``$ 25$`\%`$ between our results and those in  may well be explained by different treatments of the fragmentation processes with $`\pi ^0`$ in the final state.

The results given in Table 4 on the $`\gamma `$-ray flux from the Galactic Center are given for NFW  halo density parameterization, determined above a threshold of $`1`$ GeV within a solid angle of $`10^3`$ sr. For comparison, the flux from M31 has been also evaluated for the same energy threshold and acceptance value, assuming a NFW profile and model 3 in Table 1.

#### 3.3.2 Predictions from “wild scan” simulations

We have also performed a “wild scan” by having less restrictive conditions on the supersymmetric parameters (referred to as ”LEU” - see above). To achieve this, three thousand models have been simulated in the 0.025

<
[-0.07cm]
Ωχh2

<
[-0.07cm]
0.30.025

<
[-0.07cm]
subscriptΩ𝜒superscript2

<
[-0.07cm]
0.30.025\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\
[-0.07cm] $\sim$}}~{}\Omega_{\chi}h^{2}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\
[-0.07cm] $\sim$}}~{}0.3 region. The value of the lower limit on $`\mathrm{\Omega }_\chi h^2`$ equal to $`0.025`$ has been used with a renormalization procedure for flux estimation described below (for Ωχh2

<
[-0.07cm]
0.1subscriptΩ𝜒superscript2

<
[-0.07cm]
0.1\Omega_{\chi}h^{2}\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\
[-0.07cm] $\sim$}}~{}0.1). The low values of $`\mathrm{\Omega }_\chi h^2`$ may indicate that the main component of the dark matter is not the SUSY LSP and other contributions should be considered. The conditions of the performed simulations were for:

LEU 7 parameters:

$`10.<|\mu |<5000.`$

$`10.<|M_2|<1600.`$

$`10.<m_A<1000`$.

$`1.001<tan(\beta )<60.`$

$`50.<m_{\stackrel{~}{q}}<1000.`$

$`3.<A_t/m_{\stackrel{~}{q}}<3.`$

$`3.<A_b/m_{\stackrel{~}{q}}<3.`$

mSUGRA 5 parameters:

sign($`\mu `$) not constraint

$`50.<m_0<3000.`$

$`50.<m_{1/2}<1600.`$

$`0.1<|A_0|<2000.`$

$`3.<tan(\beta )<60.`$

Only co-annihilation with $`\chi ^+`$ or $`\chi _2^0`$ has been considered which is relevant when the higgsino component of neutralino becomes substantial. All results for the integrated gamma fluxes from M31 as a function of the $`\chi _1^0`$ mass were obtained for a NFW profile and model 3 in Table 1, and for $`\gamma `$-ray energy threshold of 30 GeV. The decay channels with $`\pi ^0`$ in the final state subsequently decaying into $`\gamma `$s, such as bb̄, $`W^+W^{}`$, $`Z^0Z^0`$, HH, …, were provided by the PYTHIA 6.1 program which is included in data format of the DarkSUSY package.

Figure 4 shows the integrated flux as a function of the neutralino relic density $`\mathrm{\Omega }_\chi h^2`$. The shape of the envelope corresponding to maximal fluxes can be understood qualitatively by using Eq. (2), Eq. (9), and Eq.(10) and maximizing $`<\sigma v>`$ in the parameter space. For instance in the $`b\overline{b}`$ channel, $`<\sigma v>`$ reaches a maximum when the lightest sbottom is degenerate with the LSP leading to $`<\sigma v>_{max}m_\chi ^2`$, that is $`\mathrm{\Omega }_\chi ^{min}m_\chi ^2`$. Combining this with the appropriate dependence of $`N_\gamma `$ on $`m_\chi `$ which can be obtained from the energy distribution

$$\frac{dN_\gamma }{dE}=a\frac{e^{b\frac{E}{m_\chi }}}{m_\chi (\frac{E}{m_\chi })^{1.5}}$$

by integrating over the full spectrum,

$$_{E_{min}}^{\mathrm{}}\frac{dN_\gamma }{dE}𝑑E=2a\{\sqrt{\frac{m_\chi }{E_{min}}}e^{b\frac{E_{min}}{m_\chi }}+(Erf[\sqrt{\frac{bE_{min}}{m_\chi }}]1)\sqrt{b\pi }\}$$
(11)
where $`E_{min}`$ is an acceptance energy cut, one recovers a behaviour similar to the envelope shown in Fig 4.

A renormalisation procedure has been used for low $`\mathrm{\Omega }_\chi `$$`h^2`$ values where the dark matter halo should not be dominated by the neutralinos. For values of $`\mathrm{\Omega }_\chi `$$`h^2`$ below 0.1 the predicted flux is renormalized by the ratio $`(\mathrm{\Omega }_\chi `$$`h^2/0.1)^2`$ to account for neutralinos only contributing partially to the dark halo. As a result of this procedure, the obtained fluxes went down by a factor of 10 at low values of $`\mathrm{\Omega }_\chi `$$`h^2`$. The effect of the renormalization procedure described above is also shown in Fig. 4.

Figure 5 presents the integrated gamma flux as a function of $`m_\chi `$ for the LEU scheme. The expected fluxes for M31 are of the order of $`10^{13}`$ $`\gamma cm^2s^1`$ and show only a weak dependence on the gaugino fraction of the neutralino. Comparing to mSUGRA predictions shown in Fig. 6, clear differences are observed:
\- the unification at GUT scale and electroweak symmetry breaking constraints lead to the exclusion of an important number of models in mSUGRA models, which are present in the LEU scheme,
\- the mSUGRA model reduces the number of viable Higgsino like models when compared to the LEU model, and
\- mSUGRA models generally show factor $`5`$ lower fluxes at low values of $`m_\chi `$ as compared to LEU models.
The benchmark model points are also presented.

Figure 7 presents mSUGRA results in case of the so-called “Focus Point” inspired scenario where $`m_0`$ values may vary up to 3 TeV. In this case there are more viable models with a neutralino of comparatively large Higgsino content than in mSUGRA and the mean fluxes reach values predicted by LEU.

The CP-odd Higgs mass resonance effect (i.e. $`m_A2m_{\chi ^0}`$) on the annihilation cross section and predicted fluxes has also been examined. Figure 8 presents an increase of the integrated flux at the $`m_A`$ pole by at least an order of magnitude compared to the region far from the resonance. This effect reflects increase of the predicted fluxes for large $`tan\beta `$ shown in Fig. 9, as the CP-odd Higgs s-channel dominates at large $`tan\beta `$ in the mSUGRA frame.

## 4 Detection of Dark Matter in M31 with
CELESTE

### 4.1 The CELESTE experiment

The “Cherenkov Low Energy Sampling and Timing Experiment” (CELESTE) is located on the site of Thémis in the French Pyrénées. It detects $`\gamma `$-rays in the high energy domain and was designed to fill the energy gap between satellites and imaging atmospheric Cherenkov telescopes (IACTs). While the former have been so far limited to $``$$`10`$GeV because of a too small effective detection area, the energy domain covered by the latter starts at $``$$`200`$GeV due to the noise arising from the night sky background. The energy window around $`50`$GeV has been opened by CELESTE during the winter 1999-2000 by the detection of the Crab nebula and its flux measurement at $`60`$GeV . Nearly at the same time, CELESTE detected the Active Galactic Nuclei (AGN) Mkn 421 during a series of flares. It was the first sub-$`100`$GeV detection of this nearby blazar and the emission observed by CELESTE showed a clear correlation with the flux recorded above $`250`$GeV by the CAT telescope, operating on the same site . Before these results, the sub-$`100`$GeV region was the only part of the electromagnetic spectrum which remained unexplored. The access to this region is of great importance for a number of key questions in high-energy gamma-ray astrophysics. For instance, the high-energy part of some AGNs can dominate the entire spectrum, giving new insight into their structure as well as the particle acceleration and cooling processes occuring in their vicinity. Pulsars could be also very promising sources : since their $`\gamma `$-ray spectrum is believed to decrease sharply at a few tens of GeV, the observation of the predicted spectral cutoff should strongly constrain the models describing the particle acceleration in their magnetosphere. Concerning galactic cosmic-rays, the current belief of their acceleration in shell-like supernovæ remnants could be also tested in a few cases by the observation of high-energy $`\gamma `$-rays. Apart from astrophysical issues, we propose here to use CELESTE for the study of dark matter and searches for the indirect signatures of neutralinos, namely their gamma-ray signal above $`50`$GeV coming from neutralino annihilations in M31. At the latitude of Thémis ($`42^{}`$ North), M31 transits very close to Zenith.

The CELESTE experimental setup is fully described in . This detector records the Cherenkov light emitted by the secondary particles produced during the development of the cosmic-ray atmospheric showers. In order to reduce the night-sky noise and to achieve a low energy threshold, CELESTE uses $`40`$ heliostats ($`53`$ since January 2002) of the former solar plant in Thémis, with a total collection area of $``$$`2000`$m2. An efficient discrimination between $`\gamma `$ and hadron-induced showers allows to extract the $`\gamma `$-ray signal from the more abundant charged cosmic-ray background (mostly protons, helium nuclei and electrons). The event analysis is based on the differences between both types of showers regarding the homogeneity and the time dispersion of the Cherenkov light pool as sampled by the heliostats. After the detector trigger and all analysis cuts, the $`\gamma `$-ray acceptance at Zenith has been parameterized as a function of the energy  :

$$𝒜(E)=A_0\left[1\mathrm{exp}\left(\frac{EE_0}{E_c}\right)\right],$$
(12)
with $`A_0=1.37\times 10^4`$m2, $`E_0=34.7`$GeV and $`E_c=56.2`$GeV.

### 4.2 Simulation results

The flux predictions of mSUGRA models were folded with the acceptance function of the CELESTE detector. Figure 10 shows the expected count rates from M31 in units of $`\gamma `$/min. The most optimistic models predicts values of the order of $`10^3`$, which should be compared with the $`4.9\gamma `$/min which are observed by CELESTE from the Crab nebula (the standard candle for high-energy $`\gamma `$-ray astronomy) . Note that these results were obtained without any signal enhancement effect expected from clumpiness of the dark matter halo or possible black hole neutralino accretion. Since the sensitivity level of atmospheric Cherenkov detectors is proportional to $`1/\sqrt{T_{\mathrm{obs}}}`$, only sources with a flux of a few tenths of the Crab nebula flux can be detected within realistic observation times. Therefore the overall enhancement due to astrophysical effects has to exceed at least a factor of 100 for detection of a $`\gamma `$-ray signal from neutralino annihilations to become feasible. The impact of such effects will be addressed in the next section.

A possibility of the shape discrimination with the energy differential flux measurement in CELESTE has been also studied. The expected shape for $`\gamma `$-rays produced in $`\chi ^0\chi ^0`$ annihilations are described by a decreasing exponential function with the exponent strongly depending on $`m_\chi `$. This sharp exponential cut-off should be compared with standard astrophysical source power law expectations.

The shape discrimination relies on the assumption of the $`\gamma `$ detection from the considered source, in our case M31. The event selection and cuts on the kinematical variables of the $`\gamma `$ candidates allow to enhance the $`\gamma `$ contribution to the measured sample with respect to the charged cosmic-ray initiated showers. In order to reduce further the charged cosmic-ray background and to optimize the signal to background ratio, various procedures such as a sliding window type method can be considered. This analysis will be addressed in details elsewhere.

Figure 11 presents predicted energy differential distributions for various $`m_\chi `$ ranges after the renormalization of the total integrated fluxes. Here an almost scaling variable equal to $`E_\gamma /m_\chi `$ was used as proposed by . No background contribution has been taken into account here. These results suggest an increase in sensitivity with $`m_\chi `$. However, due to the CELESTE acceptance varying with energy, a realistic discrimination measurements can only be considered for neutralino masses above 200 GeV. It is interesting to note that only models with substantial co-annihilation and large $`tan\beta `$ can lead to acceptable dark matter densities for $`\mathrm{\Omega }_\chi `$$`h^2`$ above 200 GeV.

## 5 Impact of Astrophysical parameters

### 5.1 Halo Clumpiness

Early simulations of the growth of CDM halos found smooth tri-axial structures, but higher resolution simulations revealed that a significant fraction of the mass may lie in substructures (sub-halos) orbiting in virialized dark matter halos . Such substructures seem to extend down to the smallest mass scale, possibly well below that of the most massive globular clusters.

The CDM scenario, however, predicts a larger number of substructures in galactic halos than the observed number of dwarf satellites in the Milky Way . The cases of gas expulsion by supernova-driven winds or the failure of gathering enough gas for stellar formation were considered in  as possible explanations. In this case there would be a large number of clumps with few stars, and some of the known high-velocity gas clouds may possibly trace them. It has been also shown recently that there is good agreement between the rotation curves of the largest satellites and those found in N-body simulations, provided star formation is less efficient in these satellites  . The possibility that the anomalous image flux ratios observed in several gravitational lenses could be explained by low mass satellites in the lensing galaxy was also explored . It was found that the required projected satellite mass fraction is of the same order as those predicted in CDM simulations, illustrating that there is general agreement between the underlying theory for structure formation and observed structures.

In order to estimate the effect of clumpiness we calculate first the contribution to the flux of a DRACO-type object with virial velocity $`v=10\mathrm{km}\mathrm{s}^1`$, and central mass $`M10^8M_{}`$ modeled by using a NFW profile with scale radius $`r_s=0.4\mathrm{kpc}`$ and density $`\rho _0=6\times 10^{24}\mathrm{g}\mathrm{cm}^3`$. This is obtained by calculating the sum of the line of sight integral $`\mathrm{\Sigma }`$, as defined in Sect. 2, over the field of view of CELESTE which we denote by $`\mathrm{\Sigma }^{10}`$ ($`\mathrm{\Sigma }^{10}=\mathrm{\Sigma }^{19}(3.5kpc)`$ defined in section 2) . It is expressed in units of $`10^{19}\mathrm{GeV}^2\mathrm{cm}^5`$ as before. Such a single DRACO-type clump orbiting around M31 would yield $`\mathrm{\Sigma }_{\mathrm{D}}^{}{}_{}{}^{10}0.03`$, to be compared with $`\mathrm{\Sigma }_{\mathrm{M31}}^{10}3`$ for M31 (cf. Sec. 2).

Now we consider the mass and the radial distribution of clumps orbiting around a halo. The sub-halo mass function is a power law close to $`dn(m)/dm(m/M_H)^{1.9}`$  where $`M_H`$ is the mass of the halo. The radial distribution of clumps in a halo, deduced from N-body simulations, behaves as $`n(r)\left[1+\left(r/r_c\right)^2\right]^{3/2}`$ . Combining these two equations we obtain the distribution of clumps having a mass $`m`$ at radial distance $`r`$:

$$n(r,m)(m/M_H)^{1.9}\left[1+\left(r/r_c\right)^2\right]^{3/2}v^{3.8}\left[1+\left(r/r_c\right)^2\right]^{3/2}$$
(13)
where $`v`$ is the virial velocity of the clump. This equation is normalized to the number of clumps of a given mass predicted by CDM simulations.

About 1000 clumps with velocities above $`v=10\mathrm{km}\mathrm{s}^1`$ are expected in a galactic halo having the mass of M31 . We truncate the clump mass distribution between $`1v10\mathrm{km}\mathrm{s}^1`$. We consider two spatial distributions of clumps: a compact distribution with core radius $`r_c16\mathrm{kpc}`$ and a more extended distribution with $`r_c30\mathrm{kpc}`$. Both are independent of the clump mass. In these profiles, around $`75\%`$ of the clumps are within the core radius.

The total integrated contribution of clumps with $`v1\mathrm{kms}^1`$ is $`\mathrm{\Sigma }2500\times \mathrm{\Sigma }_\mathrm{D}=77`$. Therefore the flux from M31 may be boosted by a clumpiness factor $`77/325`$. The clumps in the field of view of CELESTE yield respectively $`\mathrm{\Sigma }^{10}=20`$ for $`r_c16\mathrm{kpc}`$ and $`\mathrm{\Sigma }^{10}=10`$ for $`r_c30\mathrm{kpc}`$, that is a factor $`3.36.6`$ larger than the smooth profile of M31.

Globular clusters are known to continuously loose stars due to tidal forces from the halo gravitational potential. Mass is stripped from the outer regions of globular clusters and evaporates mainly at each passage through the disk. This is why it is generally argued that globular clusters do not have dark matter except perhaps in an inner bubble. The structural evolution of sub-halos was recently analyzed in , which shows that the tidal radius of distant satellites is well beyond their optical radius. This may also explain why the distant globular cluster Palomar 13 (at a distance of $`24.3\mathrm{kpc}`$) is well fitted by a NFW profile , i.e. this anomalous globular cluster could, in fact, be a Milky Way dark clump which has not been destroyed or modified by the galactic tidal field . The central mass-to-light ratio of the distant satellites Palomar 13 and DRACO (at a distance of 80 kpc) are, respectively, $`40(M/L)_{}`$ and $`80(M/L)_{}`$, whereas for nearby globular clusters it is approximately $`3(M/L)_{}`$. We assume that this is the signature of enhanced particle evaporation in the inner region of clumps. To account empirically for this effect we have divided by 2 the central density of the clumps at radial distance $`r=20\mathrm{kpc}`$, by 20 that of objects at $`r=10\mathrm{kpc}`$ and have ignored clumps at smaller radius. With this hypothesis the most massive objects within $`r=10\mathrm{kpc}`$ have a mass $`<5\times 10^6M_{}`$. In any case, a significant fraction of more massive objects is forbidden by the observations because it would dynamically heat the disk (see below).

With the hypothesis of evaporation, the integrated contribution of clumps is $`\mathrm{\Sigma }=1134`$ for the whole of M31 and to $`\mathrm{\Sigma }^{10}=0.752.2`$ in CELESTE’s field of view, indicating only very modest enhancement factors of the $`\gamma `$-ray signal.

If a large number of massive satellites orbit within the M31 halo, their passages through the disk causes disk heating  or may induce warping motions . This is confirmed by recent N-body CDM simulations of the Milky Way . In the case of a $`\mathrm{\Lambda }\mathrm{CDM}`$ cosmology, however, the clumps are less massive and are located at larger radial distance than in standard CDM, so only a small number get near the disk. In that case they do not heat the disk efficiently . We have already approximately accounted for the first effect by ignoring the most massive clumps with mass $`10^9M_{}`$ in the inner regions of galaxy. In order to account for an “extended halo distribution”, as in , we have considered the same clump distribution as before but with core radius $`100\mathrm{kpc}`$ and the same mass spectrum. With that distribution the clump flux in CELESTE’s field of view is $`\mathrm{\Sigma }^{10}=2.2`$ without evaporation, or $`\mathrm{\Sigma }^{10}=1.0`$ with the hypothesis of evaporation.

These $`\gamma `$ ray signal enhancement factors due to halo clumpiness should be compared to those determined by N-body simulations of Milky Way type objects. Whereas earlier results by one group  indicate a two order of magnitude effect, more recent simulations  find only a factor $``$ 2 enhancement. Another study performed on the clumpiness enhancement factor, in particular in case of M87, observation  predict values between 13 and 40. These calculations were extended for Milky Way or M31 type galaxy case leading to the values of the clumpiness enhancement factor of the order of 20 .

### 5.2 Black hole accretion

The central region of M31 contains a very compact object with a mass $`3.6\times 10^7M_{}`$ which is assumed to be a supermassive black hole (SBH) . The dynamical effect of a SBH on a distribution of nearby WIMPS was studied by  in the case of the Milky Way. The adiabatic growth of a SBH at the center of a halo produces a spike with density $`\rho r^\gamma `$, $`2.2\gamma 2.5`$ within a radius $`10\mathrm{pc}`$ for a SBH like that in M31. In such a spike the density reaches the annihilation density $`10^8\mathrm{M}_{}\mathrm{pc}^3`$ in a significant region surrounding the black hole producing a huge enhancement of the annihilation flux. This calculation assumes a very high accretion rate and a stable dynamical regime during an extremely long period. In fact the accretion rate for the assumed black hole adiabatic growth is much higher than the current estimates of present-day growth-rates . Several authors have discussed more realistic scenarios (i.e. ). In the CDM scenario, large halos grow through the merging of smaller building blocks. Galaxies as massive as the Milky Way or M31 have certainly experienced significant merging in the last Gyrs. There may be several SBH progenitors coming from the mergers, a SBH may be spiraling into the galactic potential. Motions of the SBH transfer energy to the particles lowering their density . Thus, even if at early times a massive BH was present in a dense environment resulting in a rapid adiabatic growth of the BH surrounded by a dense spike, there are a number of accidental events and perturbations during the life of the SBH which all tend to destroy the spike. Nevertheless a SBH will continuously accrete material. If it is coinciding with the dynamic center of a galaxy for a sufficiently long period it will built a central dark matter density cusp.

We performed N-body simulations assuming initial dark matter density profiles of $`\rho r^\alpha (1+r/r_c)^\beta `$ with $`0\alpha 1`$ and under the hypothesis of the existence of a massive central black hole - SBH. These simulations were performed with the aid of the public code Gadget . In each case, the density distribution reaches a stable limit with $`\rho r^{1.5}`$, independent of the initial profile, within dynamical time of $`10^6`$ yrs. An example of such simulations is shown in Fig. 12. This illustrates that a stable configuration which may be build in within a relatively short time.

The $`\gamma `$-ray flux enhancement due to the existence of such a dark matter spike as compared to a NFW profile is of the order of $`1+[ln(r_{max}/r_{min})1]r_{max}/R`$ where $`r_{min}`$ and $`r_{max}`$ are respectively the minimun and maximum radius of the profile in $`r^{1.5}`$ and $`R`$ is the radius of the initial region with profile in $`r^1`$. Taking $`R=3.5\mathrm{kpc}`$, $`20r_{max}100\mathrm{pc}`$, and $`r_{min}=10^5\mathrm{pc}`$ which is $`11`$ Schwarzschild radii, we get an only modest enhancement factor between 1.08 and 1.45. It is seen that the enhancement factor is small, nevertheless, the presence of the SBH in this complicated multicomponent object warrants a singular profile at small radius. This is also important with regard to evidence that the universal singular CDM profiles  are not seen in the dark halos of low surface brightness galaxies . In these objects, however, there are numerous suggestions for a hidden baryonic component  which may modify its inner profile complicating the interpretation of dwarf halos .

Finally in some cases the hypothesis of a central spike may be valid. This may happen if there is a density increase during a short period near the SBH, for example when a massive clump is falling in the SBH. As soon as the profile is steeper than $`r^{1.5}`$ the logarithmic term in the enhancement factor is replaced by a power of the inner radius. In the case of a singular isothermal (SIS) spike of $`1\mathrm{pc}`$ the total flux within $`3500\mathrm{pc}`$ is boosted by a factor 29. The impact of clumpiness and black holes are summarized in Tab. 5. The global enhancement factor on expected flux may vary between 5 and 100, depending on the clumpiness of the halo and CDM-accreting SBH. In particular, the large variation between $`\mathrm{\Sigma }`$ and $`\mathrm{\Sigma }^{10}`$ in case of an extended clump distribution (last line of Table 5) is due to the small angular size of the instrument compared to the compact distribution case (line 3 of Table 5).

## Conclusions

The putative dark matter content of M31 has been estimated and modelled from its observed rotation velocity profile, by adding a NFW dark halo to the observed disk and bulge component.

The study of the benchmark models in case of M31 allowed to qualify the DarkSUSY and SUSPECT MC programs used for gamma flux predictions in CELESTE experiment. Only models I and L (large $`tan\beta `$) could be considered as favourable for the detection. This conclusion is confirmed by a further study of the so-called “wild scan” simulation in mSUGRA scheme. Various aspects such as dependence of the results on heavy quark masses or CP-odd higgs pole contribution have been investigated.

We have taken into account the atmospheric electron and CR proton backgrounds to the M31 signal and applied the CELESTE gamma–ray energy dependent acceptance factors. Assuming a smooth neutralino halo around M31 leads to photon rates at the telescope of order a few $`10^3`$ counts per minute. The strongest signal corresponds to supersymmetric configurations where the neutralino mass $`m_\chi `$ exceeds $``$ 200 GeV. As a realisitic observation would require at least a signal of order a tenth of that of the CRAB – 4.9 photons per minute – we estimate that a gamma–ray annihilation signal from M31 is beyond the reach of an atmospheric Cerenkov detector of the CELESTE generation if the neutralinos are smoothly distributed in that galaxy. On the other hand, if the halo is made of clumps with inner profiles à la Moore or if there is strong accretion going on around the central black hole, the expected signal may be enhanced by two orders of magnitude and could become detectable. We conclude therefore that a survey of M31 with CELESTE is worth being undertaken.

## Acknowledgements

We would like to thank Jean-Loic Kneur for his help in running the SUSPECT program and building the interface between SUSPECT and DarkSUSY. We are also very much indebted to Piero Ullio for the time he spent on useful discussions about the DarkSUSY program and the physics related to the studied subjects. We thank the French GDR “PCHE” groups for supporting this study.