text:
# On the different role of protons and neutrons in antinucleon annihilations on nuclei.

## 1 Introduction.

Recently data on $`\overline{n}p`$ annihilation in the range 40-400 MeV/c (for the laboratory $`\overline{n}`$ momentum $`k`$) have been produced by the Obelix Collaboration. We would like to compare these with other data on $`\overline{p}p`$ and $`\overline{N}`$-nucleus annihilation. In particular, we are interested in: $`\overline{n}`$ annihilation on nuclei from 12C to 207Pb in the range 180-280 MeV/c; $`\overline{p}p`$ annihilation from 30 to 180 MeV/c; $`\overline{p}`$D, $`\overline{p}^4`$He and $`\overline{p}^{20}`$Ne at small momenta, down to 45 MeV/c; $`\overline{p}`$-nucleus annihilation on intermediate nuclei at larger momenta (over 200 MeV/c, see for a recollection of these data); and on the ratio between $`\overline{p}p`$ and $`\overline{p}n`$ annihilations inside nuclear targets. We will start by comparing $`\overline{p}p`$ and $`\overline{n}p`$ data, and then we will try to correlate these with the low energy nuclear data.

One of our aims is simply to compare $`\overline{n}p`$ and $`\overline{p}p`$ annihilation cross sections. This comparison is a delicate operation, because it requires subtraction of Coulomb effects, and comparison of data coming from different experiments with antiproton and antineutron beams. Due to the difficulties in calculating a 100 % reliable flux normalization, it is preferable, as far as possible, to compare data from the same experiment.

A second goal of this work is to establish which is the degree of correlation between antineutron-nucleon and antinucleon-nucleus data. Indeed, as better discussed later, $`\overline{p}p`$ data below 600 MeV/c can be well fitted via energy independent optical potentials, while $`\overline{n}p`$ data show a nontrivial energy dependence that should be reflected in the energy dependence of nuclear annihilations.

Another interesting point is the role of the different isospin channels in the low energy nuclear shadowing. It has been demonstrated that a marked nuclear shadowing characterizes low energy $`\overline{N}`$-nucleus annihilation cross sections. For this phenomenon it is possible to imagine two classes of explainations. First, general quantum mechanical processes, as described in a later section. Mechanisms of this kind do not discriminate too much between different annihilation channels. On the other side we may imagine a role for peculiar properties of single nuclear species. Indeed, the data are rather incomplete and impossible to organize in a systematic way. In particular a different $`N/Z`$ composition, or a different $`N/Z`$ distribution near the nuclear surface can be relevant, since most of the theoretical models (for reviews see e.g. ) establish peculiar properties for the different isospin channels, properties which are normally derived either from G-parity (or C-parity) transformations of the better known nucleon-nucleon potentials, or from peculiar quark diagrams. So, to summarize this point, it is relevant to establish how much of the low energy nuclear shadowing of annihilation processes can be due to general mechanisms, and how much is due to the nuclear composition in terms of protons and neutrons, and to the details of the nuclear structure.

We limit ourselves to a phenomenological analysis of the data, with standard methods of nuclear physics, without entering the debate about the underlying hadronic structure mechanisms governing the annihilation process. Of course, it is not always possible to separate completely the two levels of the analysis, i.e. nuclear and hadronic. Wide reviews on the debate concerning the annihilation mechanisms can be found in the references listed above. Experimental (scattering and annihilation) data up to 1994 are recollected in , data on antiprotonic atoms up to 1988 in , and a very recent account of the many problems related with the antinucleon experimental techniques can be found in .

We still lack a well founded method that can be properly used at low energies to link the shadowing properties (that automatically arise from an optical potential treatment) with the details of the nuclear structure and internal motions (which are taken into account in Impulse Approximation inspired treatments). The main exception in this respect is Deuteron, which has been considered by some authors and where it was possible to relate clearly the overall antinucleon-nucleus processes with the antinucleon-nucleon interactions. With more complex nuclei, optical potential analyses have been carried on to study antiprotonic atoms. At momenta over 200 MeV/c the KMT method has been used in  and the Glauber method in  (and probably in other works). In the present work our hope is that some interesting information can be separately extracted from a PWIA and an optical potential analysis.

For the sake of brevity we will indicate total $`\overline{n}`$ and $`\overline{p}`$ annihilation cross sections with TNA and TPA respectively, specifying the target. We will indicate as $`F_A`$ the ratio between $`\overline{p}n`$ and $`\overline{p}p`$ annihilations in a given nucleus with atomic number $`A`$.

Concerning the choice of an optical potential for discussing nucleon-antinucleon annihilation, it is well known that many quite different optical potentials fitted successfully pre-Obelix nucleon-antinucleon data and they can probably fit the data under discussion too, since these data do not represent such a strong contraint. Here and in previous papers we have relied on very simple optical potentials, with Woods-Saxon shape, the same parameters for each spin channel and no energy dependence in the momentum range 0-600 MeV/c (more indications are given in section 3). Although such potentials cannot reproduce the full phenomenology of nucleon-antinucleon interactions, they seem sufficient for fitting the available $`\overline{p}p`$ low energy data. With the $`\overline{n}p`$ annihilation data, as described in section 3, a higher level of sophistication is perhaps necessary.

## 2 The role of Coulomb interactions. Apparent violation of the isospin invariance.

Comparison of annihilation data with different projectiles and targets is nonsense without subtracting Coulomb effects at momenta below 200 MeV/c. The traditional estimation has undercome modifications in the last years. Qualitatively the effect of the Coulomb forces is to focus the projectile wavefunction in the annihilation region, which is a spherical shell of thickness 0.5-1 fm at the target proton/nucleus surface.

In a previous work we gave analytical expressions for the correcting $`\overline{p}`$-nucleus enhancement factors. E.g., for the $`\overline{p}p`$ annihilation cross section the enhancement factor can be reproduced within 5 % in the range 30-400 MeV/c by the function $`1+0.0003\beta ^2`$. To calculate the enhancement factor we considered some completely different optical potentials. All of them included the electrostatic potential of a spherical charge distribution and fitted the available low energy data. Then we removed the electrostatic potential and calculated again the annihilation rates. The obtained Coulomb enhancement factor was sufficiently independent of the choice of the strong part of the optical potential, in the considered momentum range.

An observation is necessary concerning apparent “violations” of the isospin invariance. Isospin symmetry suggests that: (i) $`\overline{n}p`$ cross sections are equal to $`\overline{p}n`$ ones; (ii) $`\overline{p}p`$ and $`\overline{n}n`$ are equal apart from electromagnetic effects; (iii) the four cross sections that one can imagine are actually combinations of two. In principle this is undoubtful, but in practice it does not work, and an example can clarify the point. Normally, $`\overline{n}p`$ data come from collisions between free antineutrons and protons, while $`\overline{p}n`$ cross sections are extracted from deuteron targets, with two consequences:

1) The antiproton is attracted by the deuteron charge. The range of action of the Coulomb forces is much larger than the deuteron radius, with the result that the $`\overline{p}n`$ annihilation rate is almost as much “Coulomb distorted” as the $`\overline{p}p`$ one, so it should be expected to be larger (much larger at very low energies) than the $`\overline{n}p`$ one. This problem is not present at momenta $`>>`$ 100 MeV/c, where Coulomb focusing effects can be neglected.

2) As already noticed, below 60 MeV/c shadowing effects in nuclei are so strong that the $`\overline{p}`$D annihilation cross section is smaller than the corresponding $`\overline{p}p`$ one in the laboratory frame, or approximately equal in the center of mass frame (due to the different transformations relating laboratory with center of mass variables for the cases of proton and deuteron targets). In both cases $`\sigma _{\overline{p}D}`$ is much smaller than $`\sigma _{\overline{p}p}+\sigma _{\overline{n}p}`$.

The previous example should clarify that, although it is very likely that one would find $`\sigma _{\overline{n}p}`$ $``$ $`\sigma _{\overline{p}n}`$ in an experiment on $`free`$ neutrons, comparisons involving different isospin channels for $`k`$ $`<<`$ 100 MeV/c should be performed with the greatest care as far as free neutron targets, or antiproton targets, are not available. One should therefore be aware that an isospin decomposition cannot be completely free from model dependence.

## 3 Nuclear shadowing effects.

Recent experimental data show that at antinucleon momenta (in the laboratory) below 60 MeV/c the nuclear shadowing effect is very strong. Below 60 MeV/c antiproton annihilation rates on Deuteron and 4He are smaller than on Hydrogen. The annihilation rate on 20Ne is larger but not that much. This and related phenomena has been discussed by us and other authors.

Rather independently of the mechanism underlying the annihilation process, it had been previously demonstrated that in the framework of the multiple scattering theory, of variational methods and of optical potential treatments one can predict such shadowing effects. It has been reported that also in the coupled-channel approach one can obtain the same result.

The fact that different methods lead to similar results suggested us to investigate the problem from a more general and qualitative, although less precise, point of view. In our work we have shown that due to the quantum uncertainty principle the $`\overline{n}`$-nucleus cross sections should be almost $`A`$-independent, apart for fluctuations due to nuclear surface effects. Consequently the $`\overline{p}`$-nucleus cross sections should depend on the target because of its electric charge only. The underlying argument is that most of the existing models (see the suggested references ) and analyses establish that the annihilation process takes place when the centers of mass of the antinucleon and of the target nucleus are at a relative distance $`d`$ such that $`R_{nucleus}`$ $`<`$ $`d`$ $`<`$ $`R_{nucleus}+\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ $``$ 1 fm (or smaller, depending on the model) does not depend too much on the target. So the annihilation is equivalent to a measurement of the projectile-target relative distance with uncertainty $`\mathrm{\Delta }`$ $`<`$ 1 fm, and this measurement is incompatible with a relative momentum $`<<`$ 200 MeV/c.

To see it another way, we distinguish between two classes of nuclear reactions. On one side, inelastic reactions where the entire nucleus is involved, as in compound nucleus reactions, but the underlying projectile-nucleon processes are elastic (e.g. neutron induced nuclear reactions). In this case the characteristic reaction region coincides approximately with the target nucleus. Then the uncertainty $`\mathrm{\Delta }`$ coincides approximately with the nuclear radius. On the other side, we find reactions where a strong inelasticity is present at the projectile-nucleon level. In this case reactions deep inside the nuclear volume are rare, the reaction region is a shell at the surface of the target nucleus, with thickness $`\mathrm{\Delta }`$, and $`\mathrm{\Delta }`$ is approximately the same for all the possible targets.

The consequence of the limitations imposed by the uncertainty principle is that for antinucleon momenta $`k`$ $`<<`$ $`1/\mathrm{\Delta }`$ the total reaction cross section becomes much smaller than its possible unitarity limit. This is also established by the well known low energy limit for the phase shifts: $`\delta _l`$ $``$ $`k^{2l+1}`$ for $`k`$ $``$ 0. The unitarity limit is reached when a partial wave is completely absorbed in the reaction process, which means $`exp(i\delta _l)`$ $`=`$ 0, i.e. $`Im(\delta _l)`$ $`=`$ $`\mathrm{}`$, so the unitarity limit cannot be attained at small enough $`k`$. Uncertainty considerations suggest that for $`k`$ $`>>`$ $`1/\mathrm{\Delta }`$ it is possible, for strong enough reactions, to saturate the unitarity limit, while for $`k`$ $`<<`$ $`1/\mathrm{\Delta }`$ we are in the situation where $`\delta _l`$ $`=`$ $`O(k^{2l+1})`$, whatever the strength of the reaction. A paradoxical consequence is that a smaller $`\mathrm{\Delta }`$ corresponds to what would be a stronger reaction at large energies, so that at low energies “stronger” interactions can lead to a smaller reaction rate. This fact can be verified in optical potential treatments.

On the ground that the projectile wavefunction $`\mathrm{\Psi }`$ is completely damped within a range $`\mathrm{\Delta }`$ (i.e $`|\mathrm{\Psi }|`$ is large for $`r`$ $`>`$ $`R_{nucleus}+\mathrm{\Delta }`$ and very small for $`r`$ $`<`$ $`R_{nucleus}`$) it is straightforward to demonstrate that for the scattering length $`\alpha `$ we have (approximately):

$`Im(\alpha )`$ $``$ $`\mathrm{\Delta }`$,

$`Re(\alpha )`$ $``$ $`+R_{nucleus}`$.

Indeed, the $`\mathrm{\Psi }`$ damping requirement implies for the logarithmic derivative $`|\mathrm{\Psi }^{}/\mathrm{\Psi }|`$ $``$ $`1/\mathrm{\Delta }`$. This is an obvious geometrical fact, but in more physical terms it is a consequence of the uncertainty principle. Together with the matching condition between the logarithmic derivatives of the free motion wavefunction and of the wavefunction in the annihilation region $`|\mathrm{\Psi }^{}/\mathrm{\Psi }|_{r=R_{nucleus}+\mathrm{\Delta }}`$ $`=`$ $`|kcotg\{k(R_{nucleus}+\mathrm{\Delta }\alpha )\}|`$ this leads to the previous $`\alpha `$-values in the limit $`k`$ $``$ 0.

These values of course are deduced from approximate equations, so they represent just estimates, however they suggest that the antineutron annihilation cross sections should not show a $`systematic`$ increase with the target mass number $`A`$. Such an increase could be present for antiproton annihilations, but because of Coulomb effects only. When going to any specific target nucleus, non-systematic effects could be present, especially related with the structure of the nuclear surface. An example is given in ref.. There, anomalous behaviors are related with the non-sphericity of the nucleus. Another exception should be represented by neutron-halo nuclei, because in this case the annihilation range could be much larger than 1 fm. Also the composition of the nuclear surface in terms of protons or neutrons could be important, since all models attribute a strong isospin dependence to the antinucleon-nucleon interaction.

The exposed mechanism has an interesting consequence in the case of optical potential analyses: an increase of the strength of the imaginary part of the optical potential can lead to a decrease of the consequent reaction rate at small momenta. In the above language, an increase in the potential strength leads to a decrease in the size parameter $`\mathrm{\Delta }`$, since the absorption of the projectile wavefunction takes place in a shorter range. Also modification of other parameters (radius, diffuseness, etc) leads to consequences that are not necessarily the most obvious ones. An example is given in the next section.

## 4 Comparison between $`\overline{np}`$ and $`\overline{p}p`$ total annihilation cross sections.

In fig. 1 we show the TNA and TPA on Hydrogen, together with two fits by energy-independent optical potentials. For the $`\overline{p}p`$ case the total interaction includes the electrostatic potential of a spherical charge distribution with radius 1.25 fm. This charge radius is $`\sqrt{2}`$ times the charge radius of the proton, to take into account both the proton and the antiproton extended charges. Details on the potential are given below. A third optical potential curve shows what would be the $`\overline{p}p`$ cross section in absence of electrostatic interactions, within the same optical model.

What we notice first is that TNA and TPA are reasonably similar for $`k`$ $`>`$ 200 MeV/c, but below 200 MeV/c the TNA falls clearly below the TPA, and this fall can not be justified in terms of Coulomb interactions only, although these contribute. We could speak of “low energy fall” with respect to a “background” which is similar to the $`\overline{p}p`$ case. Alternatively one can think that the TNA sums two contributions: a monotonous background that is proportional to the $`\overline{p}p`$ annihilation rate (Coulomb effects apart) but lower, and a broad peak (resonance?) at $`k`$ $``$ 200$``$300 MeV/c. This broad peak would be hidden by the isospin-0 background and by the low energy Coulomb rise in the TPA case. Some of the available data at momenta $`>`$ 400 MeV/c (see  for a more systematic review of these data) could support the “peak” idea, since at these momenta the $`\overline{n}p`$ annihilation rate seems to be lower than the $`\overline{p}p`$ one, but it is not unequivocally clear how much lower. Indeed, one $`\overline{n}p`$ point at 700 MeV/c by  is clearly below the $`\overline{p}p`$ annihilation, while data reported in  are not very different from the corresponding $`\overline{p}p`$ data in the region 300-500 MeV/c. On the other side, some data at very low energy support the “gap” idea, since they show a rise of the $`\overline{n}p`$ annihilation rate $`\beta \sigma `$ below 50 MeV/c. These data have been reported in the figure, and the one at 20 MeV/c lies on the optical potential extrapolation of the $`\overline{p}p`$ data, once the electrostatic part of the potential has been removed. The authors of ref. also report some potential predictions that support the existence of a minimum for the $`\overline{n}p`$ annihilation rate $`\beta \sigma _{ann}`$ at 50-100 MeV/c, and a set of scattering length predictions that show no great differences between $`\overline{p}p`$ and $`\overline{n}p`$ expected scattering lengths. They also report a collection of pre-Obelix $`\overline{n}`$ data.

We must remark that the normalization of antineutron fluxes is a delicate matter, with the consequence that both a comparison of antineutron data coming from different experiments, and a comparison of antineutron and antiproton data, must be taken with the greatest caution. So we will generically speak of a “gap/peak structure” of the $`\overline{n}p`$ annihilation rate compared to the $`\overline{p}p`$ one. We will show in the next section that a comparison with nuclear data supports the existence of this gap/peak structure. Whatever the interpretation, there are some important remarks:

1) The strong oscillation of the TNA with respect to the TPA (see also fig.4) is an evidence of the fact that different physical mechanisms are dominating in the two cases in the considered momentum range. If a resonance is present this is obvious. If there is no resonance, according to the partial wave analyses presented in refs. the TNA P-wave has its maximum at $`k`$ $``$ 250 MeV/c where it attains the value 120-140 mb. Our optical potential fit (which seems good for the TPA) fixes the TPA P-wave maximum to about 80 mb at 150 MeV/c. Semiclassical intuition fixes the corresponding impact parameters $`b`$ $`=`$ $`L/k`$ to 1.25 fm (TPA) and 0.8 fm (TNA). Then, purely geometrical considerations suggest a ratio, between the size of the two P-wave outcomes at their peak momentum, TNA/TPA $``$ $`(0.8/1.25)^2`$ $``$ 0.4 instead of the found value 1.6. The most obvious conclusion would be that the P-wave interaction is much more effective in the $`\overline{n}p`$ case, where it acts at smaller impact parameters (by a factor 0.6). This short range interaction does not affect the S-wave contribution, which is clearly larger in the $`\overline{p}p`$ case (even after subtracting Coulomb effects) at momenta over 30 MeV/c. At lower momenta (that means longer range) the two S-wave contributions could be equal, according to and references therein. So the geometry and strength of the isospin-1 and isospin-0 channel interactions must be different, as far as these semiclassical considerations can be trusted.

2) Absurd as it may sound, the fact that at small momenta TNA are sensibly smaller than TPA could mean that the interactions are stronger in the $`\overline{n}p`$ case than in the $`\overline{p}p`$ one. This can be immediately seen in the choice of the optical model parameters. We used a Woods-Saxon form, with all parameters, but one, equal for TPA and TNA. TPA: imaginary strength 8000 MeV, real strength 46 MeV (attractive), imaginary radius 0.52 fm, real radius 1.89 fm, real and imaginary diffuseness 0.15 fm. These values fall in the ranges used by previous authors to fit elastic data, and its imaginary part is pretty similar to other previously used optical potentials (e.g. the one of ref.). Many other sets of parameters, and also different potential shapes, can lead to similar results for the available low-energy data. Starting from this peculiar set, to obtain the TNA curve it is sufficient to increase the imaginary radius to 0.75 fm. Alternatively, one can obtain the TNA curve by leaving the imaginary radius at 0.51 fm and increasing the imaginary strength to about 12000-16000 MeV. So, to get a smaller cross section we need a “stronger” potential. This “inversion” behavior is another manifestation of the general shadowing mechanism discussed in section 3. In an optical potential model at small $`k`$ an increase, e.g., of the imaginary strength $`W`$ produces an increase of the reaction cross section for small values of $`W`$ only. After a certain threshold, further increases of $`W`$ lead to a decrease of the reaction cross section. At a qualitative level this can be explained by the uncertainty principle: an increase in the strength of the imaginary part of the potential decreases the thickness of that spherical shell, surrounding the proton/nucleus target, where annihilations are supposed to take place. This introduces larger gradients in the projectile wavefunction. These larger gradients imply a larger logarithmic derivative, and therefore a smaller (imaginary part of the) scattering length.

3) In the case of $`\overline{p}p`$ annihilations below 600 MeV/c many quite different choices of potential parameters can fit them (including a pure imaginary potential) without introducing a dependence of the parameters on energy. All these fits share some common features. They reproduce the data very well below 150 MeV/c and above 300 MeV/c. The fits remain satisfactory up to 600 MeV/c. In the region 150-300 MeV/c the fits stay slightly below the data. Alternatively, it is possible to choose the parameters so as to reproduce the data from 150 to 600 MeV/c, while at lower momenta the fits are above the data. So, we could say that an energy independent optical potential can produce a “background” around which a slight oscillation is present (positive for $`k`$ in the range 150-300 MeV/c, or negative for $`k`$ $`<`$ 150 MeV/c). In the TNA case the situation is qualitatively the same, but the oscillation is much stronger and more evident. So it seems that an energy-independent optical potential fit can be well related with the isospin-0 channel, which seemingly consists of the “background” only. A background with the same shape is also present in the isospin-1 channel, but this channel also contains an oscillation that is not reproduced by our set of optical potentials. As already remarked, we have no means to establish the nature of this oscillation.

4) We have also tried to change our optical potential shape so as to produce resonances that fit the gap/peak structure which is seen in the $`\overline{n}p`$ annihilation. We can not exclude that more attempts can lead to a good reproduction of the data in fig.1 (in the described “inversion plus resonance” regime there is no predictable relation between the changes of the potential and the shape of the output, which makes such attempts rather frustrating). We have not been able to reproduce satisfactorily these data, and perhaps more sophisticated potentials are necessary. It was however easy to get shapes that $`qualitatively`$ looked similar. We tried with two kinds of potential: (i) Woods-Saxon (with attractive elastic part) (ii) Woods-Saxon like the previous one, plus a repulsive elastic surface barreer of gaussian shape (of the kind $`exp[(rr_o)^2]`$). In both cases to get a broad peak where it must be it was necessary to enlarge the range of the elastic part to 2.5-3 fm. This is probably due to the fact that, despite the smallness of the radius of the imaginary part, most of the annihilations take place at $`r`$ $``$ 1 fm. So only a much wider attractive well can produce metastable states, within such models. Perhaps in a more sophisticated potential model, where different spin channels are clearly separated, it is possible to get a resonance in a spin channel where annihilation is less effective, so without the need for such a long range for the potential. We think that these interesting data really demand for a higher level potential analysis.

## 5 Annihilation on nuclei. Shadowing in the charge ratio.

Because of the shadowing, a complete analysis of annihilations in nuclei below 60 MeV/c is pretty complicated. A simple comparison of different sets of data reveals interesting features, however, and we can also hope to understand something by a PWIA analysis, especially at momenta over 100 MeV/c. In particular, the energy dependences of the $`\overline{n}p`$ annihilation cross section, of the $`\overline{n}`$-nucleus one, and of the ratio of $`\overline{p}n`$ to the $`\overline{p}p`$ annihilation in nuclei, show evident correlations. Of course nuclear shadowing effects are completely absent in PWIA, and this must be taken into account, taking with great care the calculations below 100 MeV/c. At larger momenta the main missed effect is the eclipse effect which implies an almost energy independent surface absorption factor $``$ 0.5. We take this factor into account by renormalizing our PWIA predictions with one datum at a specified momentum value over 200 MeV/c.

At small momenta however (below $``$ 200 MeV/c) shadowing is something more complicated than a mere eclipse effect. The eclipse effect could never bring $`\overline{p}`$ annihilation rates to decrease with increasing $`A`$. So we may only $`define`$ as “shadowing” the departures from PWIA predictions, and remark that they are consistent. As we show later, anyway, a certain part of the smallness of the low energy $`\overline{p}`$-Deuteron annihilation rate with respect to the $`\overline{p}p`$ one is not a shadowing effect, but just due to the smallness of the $`\overline{p}n`$ annihilation rate below 150 MeV/c.

In fig.2 we show the $`\overline{n}`$-nucleus data, relative to momenta between 180 and 280 MeV/c, for several nuclear species from 12C to 207Pb (see  for more details). These data are consistent with an $`A^{2/3}`$ law, expressing dominant surface absorption, $`among`$ $`them`$. However the $`A^{2/3}`$ law found in these data cannot be generalized to $`\overline{n}p`$ or $`\overline{p}p`$ (Coulomb subtracted) data in the laboratory frame, while this generalization is only approximately possible in the center of mass frame. At 180 MeV/c, the 12C or the 207 Pb TNA, divided by $`A^{2/3}`$, give about 150 mb. The corresponding $`\overline{n}p`$ and (Coulomb subtracted) $`\overline{p}p`$ cross sections are $``$ 200 mb at the same laboratory momentum. At the same center of mass momentum, which corresponds to 360 MeV/c in the laboratory, the $`\overline{n}p`$ total annihilation cross section is about 130 mb, and the $`\overline{p}p`$ one is slightly smaller than the $`\overline{n}p`$ one. So, in both cases we find deviations from the $`A^{2/3}`$ law at $`k_{cm}`$ $``$ 100 MeV/c, but these deviations are of opposite sign if the data are compared at the same laboratory or center of mass momentum. It is difficult to establish whether it is more proper to compare data relative to different nuclear species at the same laboratory or center of mass momentum. Impulse approximation based models suggest the same laboratory momentum, compound nucleus inspired models prefer the opposite choice. Indeed, in the former case direct momentum exchange is between the projectile and one of the target nucleons only, whose average initial momentum is zero in the laboratory. In the latter, the projectile momentum is given directly to the entire target nucleus.

An analogous, systematic comparison of all the available $`\overline{p}`$-nucleus data has been presented elsewhere.

Very recently a new series of $`\overline{n}`$-nucleus data down to 70 MeV/c has been presented, which show a respected $`A^{2/3}`$ law between different nuclear species. Since these data have been presented as “preliminary”, we can’t discuss them in detail.

Considering the energy behavior of the TNA on nuclei shown in fig.2, the most interesting feature is the systematic presence of a shoulder. However this phenomenon seems less evident in the just quoted preliminary data. TPA data in this region are not enough to to give us shape information. However, the optical potential calculations that we have reported in  suggest that the Coulomb enhancement at low energies can easily hide such details in the energy dependence of the annihilation data.

Different relevant mechanisms must be taken into account in the momentum range 100-300 MeV/c:

1) The tail of the low energy shadowing phenomenon. Presently we have no knowledge of TPA and TNA on heavy nuclei at $`k`$ $`<<`$ 100 MeV/c, so we do not know whether the discussed shadowing is a systematic phenomenon (as we suggest in the previous section 2) or is simply peculiar of some light nucleus. And in light nuclei we do not have enough TPA data in the region 100-300 MeV/c to understand completely the transition between large and low energy regimes. There are several sets of data, but regretfully we do not have a satisfactorily continuous set of data relative to the same projectile on the same target (with the exception of the $`\overline{p}p`$ and $`\overline{n}p`$ cases).

3) The nuclear Fermi motion. The Fermi momentum, which represents a rough cutoff of the nuclear momentum distribution, is $``$ 200 MeV/c. This means that for $`k`$ $`<`$ 200 MeV/c the statistical distribution of the relative $`\overline{n}`$-nucleon motion contains the zero-energy point, and that it is easy for a target nucleon to be faster than the projectile. So the projectile momentum 200 MeV/c represents a kind of special borderline in the physics of annihilations on nuclei.

4) The $`\overline{n}p`$ annihilation cross section presents the above discussed gap/peak structure below 300 MeV/c. The nuclear Fermi motion will soften this structure, but perhaps not enough to make it completely disappear.

A PWIA calculation can give some qualitative understanding about the roles of the gap/peak structure and of the Fermi motion in the nuclear TNA, and also about the relative role of neutrons and protons in the target.

In fig.3 we present a PWIA calculation of the TNA on 40Ca. In fig.4 we show a PWIA calculation of the $`F_A`$ ratio on Deuteron and 4He. We recall that we have defined $`F_A`$ as the ratio between $`\overline{p}`$ annihilations on neutrons and protons inside a nucleus with mass number $`A`$. In figs. 5 and 6 we show PWIA calculations of $`\overline{p}`$ annihilations on deuteron and 4He.

As a starting point we have used an equation contained in . Although this equation has been justified there by a long series of mathematical passages, its meaning is quite simple and can stand without demonstration: at PWIA level, any event rate on a nuclear target, leading to a well defined final state, is an incoherent sum of the rates of all the possible $`\overline{n}n`$ and $`\overline{n}p`$ events leading to the same final state, weighted by the nuclear momentum distribution of the target nucleon. By “rate” we mean a cross section divided by the corresponding flux of colliding particles. Under the simplifying assumption that different final states sum incoherently we simply have:

$$\beta (k,0)\sigma (k)=N(\stackrel{}{k}^{})\beta (\stackrel{}{k},\stackrel{}{k}^{})\sigma (|\stackrel{}{k}\stackrel{}{k}^{}|)d^3k^{},$$
(1)
where $`\beta (\stackrel{}{k},\stackrel{}{k}^{})`$ is the relative velocity between two colliding particles with momenta $`\stackrel{}{k}`$ and $`\stackrel{}{k}^{}`$.

This nuclear PWIA average of the antinucleon-nucleon processes clearly misses any eclipse or shadowing effect (which, e.g., produces a TNA proportional to $`A^{2/3}`$ at large energies). So, the rise of the PWIA curve in fig.3 below 150 MeV/c is unreliable (it is just proportional to the average rise of antiproton and antineutron annihilation cross sections on a proton). At momenta larger than 150 MeV/c the PWIA should be more reliable, apart for a slowly energy dependent eclipse factor, and could clarify the effects of the Fermi motion and the different role of protons and neutrons in the target.

For the 40Ca shells we have used harmonic oscillator eigenfunctions. Some attempts with different kinds of states have shown no qualitative differences, unless the nuclear size is changed to unphysical values. In addition, the most external wavefunctions do not show great differences (in the energy dependence of the outcome) with respect to the S-wave ground state. This is important, because in a more proper DWIA (distorted wave impulse approximation) treatment the internal shells would be scarcely involved in the annihilation process. For Deuteron and 4He (in fig.4) experimental single particle momentum distributions have been used. Again, within this simple model we did not find a great sensitivity of the results on the peculiarities of the momentum distribution, taking into account that many details that could be present in the numerator and in the denominator are lost in the ratios of fig.4. With Deuteron, we surely missed the angular dependence of the large momentum (D-wave dominated) part of the distribution.

To calculate the PWIA curves of figure 3 we need $`\overline{n}p`$ and $`\overline{n}n`$ annihilation rates. For fig.4 we need $`\overline{p}p`$ and $`\overline{p}n`$ annihilation rates, both subject to the effect of the Deuteron or 4He charge. The $`\overline{n}n`$ annihilation rate has been assumed to be equal to the Coulomb subtracted $`\overline{p}p`$ one (i.e. the continuous curve of fig.1). In fig.4 we use the same cross sections, so we neglect the effects of the nuclear electrostatic attraction. On the basis of the discussion in section 2 we assume that they do not affect the $`F_A`$ ratio, although they surely enhance both the $`\overline{p}p_{innucleus}`$ and the $`\overline{p}n_{innucleus}`$ rates. Actually a look at the data in fig.4 shows that (with large error bars) $`F_4/F_2`$ $``$ 0.48/0.78 $``$ 0.61 at zero energy, and 0.54/0.82 $``$ 0.66 at about 600 MeV/c. In the case of a different action of the nuclear electrostatic potential on the $`\overline{p}p`$ and $`\overline{p}n`$ cross sections, this effect would be practically absent at 600 MeV/c and would depress the $`F_4`$ ratio with respect to the $`F_2`$ one at momenta $`<<`$ 100 MeV/c. The large error bars and data fluctuations do not allow for easy conclusions, however there are no special reasons to assume a different Coulomb effect on $`\overline{p}p`$ and $`\overline{p}n`$ reactions, if both the proton and the neutron are bound in the same nucleus.

The Obelix $`\overline{n}p`$ data cover the range 60-400 MeV/c, and the $`\overline{p}p`$ ones start at 30 MeV/c. We therefore had to assume values for these cross sections for momenta outside these ranges. With $`\overline{p}p`$ data at momenta over 200 MeV/c we have relied on the data in ref., and at momenta below 30 MeV/c on our optical potential fit (after subtraction of Coulomb forces as in fig.1). For the $`\overline{n}p`$ case, in the region $`k`$ $`>`$ 400 MeV/c the $`\overline{n}p`$ data have been assumed to be equal to the $`\overline{p}p`$ ones, as suggested by the measurements in ref. and  for momenta over 300 MeV/c (see e.g. , fig.9, for a comparison of the $`\overline{n}p`$ and $`\overline{p}p`$ data in that momentum range). A look at fig.4 does not help too much. Indeed, we may say that the nuclear Fermi motion approximately averages the “free” annihilation rates within a 200 MeV/c range. Then the point of view of equal “free” $`\overline{p}p`$ and $`\overline{p}n`$ annihilation rates at momenta in the range 300-600 MeV/c can be supported by the absence of a clear decreasing trend in the set of Deuteron points over 300 MeV/c. On the contrary, the Helium point at 600 MeV/c supports the possibility of a decrease of the free $`\overline{p}n/\overline{p}p`$ ratio, as suggested by the measurement in ref.. Anyway, the scatter of data in the deuteron case, and the lack of data in the 4He case prevent us from precise conclusions. In the limit of zero energy, we have adopted the low energy parametrization of the $`\overline{n}p`$ cross section contained in references. It implies a decrease of the $`F_A`$ ratios at very low energies, which is more suitable for fitting the zero-energy 4He point. Adopting the low energy parameters suggested by refs. would on the contrary rise the zero energy $`F_A`$ value, which is more coherent with the two zero-energy Deuteron points.

The upper curve in fig.3 can be taken as a qualitative confirmation of the effect of the $`\overline{n}p`$ gap/peak structure in producing a very soft shoulder in the $`\overline{n}`$nucleus annihilations at 150-300 MeV/c In fig.3 we also show the separate contributions of the $`\overline{n}p`$ and $`\overline{n}n`$ reactions. Evidently the $`\overline{n}n`$ reaction alone would not interrupt the trend of a positive second derivative that characterizes expectance and data at both lower and larger momenta. Within PWIA the shoulder is a direct consequence of the gap/peak structure of the $`\overline{n}p`$ cross section. The Fermi motion, as predictable, softens this gap/peak structure, but not completely. We must also observe that the produced shoulder is much less pronounced than most of the ones seen in fig.2. Our attempts with different shapes for the nuclear shells show that the shoulder evidence is enhanced by a larger radius of each nuclear shell (i.e. a narrower momentum distribution), however one must arrive at slightly unrealistic nuclear dimensions in order to obtain more evident shoulders.

The strength of the rise of the PWIA cross sections below 100 MeV/c is probably unphysical, since we know that shadowing could be very active in that region (although the quoted recent preliminary data  confine shadowing to $`k`$ $`<`$ 70 MeV/c and with good approximation reproduce the predictions of fig.3). In ref. some optical potential shapes possibly corresponding to low energy $`\overline{n}`$-nucleus annihilations are reported, and they are almost energy-independent in the region 100-300 MeV/c. These optical potentials were based on Coulomb-subtracted $`\overline{p}p`$ data and on nuclear cross sections at momenta over 200 MeV/c and so have no relation with that gap/peak structure that is peculiar of the $`\overline{n}p`$ interaction. The data reported in  confirm the IA behavior reported in fig.3 down to 70 MeV/c, rather than the these optical potential predictions.

The systematic presence of the gap/peak structure in the TNA/TPA ratio is much more evident in fig.4, where we report the $`F_A`$ ratio for $`A`$ $`=`$ 1 ($`\overline{n}p`$ data coming from the Obelix experiment divided by the Coulomb subtracted $`\overline{p}p`$ annihilation cross section represented by the continuous line in fig.1), for $`A`$ $`=`$ 2 (Deuteron) and $`A`$ $`=`$ 4 (4He). These Deuteron and 4He data come from several experiments and are reviewed in ref..

Looking at fig.4 one immediately notices that all deuteron points lie above any of the 4He points, and the proton target data are, in the average, clearly over the nuclear target data (we remark that, although the proton target ratios have been calculated using the “Coulomb subtracted” $`\overline{p}p`$ annihilation cross section, the effect of the Coulomb subtraction is negligible above 100 MeV/c). The second observation concerning the data reported in fig.4 is that a certain degree of correlation between the shapes of the three groups of data is present.

Since $`Z/A`$ $`=`$ 1 both in Deuteron and 4He, and the average value of the free $`F_1`$ ratio is about 1, clearly PWIA cannot reproduce the average value of these nuclear $`F_A`$ distributions. For producing the fits contained in fig.3 we have multiplied the $`\overline{n}p`$ total annihilation cross section by a constant factor $`\alpha `$, which is 0.8 for Deuteron target, and 0.6 for 4He target. At the largest energies of the considered range the PWIA ratio should be reliable. The known eclipse effect should more or less cancel between numerator and denominator. Anomalous shadowing of the kind discussed in section 3 is not strong over 100 MeV/c. So the reason for $`\alpha _{{}_{}{}^{4}He}`$ $`<`$ $`\alpha _D`$ $`<`$ $`\alpha _p`$ $``$ 1 is unexplained. A possibility is final state rescattering. This should not be very effective in total annihilation cross sections, but in the $`F_A`$ case it can modify the number of charged particles that is used to decide whether the annihilation was on a proton or on a neutron. It is not unlikely that this $`A`$-dependent renormalization effect has this origin, however one should analyze carefully how each experimental point has been measured.

In the case of 4He the PWIA fit reproduces satisfactorily the behavior of the data. We may say that the gap/peak structure of the free $`\overline{n}p/\overline{p}p`$ ratio is reproduced by the $`F_4`$ ratio. In the Deuteron case at zero-energy the fit is not satisfactory, while over 200 MeV/c some correlation seems to be present. Somehow the presence of the PWIA peak at 280 MeV/c (which is a consequence of the gap/peak structure of the free ratio) is reflected in the presence of a maximum in the data at 300 MeV/c, but the large error bars and the spread of the data do not allow for more quantitative conclusions.

In the last two figures we show some PWIA fits of the $`\overline{p}`$ annihilation cross sections on Deuteron and 4He. To obtain these curves the same above method has been used, but the final result has been multiplied, to take the nuclear charge into account in agreement with ref., by the factor $`1+ZC\beta ^{1.4}`$, where $`Z`$ is the nuclear charge, $`\beta `$ is the $`\overline{p}`$-nucleus relative velocity and the constant $`C`$ is characteristic of the target nucleus. The exponent ($``$1.4) has no special physical meaning, since the above $`1+ZC\beta ^{1.4}`$ factor is only a fitting relation that reproduces within some percent the Coulomb enhancement factor in the range 40-400 MeV/c. The constant $`C`$ is 0.0060 for Deuteron, and 0.0040 for $`{}_{}{}^{4}He`$. To take surface absorption into account, the two PWIA fits have been renormalized by a constant factor, the best to reproduce the data (or some data) over 200 MeV/c. We remark that: (i) this Coulomb factor is not relevant from 100 MeV/c onward, (ii) it is not very different between light nuclei at the same $`\overline{p}`$ momentum in the laboratory frame, (iii) in Hydrogen, Deuteron and 4He it is near 1.5 at 50 MeV/c and 2 at 25 MeV/c (laboratory frame).

In the two figures we also report what would be the PWIA distribution if the target nucleus were composed by neutrons only or by protons only (without touching the Deuteron and 4He charges, which are taken into account by the $`1+ZC\beta ^{1.4}`$ factor in any case). The first evident observation is that the effect of the $`\overline{n}p`$ gap/peak structure is much more evident in the Deuteron case. This is due to the compactness of the 4He structure, and to the exactly opposite peculiarity of the Deuteron structure. A more compact space structure implies a broader momentum distribution. The Fermi motion tends to average out the details of the momentum dependence of the nucleon-antinucleon annihilation rates. This mechanism is more effective with broader nuclear momentum distributions. Another observation is that, if we define shadowing as the departure between the observed data and Impulse Approximation predictions, this departure becomes evident in both cases below 100 MeV/c, but more evident in the 4He case. In the Deuteron case the shadowing phenomenon is less dramatic, but would be overestimated by not taking into account that one of the two nucleons composing this nucleus is a neutron.

Concluding this section, we may say that some evident correlation exists between nuclear data and the behavior of the $`\overline{n}p`$ annihilation rate. In particular, it would be difficult to justify the shape of the $`F_4`$ and $`F_2`$ ratios in absence of the gap/peak structure of the free $`\overline{n}p/\overline{p}p`$ ratio. There is some chance that the shoulder found in the $`\overline{n}`$-nucleus annihilation cross sections is related with the same gap/peak structure, however in this case the situation is much less clear. Still to be clarified is the fact that the energy-averaged value of the $`F_A`$ ratio decreases at increasing $`A`$. The last two figures also show that the energy dependence of the nuclear annihilation rates can be partly explained in terms of the energy dependence of the $`\overline{p}`$ annihilation cross section on a free neutron. A consistent part of this energy dependence, however, cannot be explained in these terms and we may say that a strong energy-dependent shadowing is present.

## 6 Conclusions.

We have presented an analysis comparing low energy $`\overline{n}p`$ annihilation data with other available data on antinucleon-nucleon and antinucleon-nucleus annihilation. The different sets are all consistent with a gap/peak structure in the isospin-1 channel, over a regular “background” which can be reproduced by energy-independent optical potentials of simple form. This structure is not too pronounced but evident, and we have also shown that it affects the nuclear annihilation rates, whose behavior can be $`partly`$ explained in terms of the behavior of the antiproton-neutron annihilation rate.

We have also discussed the shadowing phenomena related with antinucleon-nucleus annihilation, and reported an interpretation in terms of the quantum uncertainty principle. The comparison of impulse approximation predictions with the nuclear data suggests that a strong energy-dependent shadowing is present at low energies: also after including absorption corrections, the data can not be fully explained in terms of individual antinucleon-nucleon annihilations.


format: text

text:
# HIERARCHICAL GALAXY FORMATION AND SUBSTRUCTURE IN THE GALAXY’S STELLAR HALO

## 1. Introduction

The origin of the Milky Way’s stellar halo is a long-standing astronomical problem. The poles of the debate are defined by the monolithic collapse model of Eggen, Lynden-Bell, & Sandage (1962) and the chaotic accretion model of Searle (1977) and Searle & Zinn (1978). The Searle & Zinn picture has gained currency in recent years in part because of growing recognition that the halo and bulge are distinct components that may have different formation mechanisms (see, e.g., the reviews by Wyse 1999ab) and in part because of “smoking gun” evidence that includes the tidally distorted Sagittarius dwarf galaxy (Ibata et al. 2000a and references therein) and the presence of extra-tidal stars around many dwarf spheroidal satellites (Gould et al. 1992; Irwin & Hatzidimitriou 1995; Kuhn, Smith & Hawley 1996). The Searle & Zinn scenario also bears a strong anecdotal resemblance to the hierarchical galaxy formation scenario characteristic of inflationary cold dark matter (CDM) cosmological models. In this paper, we make the connection between CDM cosmology and hierarchical stellar halo formation much more explicit, by presenting a simple but quantitative model for halo formation in the CDM framework and obtaining predictions for the degree of residual substructure in the outer halo.

Previous studies of stellar halo formation in the hierarchical framework have focused on the fossil evidence for satellite disruption preserved in phase space substructure of the halo stars (e.g., Johnston, Spergel & Hernquist 1995; Helmi & White 1999). These studies were aimed primarily at exploiting surveys of the halo in the solar neighborhood (e.g., Arnold & Gilmore 1992; Preston, Beers & Shectman 1994; Majewski, Munn & Hawley 1994, 1996). Wide-angle, deep, multi-color surveys, such as the Sloan Digital Sky Survey (SDSS; York et al. 2000), open up new avenues for studying the structure of the stellar halo. RR Lyrae stars, detected by their variability and color, can provide a three-dimensional map of the distribution of halo stars. The more general population of A-colored halo stars can be used for the same purpose; relative to RR Lyrae they have the advantages of greater numbers and detectability in a single observation epoch but the disadvantage of being less precise standard candles. Studies of RR Lyrae stars and A-colored stars in SDSS commissioning data have already revealed two large substructures in the outer halo (Ivezić et al. 2000; Yanny et al. 2000). The photometric depth of the SDSS and the intrinsic brightness of RR Lyrae and A stars allows a probe of halo structure out to large distances, $`75`$ kpc, and the restricted absolute-magnitude range of RR Lyrae and A stars prevents 3-dimensional substructure from being washed out by projection. Carbon star surveys (e.g., Ibata et al. 2000a) and surveys of giant stars (Majewski et al. 2000) offer similar prospects.

The model presented in this paper is an extension of our previous work aimed at explaining the observed abundance of dwarf satellite galaxies in the Local Group within the CDM framework (Bullock, Kravtsov & Weinberg 2000; hereafter BKW). In BKW we showed that the observed shape and amplitude of the velocity function of dwarf satellite galaxies in the Local Group can be explained if gas accretion and star formation are suppressed in low-mass dark matter clumps after intergalactic gas is reheated during the epoch of reionization. In this picture, the observed dwarf satellites around our Galaxy are those that assembled a large fraction of their mass before reionization and survived the decay of their orbit as a result of dynamical friction and avoided tidal disruption by the Milky Way potential. As shown in Fig. 1 of BKW, the number of tidally disrupted objects is similar to the number of surviving dwarf satellites. Here we show that the disrupted satellites produce a population of stars whose total mass and radial profile are consistent with observations of the Milky Way’s stellar halo. In the outer halo, where the number of contributing satellites is relatively small and the orbital times are long, the model predicts substantial substructure, which should be detectable with the SDSS and other deep, wide-angle surveys.

Our model is simple and is unlikely to be accurate in full quantitative detail. However, the qualitative predictions should be characteristic of conventional CDM cosmologies combined with straightforward assumptions about the star formation in low-mass dark matter potential wells. If the observed stellar halo is found to be radically different from these predictions, it will mean that either our star formation assumptions or the CDM predictions for hierarchical small-scale structure in the dark matter distribution are incorrect. In this sense, studies of the stellar distribution in the outer halo can play a valuable role in testing more general ideas about galaxy formation.

For convenience, we will focus on predictions for the RR Lyrae distribution. RR Lyrae are especially useful probes of the stellar halo because they are relatively easy to identify, they are luminous enough to be detected out to large distances ($`r100`$ kpc), and they are nearly standard candles and therefore yield 3-dimensional maps. RR Lyrae are believed to be good tracers of the more general halo stellar distribution, and they have often been used in kinematic studies of the halo (Hawkins 1984). Finally, while RR Lyrae are numerous enough to trace halo substructure, they are rare enough that we can construct numerical realizations that contain every individual RR Lyrae star. The stellar density fluctuations predicted by our model are far in excess of Poisson fluctuations, so it is straightforward to scale our predictions to other halo tracers like blue horizontal branch stars or carbon stars, just by putting in appropriate stellar population weights.

The remainder of the paper is organized as follows. In § 2 we will describe our model and assumptions. Specifically, we will describe our semi-analytic method for following the accretion and orbital evolution of satellite galaxies in § 2.1 and our modeling of stellar debris of disrupted satellites in § 2.2. We present our results in § 3. We finish with discussion and conclusions in § 4.

## 2. Method

### 2.1. Accretion and orbital evolution of satellites

We use a semi-analytic method to trace the accretion history and orbital evolution of satellite galaxies within a typical Milky Way-size dark halo.111In this paper, the term “halo” sometimes refers to a dark matter halo and sometimes to a stellar halo. Usually the meaning is clear from context, but we will specify “dark” or “stellar” where necessary. The term “subhalo” always refers to a dark matter halo, one that is accreted into a larger dark matter halo before redshift zero. Detailed description of the model is given in BKW. Here, we briefly review its essential aspects. The model uses the extended Press-Schechter formalism (Bond et al. 1991; Lacey & Cole 1993) to construct the accretion history for each galactic halo. The mass of a halo in this formalism is accumulated via accretion of individual subhalos of different masses, and we keep track of all the accreted subhalos down to some minimum mass. The second part of our model is a semi-analytic prescription for orbital evolution of the accreted subhalos. This prescription is used to determine whether a subhalo survives to the present day, is tidally disrupted, or is dragged into the central galaxy by dynamical friction. Only subhalos that accreted a significant fraction ($`\text{ }>\text{ }0.20.3`$) of their mass before intergalactic gas was reheated during the epoch of reionization are assumed to host luminous galaxies. In the model presented in this paper, the stellar halo of Milky Way-like galaxies is formed from the debris of those subhalos that once hosted luminous galaxies but were tidally disrupted before the present day. For our analysis, we adopt a flat CDM model with a non-zero vacuum energy and the following parameters: $`\mathrm{\Omega }_m=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,h=0.7,\sigma _8=1.0`$, where $`\sigma _8`$ is the rms fluctuation on the scale of $`8h^1`$ Mpc, $`h`$ is the Hubble constant in units of $`100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, and $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ are the density contributions of matter and the vacuum respectively in units of the critical density.

We assume that the density profile of each dark matter halo is described by the NFW profile (Navarro, Frenk, & White 1997): $`\rho _{\mathrm{NFW}}(x)x^1(1+x)^2`$, where $`x=r/r_s`$, and $`r_s`$ is a characteristic inner radius. Given a halo of mass $`M_{\mathrm{vir}}`$ at redshift $`z`$, the model of Bullock et al. (2000a) supplies the typical $`r_s`$ value and specifies the profile completely. The circular velocity curve, $`v^2(r)GM(r)/r`$, peaks at a value $`v_m`$ at a radius $`r_\mathrm{m}2.16r_\mathrm{s}`$.

We use the merger tree method of Somerville & Kolatt (1999) to construct mass growth and halo accretion histories for an ensemble of galaxy-sized dark matter halos. We start with halos of mass $`M_{\mathrm{vir}}=1.1\times 10^{12}h^1M_{}`$, at $`z=0`$, and trace subhalo accretion histories back to $`z=10`$. We record the mass growth for the primary halo, $`M_{\mathrm{vir}}(z)`$, as well as the mass of each accreted subhalo, $`M_a`$, and the time of its accretion, $`t_a`$ (or $`z_a`$). We assign the subhalo $`v_m`$ according to the mass-velocity relation at the epoch of accretion. For the results presented below, we use 100 ensembles of formation histories for galactic host halos. We obtain very similar results if this number is increased.

Each subhalo is assigned an initial orbital circularity, $`ϵ`$, defined as the ratio of the angular momentum of the subhalo to that of a circular orbit with the same energy, $`ϵJ/J_c`$. We choose $`ϵ`$ randomly in the range $`0.11.0`$ (Ghigna et al. 1998). To determine whether the accreted halo’s orbit will decay, we use Chandrasekhar’s formula to calculate the decay time, $`\tau _{DF}^c`$, of the orbit’s circular radius $`r_c`$ — the radius of a circular orbit with the same energy as the actual orbit. Each subhalo is assumed to start at a randomly assigned radius $`r_c^a=(0.40.75)R_{\mathrm{vir}}(t_a)`$, where $`R_{\mathrm{vir}}(t_a)`$ is the virial radius of the host halo at the time of accretion. We determined this distribution of circular radii by measuring the range of binding energies of subhalos in the ART simulations used by Klypin et al (1999a). Once $`\tau _{DF}^\mathrm{c}`$ is known, the decay time for the given circularity is $`\tau _{DF}=\tau _{DF}^\mathrm{c}ϵ^{0.4}`$ (Colpi et al. 1999). If $`\tau _{DF}`$ is smaller than the time left between accretion and $`z=0`$, $`\tau _{DF}t_0t_a`$, then the subhalo will merge with the central object. In our modeling of the stellar halo, we do not consider the contribution due to galaxies that subsequently merge with the central object. Due to the rapid decay of the orbits, any debris associated with these objects will likely remained confined within the radius where stripping first becomes important, typically $`r10`$ kpc. For this reason, we consider the predictions for the stellar distribution only for galactocentric radii $`r>10`$ kpc.

If $`\tau _{DF}`$ is too long for the orbit to have decayed completely ($`\tau _{DF}>t_0t_a`$), we check whether the subhalo would have been tidally disrupted. We assume that the halo is disrupted if the tidal radius becomes smaller than $`r_\mathrm{m}`$. The tidal radius, $`r_t`$, is determined at the pericenter of the orbit at $`z=0`$, where the tides are the strongest, following Klypin et al. (1999b). If $`r_tr_{\mathrm{max}}`$ we declare the subhalo to be tidally destroyed and record its mass and orbital parameters so that we may model the evolution of its tidal debris (§2.2).

The resulting average mass functions for all accreted halos and the subset of the halos that were tidally disrupted are shown in Figure 1 with the thin dashed and solid lines, respectively. The error bars represent the run-to-run dispersion over 100 realizations. For comparison, the dotted line corresponds to halos that were dragged to the halo center as a result of dynamical friction. As expected, the dynamical friction is more efficient for massive subhalos. The mass function of surviving subhalos is similar to that of the disrupted halos, but we omit it here to preserve visual clarity.

As in BKW, we assume that of all subhalos with $`v_m<30\mathrm{km}\mathrm{s}^1`$ only those that accreted a fraction $`f`$ of their gas before the redshift of reionization, $`z_{\mathrm{re}}`$, host luminous galaxies. Although our results are consistent with the observed abundance of dwarf satellites for a range of values of $`f`$ and $`z_{\mathrm{re}}`$, here we use the fiducial values of BKW, $`f=0.3`$ and $`z_{\mathrm{re}}=8`$. We expect that our results for disrupted satellites would be similar if we chose other $`f`$, $`z_{\mathrm{re}}`$ combinations that also match the observed (surviving) dwarf satellite population. For a given subhalo of mass $`M_a`$ and accretion redshift $`z_a`$, we use equation (2.26) of LC93 to probabilistically assign the redshift, $`z_f`$, at which the main progenitor of the subhalo reached mass $`M_f=fM_a`$ for the first time. The subhalo hosts a luminous galaxy only if $`z_fz_{\mathrm{re}}`$. We also assume that subhalos with $`v_m<v_l=10\mathrm{km}\mathrm{s}^1`$ do not host galaxies, since any gas that was initially accreted in these small systems would be unable to cool, and it should quickly boil out of the halo after reionization (e.g., Barkana & Loeb 1999). Our results do not change if we vary $`v_l`$ by $`50\%`$. The thick solid line in Figure 1 shows the average mass function of disrupted halos that once hosted luminous galaxies. The mass function of the surviving galaxies is similar to that of the disrupted galaxies, but

lower in amplitude by about a factor of two. Surviving halos are typically accreted later than tidally destroyed halos, and they are less likely to form before reionization and host a galaxy.

### 2.2. Modeling stellar tidal debris

Estimating the number of tidally stripped stars, and RR Lyrae stars in particular, requires several uncertain assumptions about gas cooling, star formation, and stellar population morphology. The following approximations are extremely simplified, and the estimated number of disrupted stars for each object is uncertain at the factor of $`23`$ level. This uncertainty is passed on to the overall amplitude of calculated stellar density distribution in the halo. Nevertheless, our statistical measures of substructure depend only on the ratio to the background density, and they are thus largely insensitive to the precise values we assume. Furthermore, we adopt the same assumptions that we used in BKW to obtain consistency with the observed dwarf satellite population, and this matching normalizes out some of the uncertainties in the overall stellar halo amplitude.

We estimate the luminosity of every disrupted luminous galaxy assuming that it has a baryonic mass $`fM_a(\mathrm{\Omega }_b/\mathrm{\Omega }_m)`$ (the mass of baryons accreted at $`z>z_{\mathrm{re}}`$). We assume that a fraction $`ϵ_{}`$ of this baryonic mass is converted to a stellar population with mass-to-light ratio $`M_{}/L_V`$. The total mass to light ratio of the subhalo is thus

$$\left(\frac{M_{\mathrm{vir}}}{L_V}\right)=f^1\left(\frac{\mathrm{\Omega }_m}{\mathrm{\Omega }_b}\right)\left(\frac{M_{}}{L_V}\right)ϵ_{}^1.$$
(1)
Adopting $`M_{}/L_V0.7`$, typical for galactic disk stars (e.g., Binney & Merrifield 1998), $`\mathrm{\Omega }_m/\mathrm{\Omega }_b7`$ (based on $`\mathrm{\Omega }_m=0.3,`$ $`h=0.7`$, and $`\mathrm{\Omega }_bh^20.02`$ from Burles & Tytler 1998), and $`ϵ_{}=0.5`$, we obtain $`(M/L_V)10f^133`$. We estimate the number of horizontal branch stars in each galaxy using $`N_{HB}=L_V/(540L_{})`$ (Preston, Shectman, & Beers 1991). The fraction of the horizontal branch stars that are RR Lyrae variables is strongly dependent on the metallicities and ages of the populations, and it will vary significantly from one object to another. For simplicity, we assume that all disrupted objects have $`N_{RR}=0.3N_{HB}`$. This fraction is high compared to local halo stars, but it is consistent with fractions observed for more distant ($`r10`$ kpc) globular clusters, likely reflecting the tendency for the outer halo to be younger (Preston et al. 1991; Brocato et al. 1996; Layden 1998).

For each disrupted luminous subhalo, we randomly assign a direction for its angular momentum vector. This direction fixes the plane of the orbit since we assume that the dark halo potential is spherical (i.e., no orbital precession). We follow the orbit from the time it was accreted at $`t_a`$ to $`t_0`$ using

$`{\displaystyle \frac{dr}{dt}}`$ $`=`$ $`\pm v(r_c)\sqrt{{\displaystyle \frac{2}{v^2(r_c)}}[\mathrm{\Phi }(r_e)\mathrm{\Phi }(r)]+1{\displaystyle \frac{ϵ^2r_c^2}{r^2}}},`$ (2)
$`{\displaystyle \frac{d\psi }{dt}}`$ $`=`$ $`{\displaystyle \frac{v(r_c)r_cϵ}{r^2}},`$ (3)
where $`r`$ is the distance to the galactic center, $`\psi `$ is the angle in the orbital plane, $`\mathrm{\Phi }(r)=4.6v_m^2\mathrm{ln}(1+x)/x`$ is the potential of the host, and the $`\pm `$ sign signifies whether the object is approaching its apocenter or pericenter, respectively. We will work in the approximation that the dynamical friction timescale is long compared to the time remaining for the orbital evolution and disruption: $`\tau _{DF}t_0t_a`$. This is a good approximation for $`90\%`$ of the disrupted halos — not surprisingly, since we have deliberately restricted our analysis to halos whose orbits have not decayed (i.e., long $`\tau _{DF}`$). In order to approximately account for cases where this approximation breaks down, we start integrating the satellite orbit at $`t=t_a`$ but set its starting radius equal to the circular radius it will have decayed to by $`t=t_0`$: $`r(t_a)=r_c(t_0)r_c`$. The initial value for the angle, $`\psi (t_a)`$, is chosen randomly, and we assume that the satellite is initially infalling (approaching the pericenter).

We assume that the satellite is tidally disrupted after the first passage of its orbit pericenter. At this time, the tidal debris will obtain an energy distribution from the encounter. Our model for the evolution of the debris along the tidal tail is motivated by numerical results of Johnston (1998), who showed that the following approximations provide good description of the positions of stripped particles in her simulations. The typical energy scale of the debris is set by the change in the host halo potential energy over the size of the tidally disrupted object,

$`\epsilon =r_t{\displaystyle \frac{d\mathrm{\Phi }}{dr}}v_m^2\left({\displaystyle \frac{r_t}{r_p}}\right),`$ (4)
where $`r_p`$ is the pericenter radius of the orbit. The last approximation is exactly true for a logarithmic potential (a singular isothermal density profile). We assume that the satellite is completely disrupted after the first passage, and that the energy of the debris is evenly distributed over the energy range $`\epsilon >dE>\epsilon `$ (Evans & Kochanek 1989). This assumption ignores the possibility of disruption over several orbits, but we find that our results are robust to the choice of distribution and do not change significantly if the energy range is altered by $`50\%`$. Using intuition gained from a circular orbit within a singular isothermal density background, we may estimate how a change in energy from $`E`$ to $`E+dE`$ affects the orbit of a particle. In this approximation, the azimuthal period, $`T_\psi `$, and radial period, $`T_r`$, depend only on the orbit energy, and they are both increased or decreased depending on the sign of the deposited energy: $`T_{\psi ,r}(E+dE)=\tau T_{\psi ,r}(E)`$, where

$`\tau =\mathrm{exp}\left({\displaystyle \frac{dE}{v_m^2}}\right).`$ (5)
This result allows us to map the orbital trajectory of the initial object with energy $`E`$ (Eqs. 4 and 5) to that of a debris particle with energy $`E+dE`$ via \[$`r(t)`$,$`\psi (t)`$\] $``$ \[$`r(t/\tau )`$,$`\psi (t/\tau )`$\].

For each RR Lyrae star in the disrupted galaxy, we assign a change in energy $`dE`$ and integrate the orbital equations to determine its position at $`z=0`$. Since the disrupted galaxy will have some finite spherical extent, we add a random offset to this calculated central orbit position. The magnitude of this offset is a Gaussian deviate with dispersion $`2`$ kpc, the typical optical radius for a dwarf galaxy (Mateo 1998), and the direction is random.

## 3. Results

Figures 2 and 3 show two realizations of the RR Lyrae distribution from disrupted satellites in sky-projected galactic coordinates. The panels in each figure correspond to the indicated radial bins in galacto-centric radius. Qualitatively, it is clear that substructures become more pronounced at larger radii. This radial trend reflects the smaller number of disrupted satellites with large apocenters and the longer periods of their orbits, which reduces the extent of debris spreading along the orbit. A comparison of the maps in Figures 2 and 3 illustrates the differences between different merger histories of the host halo. One can see that the stochastic variations in merger history at fixed host mass lead to substantial variations in the appearance and abundance of substructures.

In light of the SDSS results referred to in § 1, the most interesting predictions of the model are the radial number density profile of halo stars and the clumpiness and spatial extent of the stellar distribution. Figure 4 shows radial number density profiles of the halo RR Lyrae stars. The long dashed lines in each panel represent the power law ($`n_{}r^3`$) RR Lyrae density profile derived by Wetterer & McGraw (1996), based on their large compilation of RR Lyrae. The solid points show the profile computed

by Ivezić et al. (2000) for their sample of RR Lyrae candidates obtained from SDSS commissioning data, which covers roughly a one-degree wide, 100 square degree strip of sky. Note that the SDSS and Wetterer & McGraw profiles agree well at $`\text{ }<\text{ }35`$ kpc. At larger radii, however, the SDSS sample shows two significant deviations from the smooth $`n_{}r^3`$ profile: a “bump” in number density at $`r40`$ kpc and a sharp drop at $`r\text{ }>\text{ }50`$ kpc. As noted by Ivezić et al., this structure in the radial profile likely indicates significant clumpiness of the stellar halo at these galactocentric radii, and the bump in particular is associated with an identifiable coherent structure containing $`70`$ RR Lyrae within the observed region. The model predictions are shown in Figure 4 by thick and thin solid lines. The thick solid lines represent the computed RR Lyrae number density profile averaged over all merging history realizations and the full sky. The error bars in the upper left panel show the dispersion from realization to realization around this average, demonstrating that stochastic variations in merger histories lead to a factor of $`2`$ rms variation in the overall normalization of the predicted halo density profile. In the remaining three panels, the thin solid lines show examples of density profiles derived from a single host halo realization, viewed through three randomly chosen strips similar in solid angle and geometry to the strips used to derive the SDSS sample. The first remarkable feature of Figure 4 is the agreement of the predicted average profile with the slope and amplitude found by Wetterer & McGraw (1996). As discussed in §2.2, the predicted amplitude is uncertain by a factor of several, and even the statistical fluctuations from one galaxy to another are significant, so the degree of agreement must be somewhat fortuitous. However, it is worth noting that we did not adjust any parameters to fit the observed halo profile but chose “best guess” values based on other considerations — in particular, the requirement of matching the observed dwarf satellite population. Figure 4 suggests that disruption of accreted satellites can produce not just substructure in the stellar halo but the entire stellar halo itself. If we have overestimated the number of RR Lyrae per unit dark matter mass, then there is room for another physical mechanism that creates a smooth underlying halo, but it seems that no such additional mechanism is necessary. The second remarkable feature of Figure 4 is the jaggedness of the observed and predicted profiles of individual strips, which becomes especially pronounced at radii $`r40`$ kpc. These large fluctuations reflect the substructure visible in Figures 2 and 3. The steep drop in the SDSS counts between 60 and 70 kpc suggests detection of an “edge” of the stellar halo (Ivezic et al. 2000). However, the second of our numerical realizations shows an equally sharp edge, even though the average halo profile is smooth. Our model predicts a gradual steepening of the halo profile at $`r60`$ kpc, but although surveys in small solid angles should show large count fluctuations, the profile averaged over the full sky should not cut off sharply. Figures 5 and 6 present a different view of the structure associated with disrupted satellites, in a form more comparable to the plots of Yanny et al. (2000). Each figure shows a one degree wide, randomly oriented, great circle slice through a realization of the RR Lyrae distribution. The two figures show stellar halos for different Monte Carlo accretion histories, one with a total of about $`60`$ tidally destroyed galaxies (Figure 5) and one with a more quiescent accretion history and about $`20`$ tidally destroyed galaxies (Figure 6). This range roughly covers the scatter in the number of disruption events expected from galaxy to galaxy. The three concentric circles indicate galactocentric radii of 25, 50, and 75 kpc. In both figures, structures associated with the individual disrupted objects become more easily identifiable at larger radii. To make our predictions more quantitative, we present two simple statistical measures of the halo clumpiness. Figure 7 shows the probability distribution of model RR Lyrae counts in solid angle cells of different sizes and for different ranges of galactocentric radii. The cells are roughly square on the sky; they are defined by dividing the sky map for each realization into patches of a given size. In order to take out the uncertainty in the overall amplitude of the density distribution, and to factor out the variation in the overall amplitude from one halo realization to another, the counts are presented in units of the average expected number, $`N`$, of RR Lyrae in patches of the chosen size for each realization. Figure 7 shows that the amplitude of patch-to-patch count fluctuations is higher for smaller solid angles and for larger galactocentric radii. Figure 8 shows the rms dispersion in RR Lyrae counts, $`\sigma (N/N)`$, as a function of galactocentric radius for cells of different solid angle and geometry. The solid points/line show counts in square cells, while open points represent the counts in one degree wide strips. The error bars for the solid points show the standard deviation of the dispersion for different merger history realizations; the error bars are similar for the open points. These error bars reinforce the point made in our illustrative figures above, namely that the degree of surviving substructure is quite variable from one realization of the stellar halo to another.The dashed lines show the expected amplitude of Poisson fluctuations based on the average number of stars expected within the given solid angle and radial bin. Note that the predicted fluctuations are always larger than Poisson fluctuations, especially at large radii, because they are dominated by fluctuations in the number of debris streams rather than $`\sqrt{N}`$ fluctuations in the number of RR Lyrae. Again, the fluctuations are larger for cells of smaller solid angle and for larger radii. They are also larger for square patches than for narrow strips of the same solid angle, reflecting the fact that a narrow strip is less likely to enclose an entire disrupted object and instead encloses fragments of multiple debris streams.

## 4. Discussion and conclusions

We have presented a model in which stellar halos of Milky Way type galaxies are built via the accretion and tidal disruption of satellite galaxies. The model is based on the CDM structure formation scenario, with the crucial assumption that only a fraction of dark matter halos with circular velocities $`\text{ }<\text{ }30\mathrm{k}\mathrm{m}\mathrm{s}^1`$ are luminous and host a sizable stellar system222All DM halos with circular velocities $`>30\mathrm{k}\mathrm{m}\mathrm{s}^1`$ are assumed to host a luminous galaxy.. These luminous halos are those that collapsed and accreted a substantial fraction of their mass prior to the epoch of reionization, $`z_{\mathrm{re}}`$. Only the accretion and tidal disruption of luminous dwarf galaxies contribute to the build up of the stellar halo in our model. The model predicts an average density profile for the stellar halo of $`n_{}r^\alpha `$, with $`\alpha 3`$ in the range $`r1050`$ kpc, in good agreement with observations. The density profiles of individual realizations of the host halo merger histories, however, exhibit a substantial scatter around this average, with the power law varying from $`\alpha 2.53.5`$. The reason why the stellar halo density profile is steeper than that of the dark matter ($`\alpha _{DM}2`$, which is indeed typical of that of the surviving subhalos in our model), is that central satellites are more likely to be destroyed than those at large radii, since satellites far from the galactic center must have extremely eccentric orbits in order to pass close enough to be tidally destroyed. The average predicted profile falls off more steeply at $`r50`$ kpc, but it does not have a sharp outer boundary. The amplitude of the mean profile depends on several uncertain factors, such as the mean stellar luminosity to dark mass ratio of the disrupted dwarfs and the mean number of RR Lyrae per solar luminosity. There are also significant (factor of two) statistical fluctuations in the amplitude from one realization to another. The amplitude of the predicted RR Lyrae profile is therefore uncertain by a factor of a few. Nonetheless, with “best guess” parameter values chosen on the basis of other considerations, the mean RR Lyrae profile agrees very well with that determined by Wetterer & McGraw (1996), in slope and amplitude. It therefore appears that disruption of dwarf satellites is a plausible mechanism for producing the entire stellar halo within the CDM framework, though it could also co-exist with some other mechanism. The uncertainty in our predicted normalization is reduced by the fact that we require the model to self-consistently reproduce the velocity function of observed dwarf satellites, a point that we will return to shortly. The main qualitative prediction of the model is the presence of significant clumpiness in the outer regions ($`\text{ }>\text{ }30`$ kpc) of the Galaxy’s stellar halo. This clumpiness is due to the surviving tidal debris of dozens of satellite galaxies disrupted during evolution of their host. In the inner regions of the stellar halo ($`r\text{ }<\text{ }30`$ kpc), the density distribution is relatively smooth. At larger radii, however, the clumpiness of the stellar halo manifests itself when viewed through fixed solid angle patches in the sky. For typical modeled stellar halos, RR Lyrae profiles of the type observed by Ivezić et al. (2000) in the SDSS commissioning data (with a coherent structure at $`r50`$ kpc) are not uncommon. We have quantified our predictions by presenting some statistical measures of the “clumpiness” in our modeled stellar halos. First, we measured the probability distribution (see Fig. 7) of RR Lyrae counts, $`N/N`$, in solid angle cells of different sizes and for different galactocentric radial bins (here $`N`$ is the average number of stars expected in a cell). We also presented the rms width of this probability distribution as a function of galactocentric radius for solid angle cells of different geometries (Fig. 8). These statistics show that the variance in the stellar counts (i.e., clumpiness of the stellar halo) increases with increasing galactocentric radius and decreasing solid angle. Although current observational data sets are not sufficiently large to derive similar statistics for comparison, future samples of RR Lyrae and A-stars from the SDSS and 2dF surveys should make such comparisons possible. For cell sizes of 36 deg2 or larger, predicted fluctuations in RR Lyrae counts are much larger than Poisson fluctuations, at all galactocentric radii considered. These quantitative predictions depend on some of the model assumptions. For instance, one of the implicit assumptions in our analysis is that the dark matter distribution in the host halo is nearly spherically symmetric. In particular, we do not include the possible precession of satellite orbits. In the case of the significantly oblate halo, such precession can erase signatures of the tidal tails (at least in configuration space), so it would reduce the predicted clumpiness of the outer halo (e.g., Ibata et al. 2000a). However, the observed narrowness of the tidal tail of the Sagittarius dwarf galaxy implies a nearly spherical halo potential at $`r\text{ }<\text{ }60`$ kpc (Ibata et al. 2000a). A nearly spherical mass distribution in the inner region of a galaxy-mass halo is consistent with predictions of CDM models (Bullock et al. 2000b). We have also neglected the effects of the disk component on the background potential. Including this non-spherical central component would induce precession in the tidal orbits (Helmi & White 1999) and smear out residual structure at small galactic radii. But since this effect would only be important in the central regions ($`r\text{ }<\text{ }20`$ kpc), the net effect would be to strengthen the trend of increased variance in star counts with radius. The assumption that affects our predictions the most, however, is that only a small fraction of disrupted subhalos host a stellar system and thereby contribute to the build up of the stellar halo. In other words, the predicted properties of the stellar halo depend crucially on the way that we have solved the dwarf satellite problem. This problem, namely that the predicted number of dark halos with circular velocities $`\text{ }<\text{ }30\mathrm{km}\mathrm{s}^1`$ is much larger than the observed number of dwarf satellites within the virialized dark halos of the Milky Way and M31 (Kauffmann et al. 1993; Klypin et al. 1999a; Moore et al. 1999), is one of the few outstanding problems of the conventional inflationary CDM model, which, with the inclusion of a cosmological constant, accounts well for a wide variety of other observational data. Other proposed solutions to the dwarf satellite problem include modifying the inflationary fluctuation spectrum (Kamionkowski & Liddle 2000) or modifying the properties of dark matter by making it warm (WDM; e.g., Hogan & Dalcanton 2000) or self-interacting (SIDM; e.g., Spergel & Steinhardt 2000). In the SIDM model, halos collapse and accrete mass in a similar manner to conventional CDM halos. However, the number of surviving subhalos within a Milky Way mass halo is reduced because lower concentration of the SIDM halos makes them easier to disrupt and because interactions lead to “ram pressure” stripping of the dark halos. Relative to our reionization solution, the SIDM solution seems to predict a more massive and more extended stellar halo built from a larger number of tidal streams, since the abundance of dwarf satellites is reduced by a higher efficiency of satellite disruption rather than suppression of star formation in low-mass systems. Indeed, the SIDM model seems at some risk of overproducing the stellar halo. However, if the predicted stellar halo were normalized to the mean stellar density of the observed one, the SIDM model would probably predict less substructure than the model presented here because of the higher density of independent tidal streams. The predictions of WDM and broken scale-invariance models are less clear, since they can reduce dwarf satellite numbers both by suppressing their formation in the first place (White & Croft 2000) and by making them less concentrated and therefore easier to disrupt (Colín et al. 2000). If the second effect dominates, the predictions might be closer to those of SIDM; if the first effect dominates, they might be closer to those of the reionization model. At present, we lack quantitative predictions on this matter from SIDM and WDM, and we lack detailed observational constraints, so we cannot draw conclusions about which model fares best. However, it appears that studies of the mass, radial profile, and substructure of the stellar halo might provide useful constraints on these models in the near future. It is interesting that all of the substructures detected in the nearly all-sky Carbon star survey and the SDSS commissioning data could be produced by a single tidal stream of the Sagittarius dwarf galaxy (Ibata et al. 2000b). In particular, the observed excess in the RR Lyrae number density profile at $`r45`$ kpc is caused by a clump of stars near the apocenter of the Sagittarius orbit. If most of the RR Lyrae in the surveyed strip belong to the Sagittarius stream, then no stars are expected beyond the apocenter of the Sagittarius orbit, and this would naturally explain the drop in the RR Lyrae number density at $`r\text{ }>\text{ }60`$ kpc detected in the SDSS data. However, the presence of only a single stream in a large volume of sky would be both puzzling and intriguing, since the CDM models considered in this paper predict at least $`20`$ (and typically more) tidal streams in the halo of a Milky Way-size galaxy. The SIDM and WDM stellar halos should be built from even larger number of tidal streams. Carbon stars are relatively young and rare, and therefore older tidal streams may not be revealed by their distribution. We predict, however, that many more substructures not associated with the Sagittarius tidal stream should be detected in the future as the SDSS covers a larger area of the sky. Absence of such detections would spell serious trouble for CDM models, and possibly even more so for the variants of CDM that we have discussed. Several past theoretical studies of tidal stripping and disruption of galactic satellites have shown that tidal streams can be used as powerful probes of both the present potential of the Milky Way (Johnston et al. 1999; Ibata et al. 2000a) and its accretion history (e.g., Helmi & White 1999; Helmi et al. 1999; Helmi & de Zeeuw 2000). For example, Johnston et al. (1999) show that parameters characterizing the mass distribution in the Milky Way halo can be determined with an accuracy of a few percent using a single tidal stream if accurate phase-space information is available for as few as 100 stream stars. Helmi et al. (1999) show how phase-space information can be used to recover fossil remnants of the satellites accreted and disrupted early in the Milky Way evolution. These techniques can be used to construct a detailed formation history of the Galaxy. In this study, we have focused on the spatial distribution of stars in the Milky Way halo. The radial density behavior of stars may not prove to be as sensitive a probe of the galactic potential as the accurate phase-space mapping of tidal streams (Johnston et al. 1999), but it should be very useful in recovering details of the Milky Way accretion history. In some respects, the number of surviving tidal streams or the degree of clumpiness of the stellar distribution can provide better constraints on the variants of the CDM scenario than the abundance of galactic satellites. An advantage here is that the spatial distribution of halo stars should be possible to map in the very near future, when large samples of halo stars from the SDSS and other surveys become available. Accurate, large-scale, phase-space mapping of tidal streams, on the other hand, will become possible only after launch of the next generation astrometric space missions (i.e., at the end of the decade). On the theoretical side, our model and predictions can be improved upon by combining mass accretion histories typical for the galaxy-size halos formed in CDM models with more sophisticated numerical models of orbital evolution of satellites and their tidal debris.

We thank Amina Helmi and Kathryn Johnston for useful discussion and suggestions. This work was supported in part by NASA LTSA grant NAG5-3525 and NSF grant AST-9802568. Support for A.V.K. was provided by NASA through Hubble Fellowship grant HF-01121.01-99A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.