text:
# Vortices in the presence of a nonmagnetic atom impurity in 2D XY ferromagnets

## Abstract

Using a model of nonmagnetic impurity potential, we have examined the behavior of planar vortex solutions in the classical two-dimensional XY ferromagnets in the presence of a spin vacancy localized out of the vortex core. Our results show that a spinless atom impurity gives rise to an effective potential that repels the vortex structure.

The easy-plane Heisenberg ferromagnet in two dimensions and continuum limit supports non-linear pseudo-particles with a vortex structure. These excitations are of paramount importance in the understanding of static and dynamical properties of magnetism. For example, the vortex unbinding is responsible for a phase transition known as Kosterlitz-Thouless transition 1 . Besides, these may be responsible for a central peak in the dynamical correlation function 2 ; 3 ; 4 ; 5  observed in Monte-Carlo simulations 6 ; 7  and experiments 8 . The simplest vortex configuration, referred as the planar vortex, occurs when the anisotropy is large, resulting in spin confinement to the lattice plane 9 ; 10 .

The interaction of vortices with spatial inhomogeneities is of considerable importance from both the purely theoretical and applied point of view. Impurities and/or defects are present even in the purest of material samples and their effect on the motion or structure of nonlinear excitations must be considered when the dynamics or configurations of such solutions are important in the problem at hand. Recently, Zaspel, McKennan and Snaric 11  investigated, using the discrete lattice, the instability of planar vortices and concluded that they will be stable at a large range of anisotropies if there is a nonmagnetic impurity such as Cd or Zn at the center of the vortex. In this paper we study, using the continuum approximation, the interaction between a planar vortex and a nonmagnetic impurity localized out of the vortex center. To this end, we start defining the classical XY ferromagnet model, which is given by the following Hamiltonian:

$`H=J{\displaystyle \underset{m,n}{}}(S_m^xS_n^x+S_m^yS_n^y),`$ (1)
Where $`J`$ is a coupling constant, the classical spin vector has three components $`\stackrel{}{S}=(S^x,S^y,S^z)`$ and the summation is taking over the nearest-neighbor square lattice sites. This model is one of the most studied in statistical physics and has been found to describe a wide variety of systems with complex scalar order parameters, including superconductor films, Josephson junction arrays, and superfluid He4 films. The choice of the XY model is arbitrary for our purpose, since the results can be used, without modifications, to any other model with XY symmetry such as the classical easy-plane ferromagnetic model. It is convenient to parametrize the spin field in terms of spherical coordinates as follows: $`\stackrel{}{S}=S(\mathrm{cos}\theta \mathrm{cos}\varphi ,\mathrm{cos}\theta \mathrm{sin}\varphi ,\mathrm{sin}\theta ).`$ By a straightforward generalization of arguments used to obtain the continuum limit of Heisenberg Hamiltonians in the case of one dimension 12 , we can write the continuum version of the Hamiltonian (1) as

$$H=\frac{J}{2}d^2r\left[\frac{m^2(\stackrel{}{}m)^2}{1m^2}+(1m^2)(\stackrel{}{}\mathrm{\Phi })^2+\frac{4}{a_0^2}m^2\right],$$
(2)
where $`m=\mathrm{sin}\theta `$, $`a_0`$ is the lattice constant and we have taken $`S^2=1`$. One can obtain the motion equations, $`m/t=\delta H/\delta \mathrm{\Phi },\mathrm{\Phi }/t=\delta H/\delta m,`$ for this field theory in the usual way using the pair of canonically conjugated variables $`m`$ and $`\mathrm{\Phi }`$.

Now, consider that the system contains a nonmagnetic impurity concentration in the plane. For simplicity, we consider only one nonmagnetic atom present. If we remove a spin from the lattice, the nearest neighbors of that spin will have a coordination number of three, instead of bulk-spin coordination number of four. Therefore, such boundary spins would have larger fluctuations than the bulk spins and it is conceivable that the non-linear configurations such as a vortex, would preferentially nucleate around this vacancy. In fact, in this circumstance, the vortex energy is lowered, since the nonmagnetic impurity at the vortex center will remove the nearest-neighbor exchange bonds at the impurity in a radially symmetric way without modifying the symmetric configuration of the vortex, while its energy in the region without the impurity remains the same. However, the vortex energy increases logarithmically with the system size $`L`$, and in an infinitely extended system, this energy would diverges as $`L\mathrm{}`$ so that we should not expect that a single vortex could nucleate around the spin vacancy. In fact, vortices are always created in pairs of vortex-antivortex having finite energy and separated by a few lattice constants. Then, if one member of the pair nucleates around the impurity, the other member will be near it. Since a vortex pair does not have a cylindrical symmetry, the energy of this system does not necessarily decrease, although there are less nearest-neighbor exchange bounds, because the spin vacancy may deform the pair configuration, increasing its energy. Nevertheless, the interest of this paper is to study the behavior of a vortex in the system, that is not nucleated at the position of the vacancy. In this case things must change considerably because the spin vacancy may also deform the single vortex configuration. To take into account the nonmagnetic impurities we consider the following modified XY ferromagnetic Hamiltonian in the continuum limit

$`H_I={\displaystyle \frac{J}{2}}{\displaystyle d^2r\left[\frac{m^2(\stackrel{}{}m)^2}{1m^2}+(1m^2)(\stackrel{}{}\mathrm{\Phi })^2+\frac{4}{a_0^2}m^2\right]}`$
$`\times V(\stackrel{}{r}),`$ (3)
where $`V(\stackrel{}{r})`$ is a nonmagnetic impurity potential given by

$`V(\stackrel{}{r})=\{\begin{array}{cc}\hfill 1& \text{if}|\stackrel{}{r}\stackrel{}{r}_0|b,\\ \hfill 0& \text{if}|\stackrel{}{r}\stackrel{}{r}_0|<b.\end{array}`$ (6)
Here the impurity is centralized at the point $`\stackrel{}{r}_0`$ and has the form of a circle with diameter equal $`2b`$. There is a circular region in the plane, around the point $`\stackrel{}{r}_0`$, without any magnetic interaction.

Substituting Eq. (6) into the equations of motion we get

$`{\displaystyle \frac{1}{J}}{\displaystyle \frac{\theta }{t}}=\mathrm{cos}\theta V(\stackrel{}{r})^2\mathrm{\Phi }2\mathrm{sin}\theta V(\stackrel{}{r})\stackrel{}{}\theta \stackrel{}{}\mathrm{\Phi }`$
$`+\mathrm{cos}\theta \stackrel{}{}V(\stackrel{}{r})\stackrel{}{}\mathrm{\Phi },`$ (7)
$`{\displaystyle \frac{1}{J}}{\displaystyle \frac{\mathrm{\Phi }}{t}}=\mathrm{tan}\theta \mathrm{sin}\theta V(\stackrel{}{r})^2\theta \mathrm{sin}\theta V(\stackrel{}{r})(\stackrel{}{}\theta )^2`$
$`+\mathrm{sin}\theta V(\stackrel{}{r})[4/a_0^2(\stackrel{}{}\mathrm{\Phi })^2]\mathrm{tan}^2\theta \mathrm{sin}\theta \stackrel{}{}V(\stackrel{}{r})\stackrel{}{}\theta .`$ (8)
Our interest is in planar and static solutions in the presence of this nonmagnetic impurity potential. Hence, we take $`\theta /t=\mathrm{\Phi }/t=0`$ and $`m=\mathrm{sin}\theta =0`$ in Eqs. (7) and (8), obtaining only one and simpler equation to be solved

$`V(\stackrel{}{r})^2\mathrm{\Phi }=\stackrel{}{}V(\stackrel{}{r})\stackrel{}{}\mathrm{\Phi }.`$ (9)
One point to note in this equation is its dependence on the spin field around the position of the impurity. If the vacancy is localized in a region where the spin configuration consists of aligned spins like a domain with all spins aligned along some direction $`(|\stackrel{}{}\mathrm{\Phi }|0)`$, the spin field practically does not feel the presence of the impurity. However, an impurity placed in a region where the spin directions vary considerably $`(|\stackrel{}{}\mathrm{\Phi }|>>1/a_0)`$ may have a strong coupling with the spin field and may modify the initial spin configuration for large distances.

In polar coordinates the vectors $`\stackrel{}{r}`$ and $`\stackrel{}{r}_0`$ are written as $`(r,\varphi ),`$ and $`(r_0,\varphi _0)`$ respectively. To solve Eq.(7), we note first that the gradient of the impurity potential can be expressed as

$`\stackrel{}{}V(\stackrel{}{r})=a_0[\widehat{r}\mathrm{cos}(\alpha |\varphi \varphi _0|)+\widehat{\varphi }\mathrm{sin}(\alpha |\varphi \varphi _0|)]`$
$`\times \delta (\stackrel{}{r}\stackrel{}{r}_0\stackrel{}{b}),`$ (10)
where $`\delta `$ is the Dirac delta function and $`\alpha `$ is the angle that the vector $`\stackrel{}{b}`$, with origin at the point $`\stackrel{}{r}_0`$ and end at a point on the circumference of the potential, makes with the vector $`\stackrel{}{r}_0`$. As we are interested in a local impurity with atomic dimensions, we make $`b0`$ in the continuum limit (to be more precise, we should make $`ba_0`$) indicating that the impurity is an atom (such as Zn, Mg or Cd for example). In this case, $`\stackrel{}{r}\stackrel{}{r}_0,\varphi \varphi _0,`$ and we rewrite Eq. (10) as

$`\stackrel{}{}V(\stackrel{}{r})a_0[\widehat{r}\mathrm{cos}(\alpha )+\widehat{\varphi }\mathrm{sin}(\alpha )]\delta (\stackrel{}{r}\stackrel{}{r}_0),`$ (11)
where we can interpret $`\mathrm{cos}(\alpha )`$ and $`\mathrm{sin}(\alpha )`$ as anisotropic coupling constants. This coupling depends on the direction one looks, if the observer center is placed on the impurity position.

Considering Eq. (9) with $`V(\stackrel{}{r})=1`$ at the left side (this fails only at the point $`\stackrel{}{r}_0`$, since the impurity is local) and supposing that the vortex structure is modified by the presence of the nonmagnetic impurity, we write $`\mathrm{\Phi }=\mathrm{\Phi }_0+\mathrm{\Phi }_1`$, where $`\mathrm{\Phi }_0=\mathrm{arctan}(y/x)`$ is the traditional single vortex solution for a vortex with center localized at the origin and $`\mathrm{\Phi }_1`$ is the deformation caused by the spinless impurity localized at $`\stackrel{}{r}_0`$. Thus, Eq. (9) with the above considerations can be written as

$`^2(\mathrm{\Phi }_0+\mathrm{\Phi }_1)=a_0\stackrel{}{}(\mathrm{\Phi }_0+\mathrm{\Phi }_1)`$
$`[\widehat{r}\mathrm{cos}(\alpha )+\widehat{\varphi }\mathrm{sin}(\alpha )]\delta (\stackrel{}{r}\stackrel{}{r}_0)`$ (12)
Using the fact that $`^2\mathrm{\Phi }_0=0`$ and taking $`\stackrel{}{}(\mathrm{\Phi }_0+\mathrm{\Phi }_1)\stackrel{}{}\mathrm{\Phi }_0=\frac{1}{r}\widehat{\varphi }`$ near the point $`\stackrel{}{r}_0`$, Eq. (12) can then be approximated by

$`^2\mathrm{\Phi }_1={\displaystyle \frac{a_0}{r_0}}\mathrm{sin}(\alpha )\delta (\stackrel{}{r}\stackrel{}{r}_0),`$ (13)
or

$`^2\left[{\displaystyle \frac{2\pi r_0\mathrm{\Phi }_1}{a_0\mathrm{sin}(\alpha )}}\right]=2\pi \delta (\stackrel{}{r}\stackrel{}{r}_0).`$ (14)
This is easily solved using the fact that in two dimensions, $`^2\mathrm{ln}(r)=2\pi \delta (\stackrel{}{r})`$. We get

$`\mathrm{\Phi }_1(\stackrel{}{r})={\displaystyle \frac{a_0\mathrm{sin}(\alpha )}{2\pi r_0}}\mathrm{ln}\left({\displaystyle \frac{|\stackrel{}{r}\stackrel{}{r}_0|}{a_0}}\right).`$ (15)
Writing the anisotropic coupling constant along the $`\alpha `$-direction in terms of $`r`$, $`\varphi `$, the vortex structure with its center at the origin in the presence of a nonmagnetic impurity localized at $`\stackrel{}{r}_0`$ is given by

$`\mathrm{\Phi }=\mathrm{arctan}(y/x){\displaystyle \frac{a_0}{2\pi r_0}}{\displaystyle \frac{r\mathrm{sin}(\varphi \varphi _0)}{|\stackrel{}{r}\stackrel{}{r}_0|}}\mathrm{ln}\left({\displaystyle \frac{|\stackrel{}{r}\stackrel{}{r}_0|}{a_0}}\right).`$ (16)
The configuration of this deformed vortex is shown in figures 1 and 2. Although the continuum theory cannot be applied near the vortex core, in Fig. 1 we have considered the impurity one lattice sapcing from the vortex center just to emphasize the vortex deformation as the vortex core approaches the impurity.We notice that if $`r_0`$ is large (the vortex center is far away from the spin vacancy) the vortex practically keeps the same original form, but for small $`r_0`$ the vortex configuration suffers a severe modification, mainly in the region in which the impurity is located. This is due to the fact that the gradient of the spin field is small in the region of the impurity if it is far way from the vortex center and large in the region of the impurity if it is near the vortex center. When animpurity is near the vortex core, Eq. (16) implies (as it can also be seen partially in figure 1) that a large domain with all spins aligned along the direction perpendicular to $`\stackrel{}{r}_0`$ will be formed in a region located after the impurity position $`\stackrel{}{r}_0`$.

In order to calculate the energy of this planar solution we consider the Hamiltonian (3) with $`m=0`$, obtaining $`E_I=(\stackrel{}{}\mathrm{\Phi })^2V(\stackrel{}{r})d^2r`$. As we have seen, the field $`\mathrm{\Phi }`$ describes a single vortex at the origin in the presence of an impurity at distance $`r_0`$ away. The effective potential experienced between the two defects (one defect in the spin field and the other in the lattice structure) is defined as

$`U_{eff}(r_0)=E_IE_\nu ,`$ (17)
where $`E_\nu =\pi J\mathrm{ln}(L/a_0)`$ is the energy of a single vortex in the absence of impurities. Making suitable approximations, we find that such effective potential results in a repulsive central interaction with a dominant term given by

$`U_{eff}(r_0){\displaystyle \frac{a_0^2E_\nu ^3}{24\pi ^4J^2}}{\displaystyle \frac{1}{r_0^2}}.`$ (18)
We see, therefore, that the presence of a nonmagnetic impurity increases the vortex energy as the distance between the impurity and the vortex decreases. In a ferromagnet with a size of the order of $`L10^8a_0`$ ( a few centimeters), a spinless atom impurity situated about $`2a_0`$ from the vortex core would increase the vortex energy about 36%. Note that the effective potential barrier becomes infinity as $`r_00`$ and it is energetically favorable that vortices and impurities become far apart. But, if the calculations were taken considering that the spin vacancy is localized at the vortex center, we would have $`\stackrel{}{}V(\stackrel{}{r})=a_0\delta (\stackrel{}{r})\widehat{r}`$ and near the vortex core $`\stackrel{}{}\mathrm{\Phi }(1/a_0)\widehat{\varphi }`$, leading to $`V(\stackrel{}{r})^2\mathrm{\Phi }=0`$. As a consequence, in the region without the spinless impurity, where $`^2\mathrm{\Phi }=0`$, one gets the same typical solutions and the vortex structure does not suffer any alteration. Hence, the only one effect of a central nonmagnetic impurity, is to make the vortex energy to decrease, because of the nonexitence of nearest-neighbor exchange bonds at the impurity. Nevertheless, as we suggested earlier, a single vortex with infinite energy may not nucleate by itself around the impurity, since these are created in pairs. Besides, Eq. (18) shows that an infinite potential barrier has to be exceeded by the vortex core in order that it might reach the nonmagnetic impurity and the minimum of energy. Then, one should not expect to find a single vortex with a spin vacancy localized in its center.

In Summary, vortices prefer to stay far way from nonmagnetic impurities and hence, the spin dynamics must be affected by these lattice defects. Our calculations could also be taken for two-dimensional easy-plane antiferromagnets. It would be carried out in essentially the same way, leading to similar results. The structure and motion of vortices in two-dimensional magnets may be driven by the presence of spinless impurities due to the repulsive effective potential. Since the dynamical structure factor is the Fourier transform of the vortex spatial and temporal configuration, we expect that nonmagnetic impurities may cause changes in the central peak 2 ; 3 ; 4 ; 5  and also in the electron paramagnetic resonance linewidth 13 , which must be seen in neutron scattering and resonance experiments. However, much work has to be done in order to see these effects. Besides, since the vortex energy is modified, the Kosterlitz-Thouless temperature may also be affected by the presence of impurities. In fact this theory holds also for temperatures below the Kosterlitz-Thouless temperature $`T_{KT}`$, where vortices are bound in pairs. However, the problem of vortex pairs interacting with nonmagnetic impurities and its influences on $`T_{KT}`$ will be treated in a future paper.We also suggest that the above calculations may have some relevance to high-Tc superconductors, because a common feature of all high-transition-temperature cuprates is the proximity between antiferromagnetic and d-wave superconducting phases controlled by the doping. The effect of impurities on superconductors has been of theoretical and experimental interest even in its own right for a long time. Recent nuclear magnetic resonance measurements have shown that when a $`Cu^{2+}`$ in the $`CuO`$ plane is substituted by a strong nonmagnetic impurity, such as $`Zn^{2+}`$, an effective magnetic moment can be induced on the $`Cu`$ sites around the impurity site 15 ; 16 . The physical picture implied by these experiments is that antiferromagnetic correlations are enhanced, not destroyed, around impurities in these cuprates.

###### Acknowledgements.

This work was partially supported by CNPq and FAPEMIG (Brazil).


format: text

text:
# Ring Nebulae around Massive Stars throughout the HR Diagram

## 1. Introduction

Since Johnson & Hogg (1965) reported the first three ring nebulae around Wolf-Rayet (WR) stars, nebulae of various shapes, sizes, and ionization conditions have been observed around massive stars of different spectral types, such as luminous blue variables (LBVs), blue supergiants (BSGs), and red supergiants (RSGs). As these various spectral types in the upper part of the Hertzsprung-Russell (HR) diagram are strung together by the evolutionary tracks of massive stars, their surrounding nebulae must be evolutionarily related. This relationship was illustrated clearly by García-Segura, Langer, & Mac Low (1996) and García-Segura, Mac Low, & Langer (1996) in their hydrodynamic simulations of the formation of WR ring nebulae taking into account the evolution and mass loss history of the central stars.

In this review, I will cover three topics. First, I will present a gallery of ring nebulae around massive stars throughout the HR diagram and compare them to García-Segura et al.’s framework of hydrodynamic evolution of circumstellar gas around massive stars. Second, I will report an X-ray view of the hot gas in the interiors of bubbles blown by massive stars, using ROSAT, ASCA, Chandra, and XMM-Newton observations. Finally, I will compare the circumstellar nebulae observed around young supernovae (SNe) to those around massive stars and use the nebular properties to diagnose the spectral types of SN progenitors.

## 2. Gallery of Ring Nebulae around Massive Stars

The hydrodynamic models of ring nebulae around massive stars by García-Segura et al. considered two possible stellar evolution tracks from the main sequence (MS) to the WR phase: MS O $``$ LBV $``$ WR, for a 60 M star; MS O $``$ RSG $``$ WR, for a 35 M star. Stars at these different evolutionary stages expel stellar material in the following ways: MS O stars possess tenuous, fast winds with terminal velocities of 1,000–2,000 km s-1; RSGs have copious, slow winds with wind velocities typically a few 10’s km s-1; LBVs have outbursts or winds at modest velocities that depend on the stellar luminosities; WR stars have the most powerful stellar winds with both large mass loss rates, $`10^5`$ M yr-1, and high wind velocity, up to 3,000 km s-1. These stellar winds interact with the ambient medium and form ring nebulae with distinct characteristics for each type of central stars.

### 2.1. Main Sequence O Stars – Interstellar Bubbles

The fast stellar wind of a MS O star sweeps up the ambient interstellar medium (ISM) to form an interstellar bubble (Weaver et al. 1977), which consists of a dense shell of interstellar material. Intuitively, we would expect around most O stars an interstellar bubble similar to the Bubble Nebula (NGC 7635) to be visible; however, hardly any O stars in H ii regions have ring nebulae, suggesting that these interstellar bubbles are rare. This puzzle is recently solved by observations of the young H ii region N11B in the Large Magellanic Cloud (LMC). While no ring nebulae can be identified morphologically, long-slit high-dispersion spectra have revealed expanding shells around single or groups of O stars (see Fig. 1a,b; Nazé et al. 2001). Their expansion velocities, typically 10–15 km s-1, are not much higher than the isothermal sound velocity of the 104 K ionized gas in the H ii region, and consequently cannot produce strong shock compression. Lacking large density enhancements, these slowly expanding interstellar bubbles do not have visible ring nebulae.

Interstellar bubbles blown by MS O stars may become morphologically identifiable at later evolutionary stages when the central stars can no longer ionize the bubble shells and the surrounding ISM. A 10 km s-1 expanding shell in a 100 K neutral H i medium (isothermal sound velocity $``$ 1 km s-1) is highly supersonic, thus can generate strong compression and produce an identifiable shell morphology. Such recombined interstellar bubbles have been observed as H i shells around WR stars (Cappa 2002, and references therein). A beautiful H i shell around WR 134 in our Galaxy reported by Gervais & St-Louis (1999) is reproduced in Fig. 1c.

### 2.2. Red Supergiants – Circumstellar Nebulae

Evolved massive stars at the RSG phase lose mass via slow winds and form circumstellar nebulae. Because of the high dust content, these circumstellar nebulae are best observed through the optical continuum scattered by dust or the thermal IR continuum emitted by dust. This is illustrated by VY CMa (M5e Ia) in Fig. 2. The optical reflection nebula of VY CMa, visible to a radius of 9000 AU (for a distance of 1.5 kpc), is asymmetric about the central star. The presence of bright arcs indicates that localized ejection of stellar material may have occurred on the stellar surface (Smith et al. 2001). More details on mass loss and circumstellar nebulae of RSGs are reviewed by Humphreys (2002).

### 2.3. Luminous Blue Variables – Circumstellar Bubbles

The name LBV has been used for both the extreme case of $`\eta `$ Car and the less luminous S Doradus variables. The Ofpe/WN9 stars have also been loosely called LBVs. The circumstellar nebulae of LBVs have been reviewed by Nota et al. (1995). Using the heavy element abundances of LBV nebulae, Smith et al. (1998) suggest that LBVs have gone through a brief RSG phase, while Lamers et al. (2001) demonstrate that the nebulae are ejected during the blue supergiant (BSG) phase and that the LBVs have never gone through a RSG phase.

Most circumstellar nebulae around LBVs are small, $`<`$2 pc in diameter, expanding shells with Vexp of few 10’s km s-1 (Nota et al. 1995; Weis 2002). Many LBV nebulae in the Galaxy and in the Magellanic Clouds have been studied recently by, e.g., Nota et al. (1996, 1997), Pasquali et al. (1997, 1999), and Weis (2002, and references therein). A large compilation of images of LBV nebulae were presented in the poster by Weis in this meeting. Fig. 3 shows the intriguing runaway LBV S119 and its nebula in the LMC; the interstellar absorption lines in the FUSE spectrum of S119 have established unambiguously that it is in the LMC (Danforth & Chu 2001).

### 2.4. Wolf-Rayet Stars – Interstellar/Circumstellar Bubbles

In the last decade, sensitive surveys of WR ring nebulae have been made using CCD cameras with interference filters for the Galaxy (Miller & Chu 1993; Marston et al. 1994a,b; Marston 1997) and the Magellanic Clouds (Dopita et al. 1994). The identification criteria in these surveys are similar to those used in the 70’s and 80’s. Based on nebular dynamics, Chu (1981) find three types of WR ring nebulae: R-type nebulae are not dynamically shaped by the WR stars, W-type nebulae are bubbles blown by the WR stars, and E-type nebulae are ejecta nebulae.

García-Segura et al. (1996a,b) show that a WR star is surrounded by an inner circumstellar bubble and an outer interstellar bubble. The interstellar bubble was blown by the progenitor of the WR star at the MS phase, while the circumstellar bubble is blown by the WR star in the circumstellar nebula ejected by the progenitor of the WR star at the RSG or LBV phase. Both the E-type and W-type WR ring nebulae refer to the circumstellar bubbles. If the circumstellar bubble is too tenuous to be detected, then the outer interstellar bubble may be ionized and become a R-type WR ring nebula.

It is thus clear that the three types of WR ring nebulae defined by Chu (1981) are not mutually exclusive. In particular, there is no physical distinction between W-type and E-type nebulae. Bearing in mind that the optically identified ring nebulae around WR stars may refer to either the circumstellar bubbles or the interstellar bubbles, it is thus meaningless to make statistical statements about “the percentage of WR stars in visible ring nebulae,” unless distinction is made between the circumstellar and interstellar cases.

WR ring nebulae have been studied in multiple wavelengths. Using the differences in H$`\alpha `$ and \[O iii\] images, Gruendl et al. (2000) presented a morphological diagnostic for dynamical evolution of WR bubbles. The most dramatic contrast between H$`\alpha `$ and \[O iii\] morphology is illustrated in the HST WFPC2 images of NGC 6888 in Fig. 4 (Moore, Hester, & Scowen 2000). H i 21-cm line observations have revealed neutral gas shells of diameters a few 10’s pc around WR stars; the properties of these shells are consistent with those expected in the interstellar bubbles blown by the MS progenitors (see review by Coppa 2002).

Observations of molecular gas associated with WR ring nebulae have yielded many intriguing results: (1) circumstellar CO associated with WR 16 (Marston et al. 1999) and possibly NGC 6888 (Rizzo et al. 2002); (2) complex molecules and shock-excited interstellar H2 and NH3 in NGC 2359 (St-Louis et al. 1998; Rizzo, Martín-Pintado, & Henkel 2001); (3) interstellar molecular gas around NGC 3199 (Marston 2001, 2002).

### 2.5. Blue Supergiants and Massive Binaries – Circumstellar Nebulae

The BSG phase has no unique meaning. Any massive star that is blue and luminous is called a BSG. If it varies, then it is called a LBV. The existence of a close binary companion can make the nature of a BSG even more mysterious. A small number of circumstellar nebulae around BSGs have been identified. Ironically, almost all of them have something peculiar about either the spectral properties of the stars or the physical structure of the nebulae. Below are a few examples.

NGC 6164-5 around HD 148937 consists of a small (2.6 $`\times `$ 1.9 pc) S-shaped circumstellar nebula and a filamentary interstellar ring nebula ($``$ 10 pc in diameter) inside a large dusty cavity structure (Leitherer & Chavarría-K. 1987). Clearly, HD 148937 is an evolved star, but its spectral type, O6.5f?p, is peculiar and it shows well-marked emission in C iii and N iii lines that are not seen in normal Of stars (Walborn 1972). What is the evolutionary status of HD 148937?

Sk$``$69 202 (B3 I), the progenitor of SN 1987A, and Sher 25 (B1.5 I) are both BSGs and have circumstellar ring nebulae of similar size and morphology (Brandner et al. 1997a, b). The hourglass morphology of these two rings are uncommon among known ring nebulae around massive stars; in fact, no other ring nebulae have such a morphology. Do these ring nebulae provide circumstellar or circumstantial evidence for a similar evolutionary status between these two BSGs?

Podsiadlowski, Joss, & Hsu (1992) have proposed that Sk$``$69 202 is a product of a massive interacting binary; their model explains the BSG spectral type of the star and the formation of the equatorial ring in the circumstellar nebula. Interestingly, a compact circumstellar nebula has been observed around the massive eclipsing binary RY Scuti (Smith et al. 1999); another circumstellar nebula has been detected around the B\[e\] supergiant R4 and is suggested to be ejected during a binary merger process (Pasquali et al. 2000). These circumstellar nebulae clearly have different formation mechanisms from those around single massive stars. The circumstellar nebula of Sher 25, resembling that around Sk$``$69 202, most likely has a similar formation mechanism, which then requires Sher 25 to involve a massive interacting binary. Unless observations of Sher 25 can exclude this possibility, the ring nebula of Sher 25 provides circumstellar evidence that it is similar to Sk$``$69 202.

### 2.6. Observations vs. Expectations

Observations of ring nebulae around single massive stars with known evolutionary status generally agree with expectations:
(1) MS O stars blow interstellar bubbles of sizes $``$ a few $`\times `$10 pc,
(2) RSGs are surrounded by circumstellar nebulae of sizes a few $`\times `$103 AU,
(3) LBVs have small circumstellar bubbles with radius $``$1 pc,
(4) WR stars have larger circumstellar bubbles, a few pc in size, surrounded by
interstellar bubbles blown by their progenitors at MS stage.
The circumstellar nebulae, consisting of material ejected by the stars, are N-enriched and can be diagnosed by a high \[N ii\]/H$`\alpha `$ ratio. LBVs may evolve into WR stars, and the sizes of their circumstellar bubbles (Chu, Weis, & Garnett 1999) support this evolutionary sequence.

On the other hand, ring nebulae around massive stars with uncertain evolutionary status, such as BSGs, show complex structures that are not completely understood. This is an area that needs a lot of future work. Ring nebulae around close massive binaries also need to be studied.

## 3. X-ray Views of Ring Nebulae

Bubbles blown by massive stars ought to be filled by shocked fast stellar wind at X-ray-emitting temperatures (Weaver et al. 1977). The first detection of diffuse X-ray emission from the hot gas in a bubble interior was made by Einstein IPC observations of NGC 6888 (Bochkarev 1988). ROSAT PSPC observations of several WR ring nebulae detected only NGC 6888 and S 308 (Chu 1994). The X-ray emission from NGC 6888 is brighter toward the dense shell rim, and the X-ray spectrum indicates a plasma temperature of 1.6$`\times `$106 K (Wrigge, Wendker, & Wisotzki 1994). ROSAT PSPC observations of S 308 is severely affected by the occultation ring, making the extraction of X-ray spectrum difficult and unreliable (Wrigge 1999). ASCA observations of NGC 6888 show an additional faint component at 8$`\times `$106 K (Wrigge et al. 1998).

The advent of Chandra and XMM-Newton finally makes it possible to study the distribution and physical conditions of the hot gas in bubble interiors. To date, NGC 6888 has been awarded a 100 ks Chandra ACIS-S observation in Cycle 4, and S 308 has recently been observed with XMM (Chu et al. 2002, in preparation). Figure 6 shows the X-ray image of the northwest quadrant of S 308. For the first time, the relative locations of the interior hot gas and the cool shell are resolved! A gap of 0.5 pc between the boundary of hot gas and the leading edge of the cool, dense bubble shell is observed; this gap would correspond to the interface region where heat conduction takes place.

It is interesting to compare the hot interior of S 308 to those of planetary nebulae (PNe). The progenitors of PNe have initial masses up to $``$10 M, which corresponds to mid-B type stars, which are qualified as “massive stars.” PNe are believed to be formed by the current fast stellar winds interacting with the former slow winds at the AGB phase (Kwok 1983), a mechanism almost identical to that for the formation of a WR bubble with a RSG progenitor (García-Segura et al. 1996). Diffuse X-ray emission has been detected in four PNe (Chu, Guerrero, & Gruendl 2002, and references therein). Figure 7 shows that the hot gas in PN interiors has temperatures of 2–3$`\times `$106 K, while S 308 has the coolest interior, only $``$1$`\times `$106 K. These temperatures are much lower than the post-shock temperature for a 1000–3000 km s-1 fast wind. Clearly, mixing with the cool nebular material has taken place (e.g., Pittard et al. 2001a,b). More X-ray observations of ring nebulae around massive stars are needed to constrain bubble models.

## 4. Ring Nebulae after the Supernova Explosion

One of the most intriguing questions about massive stellar evolution is “At what locations in the HR diagram do massive stars explode as SNe?” Only two SNe’s progenitors have adequate data to estimate their spectral type: B3 I for SN 1987A (Walborn et al. 1989) and K0 I for SN 1993J (Aldering, Humphreys, & Richmond 1994). Mass limits of SN progenitors have been estimated from photometric measurements using pre-explosion images (Smartt et al. 2001, 2002; Smartt 2002).

An alternative way to estimate a SN progenitor’s spectral type is through the physical properties of the circumstellar nebulae, which can be identified by narrow emission lines in high-dispersion spectra of SNe. The expansion of the nebula can be derived from the widths of line profiles, and the nebular density and temperature can be derived from the \[O iii\] and \[N ii\] lines. A variety of circumstellar nebulae have been observed.

SN 1997ab and SN 1997eg both showed a narrow P Cygni H$`\alpha `$ line, implying the existence of circumstellar nebulae with densities $`>`$107 cm-3 and expansion velocities 100–150 km s-1 (Salamanca et al. 1998; Salamanca, Terlevich, & Tenorio-Tagle 2002). These narrow P Cygni lines disappeared in observations made in March 1999, indicating that they are small (Gruendl et al. 2002). These circumstellar nebulae cannot be the RSG wind because the expansion velocities are too high. It is possible that these dense circumstellar nebulae are produced by a very dense SN progenitor wind, as the “deathbed ejecta” (Chu 2001).

SN 1978K (Fig. 8) shows narrow H$`\alpha `$ and \[N ii\] lines superposed on a broad H$`\alpha `$ line from the SN ejecta (Chu et al. 1999; Gruendl et al. 2002). A lower limit on the nebular size can be derived from the age of the SN and the expansion velocity of the SN ejecta, while an upper limit of 2.2 pc is derived from the unresolved HST WFPC2 image. The observed auroral to nebular line ratios of the \[O iii\] and \[N ii\] lines indicate a high density, 3-12$`\times `$105 cm-3 (Ryder et al. 1993). Such size and density are consistent with those of LBV nebulae. The progenitor of SN 1978K might have exploded at the LBV phase.

SN 1998S showed rapid variation in the narrow nebular line profiles during the first month after the SN explosion, indicating that the SN ejecta was plowed through a circumstellar nebula with a range of densities (Fassia et al. 2001). A year after the SN explosion, a narrow \[O iii\] line persisted; the FWHM of this line indicates an expansion velocity $`<`$25 km s-1 (Gruendl et al. 2002). The lack of \[N ii\] counterpart suggests an interstellar origin. Thus, this narrow \[O iii\] line might originate from an interstellar bubble.

Monitoring observations of SN 1987A have proven useful in revealing the interaction between the SN ejecta and the circumstellar nebula. Likewise, monitoring observations of extragalactic SNe with narrow nebular lines would also be very useful in determining the size of the circumstellar nebula. By comparing the density and size of a SN’s circumstellar nebula to those of ring nebular around massive stars, it is then possible to gain insight into the nature of the SN progenitor.

## References

Aldering, G., Humphreys, R. M., Richmond, M. 1994, AJ 107, 662

Bochkarev, N. G. 1988, Nature 332, 518

Brandner, W., Chu, Y.-H., Eisenhauer, F., et al. 1997a, ApJ 489, L153

Brandner, W., Grebel, E. K., Chu, Y.-H., Weis, K. 1997b, ApJ 475, L45

Cappa, C. E. 2002, in this volume

Chu, Y.-H. 1981, ApJ 249, 195

Chu, Y.-H. 1994, in AIP Conf. Proceedings “The Soft X-ray Cosmos,” 154

Chu, Y.-H. 2001, in AIP Conf. Proceedings “Young Supernova Remnants,” 409

Chu, Y.-H., Caulet, A., Montes, M. J., et al. 1999, ApJ 512, L54

Chu, Y.-H., Guerrero, M. A., Gruendl, R. A. 2002, in IAU Symposium 209 “Planetary Nebulae,” in press (astro-ph/0202509)

Chu, Y.-H., Weis, K., Garnett, D. R. 1999, AJ 117, 1433

Danforth, C. W., Chu, Y.-H. 2001, ApJ 552, L155

Dopita, M. A., Bell, J. F., Chu, Y.-H., Lozinskaya, T. A. 1994, ApJS 93, 455

Fassia, A., et al. 2001, MNRAS 325, 907

García-Segura, G., Langer, N., Mac Low, M.-M. 1996, A&A 316, 133

García-Segura, G., Mac Low, M.-M., Langer, N. 1996, A&A 305, 229

Gervais, S., St-Louis, N. 1999, AJ 118, 2394

Gruendl, R. A., Chu, Y.-H., Dunne, B. A., Points, S. D. 2000, AJ 120, 2670

Gruendl, R. A., Chu, Y.-H., Van Dyk, S. D., Stockdale, C. J. 2002, AJ 123, 2847

Humphreys, R. M. 2002, in this volume

Johnson, H. M., Hogg, D. E. 1965, ApJ 142, 1033

Kwok, S. 1983, in IAU Symposium 103 “Planetary Nebulae,” 293

Lamers, H. J. G. L. M., Nota, A., Panagia, N., Set al. 2001, ApJ 551, 764

Leitherer, C., Chavarría-K., C. 1987, A&A 175, 208

Marston, A. P. 1997, ApJ 475, 188

Marston, A. P. 2001, ApJ 563, 875

Marston, A. P. 2002, in this volume

Marston, A. P., Chu, Y.-H., García-Segura, G. 1994a, ApJS 93, 229

Marston, A. P., Welzmiller, J., Bransford, M. A., et al. 1999, ApJ 518, 769

Marston, A. P., Yocum, D. R., García-Segura, G., Chu, Y.-H. 1994b, ApJS 95, 151

Miller, G. J., Chu, Y.-H. 1993, ApJS 85, 137

Moore, B. D., Hester, J. J., Scowen, P. A. 2000, AJ 119, 2991

Nazé, Y., Chu, Y.-H., Points, S. D., et al. 2001, AJ 122, 921

Nota, A., Livio, M., Clampin, M., Schulte-Ladbeck, R. 1995, ApJ 448, 788

Nota, A. et al. 1996, ApJS 102, 383

Nota, A., Smith, L., Pasquali, A., Clampin, M., Stroud, M. 1997. ApJ 486, 338

Pasquali, A., Nota, A., Clampin, M. 1999, A&A 343, 536

Pasquali, A., Nota, A., Langer, N., et al. 2000, AJ 119, 1352

Pasquali, A., Schmutz, W., Nota, A., Origlia, L. 1997, A&A 327, 265

Pittard, J. M., Dyson, J. E., Hartquist, T. W. 2001a, A&A 367, 1000

Pittard, J. M., Hartquist, T. W., Dyson, J. E. 2001b, A&A 373, 1043

Podsiadlowski, P., Joss, P. C., Hsu, J. J. L. 1992, ApJ 391, 246

Rizzo, J. R., Martín-Pintado, J., Desmurs, J.-F. 2002, in this volume

Rizzo, J. R., Martín-Pintado, J., Henkel, C. 2001 ApJ 553, 181

Ryder, S., Staveley-Smith, L., Dopita, M. A., et al. 1993, ApJ 416, 167

Salamanca, I., et al. 1998, MNRAS 300, L17

Salamanca, I., Terlevich, R., Tenorio-Tagle, G. 2002, MNRAS 330, 844

Smartt, S. J. 2002, in this volume

Smartt, S. J., Gilmore, G. F., Tout, C. A., Hodgkin, S. T. 2002, ApJ 565, 1089

Smartt, S. J., Gilmore, G. F., Trentham, N., et al. 2001, ApJ 556, L29

Smith, L. J., Nota, A., Pasquali, A., Leitherer, C., et al. 1998, ApJ 503, 278

Smith, N., Humphreys, R. M., Davidson, K. et al. 2001, AJ 121, 1111

Smith, N., Gehrz, R. D., Humphreys, R. M., et al. 1999, AJ 118, 960

St-Louis, N., Doyon, R., Chagnon, F., Nadeau, D. 1998, AJ 115, 2475

Walborn, N. R. 1972, AJ 77, 312

Walborn, N. R., et al. 1989, A&A 219, 229

Weaver, R., McCray, R., Castor, J., Shapiro, P., Moore, R. 1977, ApJ 218, 377

Weis, K. 2002, in this volume

Wrigge, M. 1999, A&A 343, 599

Wrigge, M., Chu, Y.-H., Magnier, E. A., Kamata, Y. 1998, in Lecture Notes in Physics, vol.506, 425

Wrigge, M., Wendker, H. J., Wisotzki, L. 1994, A&A 286, 219


text:
# The phase-space structure of cold dark-matter halos: Insights into the Galactic halo

## 1 Introduction

Over the last twenty years, the hierarchical paradigm has emerged as the standard model to describe the formation of structure in the Universe. As embodied in the current “concordance” $`\mathrm{\Lambda }`$CDM model it appears to be consistent with a very wide range of cosmological data ranging from fluctuations in the Cosmic Microwave Background through the structure of Ly$`\alpha `$ forest absorption in QSO spectra and the gravitational shear induced by dark-matter structures to the observed large scale structure in the galaxy distribution. An important characteristic of such models is that they are based on a set of well-defined and testable assumptions. This renders possible the detailed modelling of the formation and evolution of galactic systems, and a later comparison to observations of the properties of these systems as a function of environment or redshift (e.g. Diaferio et al. 2001; Benson et al. 2001; Somerville, Primack & Faber 2001).

It is also possible to test the hierarchical paradigm on our Galaxy (e.g. Hernández, Avila-Reese & Firmani 2001). Several groups (Moore et al. 1999; Klypin et al. 1999; Klypin, Zhao & Somerville 2002) have focused on the properties of dark halos, hoping to constrain the nature of dark-matter. These groups performed high resolution simulations of a galactic size halo in CDM cosmologies. They confirmed earlier analytic claims (Kauffmann, White & Guiderdoni 1993) that the predicted number of satellites exceeds the known population in the Local Group by a factor of ten. Some attempts have been made to account for the disagreement, by changing the nature of the dark-matter (Spergel & Steinhardt 2000; Bode, Ostriker & Turok 2001), by modifying the initial power spectrum of density fluctuations (Kamionkowski & Liddle 2000) or by taking into account the effects of a reionising background which may inhibit star formation in the smallest mass halos (Kauffmann et al. 1993; Bullock, Kravtsov & Weinberg 2001; Benson et al. 2002). The recent results by Kleyna et al. (2002) on the mass distribution in the Draco dSph (see also Mateo 1997 and Lokas 2001 for a similar study on Fornax) favour an astrophysical explanation since the actual circular velocities of the other satellite galaxies of the Milky Way are in fact several tens of km s-1 larger than previously thought, and agree with those expected for the most massive substructures in a $`\mathrm{\Lambda }`$CDM universe (Stoehr et al. 2002).

Broadly speaking, the hierarchical paradigm predicts that the Milky Way formed through mergers of smaller systems (White & Rees 1978). These systems would contribute to the dark halo, the spheroid (the bulge and the stellar halo) and to the Galactic gas reservoir. It may be very difficult to determine the relative gas contribution of these progenitor objects to the present Galaxy, since gas “easily forgets” its site of origin. However for collisionless stars and dark-matter the situation can be quite different. If the dynamical mixing timescales are sufficiently long (i.e. longer than the age of the Universe) it may be possible to “break-up” the Galactic spheroid (stars and may be even dark matter particles) into coherent structures in phase-space directly related to the systems that merged to form the Milky Way we observe today.

A first attempt at determining whether the merging history of the Milky Way may be imprinted in the phase-space structure of nearby halo stars, and thus be recovered, was made by Helmi & White (1999; hereafter HW). They studied the infall of satellites onto a fixed Galactic potential, and the evolution of the debris in phase-space. They found that after 10 Gyr stars having a common origin are distributed smoothly in space, but appear very clumped in velocity space, where they define streams with very small velocity dispersions. The expected number of such streams scales with the initial size $`r`$, velocity dispersion $`\sigma `$ and orbital period $`P`$, of the disrupted object:

$$N_{\mathrm{stream}}10\left(\frac{r}{1\mathrm{kpc}}\right)^2\left(\frac{\sigma }{15\text{km s}\text{-1}}\right)\left(\frac{P}{0.23\mathrm{Gyr}}\right)^3.$$
(1)
The total number of stars associated with the object is $`N_{}r\sigma ^2`$ (from the virial theorem) while the volume $`V`$ over which they are spread scales with the cube of the size of the orbit, and so approximately as $`P^3`$. Hence the number of stars per stream in the Solar neighbourhood scales as $`N_{}/VN_{\mathrm{stream}}\sigma /r`$; objects with large initial velocity dispersion and small initial size should produce the most easily detectable streams with little dependence on initial period. Such arguments suggest that the Solar neighbourhood velocity ellipsoid is composed of $`300500`$ kinematically coherent structures which originated in past merger and accretion events. A pair of halo streams that can perhaps be directly linked to a disrupted satellite were detected in the Solar neighbourhood by Helmi et al.(1999). The progenitor of these two streams was probably similar to the dwarf galaxy Fornax. Substructure in the outer halo also appears to be ubiquitous, and has been found by several surveys over the last few years (e.g. Ivezic et al. (2000) and Yanny et al. (2000) for the SDSS; Dohm-Palmer et al. (2001) for the SPS; Vivas et al. (2001) for QUEST). Most of these recently discovered structures can be associated to just one of Milky Way’s satellites: the Sagittarius dwarf galaxy which is in the process of being completely disrupted (Ibata et al. 2001; Martínez-Delgado et al. 2001; Helmi & White 2001).

A weak point of the HW analysis and of similar studies (e.g. Johnston, Hernquist & Bolte 1996; Johnston 1998), is the assumption of a fixed, smooth potential onto which galaxies are accreted. In hierarchical clustering, galaxy potentials are constantly changing, and can vary very violently during mergers. Large numbers of clumps orbit the centre of even a “virialised” halo. These clumps may have substantial effects on the structure of debris streams (e.g Johnston, Spergel & Haydn 2002; Ibata et al. 2002; Mayer et al. 2002).

The main goal of the present paper is to understand the phase-space structure of cold dark-matter halos. In particular, we want to study the evolution of satellite debris, and to quantify the expected amount of substructure. We also want to determine to what extent previous results are valid in the truly hierarchical regime of the build–up of a galaxy. We tackle these problems by scaling down to galactic size a high-resolution simulation of the formation of a cluster in a $`\mathrm{\Lambda }`$CDM cosmology (Springel et al. 2001).

The paper is organised as follows. In Sec.2 we describe the simulations, in Sec.3 we follow the evolution in phase-space of debris streams, and compare to the analytic model of HW in Section 3.2. Sec.4 describes the mass-growth history of the simulated dark-matter halo, and in Section 5 we determine the number of streams and their internal properties as function of distance from the dark-matter halo centre. We leave the summary and discussion of our results for Section 6.

## 2 Methodology

### 2.1 The simulations

The simulation we analyse here was carried out using a parallel tree-code (Springel, Yoshida & White 2001) on the Cray T3E at the Garching Computing Centre of the Max Planck Society. Its initial conditions were generated by zooming in and re-creating with higher resolution a particular galaxy cluster and its surroundings formed within a cosmological simulation (as in Tormen, Bouchet & White 1997). The original parent $`\mathrm{\Lambda }`$CDM cosmological simulation (from Kauffmann et al. 1999) had parameters $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.7`$ and $`\sigma _8=0.9`$. The cluster selected for re-simulation was the second most massive cluster in this simulation, having a virial mass of $`8.4\times 10^{14}h^1\mathrm{M}_{}`$. The particles that end up in the final cluster of the cosmological simulation and in its immediate surroundings (defined by a comoving sphere of $`70h^1`$ Mpc radius) were traced back to their Lagrangian region in the initial conditions for re-simulation. The initial mass distribution between $`21`$ and $`70h^1`$ Mpc was represented by $`3\times 10^6`$ particles. In the inner region, where the original simulation had $`2.2\times 10^5`$ particles, new initial conditions were created using $`6.6\times 10^7`$ particles. Small scale power was added in accordance with the better k-space sampling allowed by the larger number of particles. The original force softening was also decreased to obtain better spatial resolution. The new simulation was run from high redshift until $`z=0`$, and was analysed in considerable detail in Springel et al. (2001). In this high resolution simulation there are about 20 million particles within the virialised region of the cluster halo.

### 2.2 Scaling to a Milky Way halo

In Figure 1 we show the circular velocity profile of the cluster at $`z=0`$. We determine the cluster centre by successively refining a mesh located on the cluster, and determining the cell containing the largest number of particles. This process is repeated until the largest number of particles in a given (now small size) cell is sufficiently small to determine by simple counts which particle has the largest number of neighbours. Its position corresponds to the point of maximum density and it is defined as the cluster centre. This determination is robust against changes in the mesh shape and size, and we estimate the error in the final position of the centre of the cluster to be of the order of 0.7 kpc, comparable to the scale of the gravitational softening used in the simulation. The circular velocity is then determined after spherically averaging the mass distribution around the centre of the cluster and is derived from $`V_c(r)=\sqrt{GM(r)/r}`$.

By fitting a Navarro, Frenk & White (NFW; 1996) profile, we find that the concentration of the cluster halo is $`c_{\mathrm{NFW}}=7.3`$, the scale radius $`r_s=285`$ kpc and the virial radius $`r_{200}=2090`$ kpc.

We can scale the cluster to a “Milky Way” halo, by requiring that its maximum circular velocity be equal to 220 km s-1. The required scaling factor $`\gamma `$ is given by

$$\gamma =\frac{V_c^{\mathrm{cl}}}{V_c^{\mathrm{MW}}}=\frac{r_{200}^{\mathrm{cl}}}{r_{200}^{\mathrm{MW}}}9.18.$$
(2)
With this scaling, the virial radius of our simulated Milky Way dark-matter halo is 228 kpc. Its virial mass is $`1.28\times 10^{12}\mathrm{M}_{}`$, the mass of an individual particle is $`8.66\times 10^4\mathrm{M}_{}`$, and the gravitational force in the final object has an equivalent Plummer softening of 0.11 kpc.

Our argument that with the simple scaling of Eq.(2) this simulation is a fair representation of the growth of the Milky Way relies on both theoretical and numerical results (Lacey & Cole 1993; Moore et al. 1999). Numerical simulations by Moore et al. (1999) have shown that galaxy and cluster halos have similar properties in terms of final structure (density profile), number of satellites, etc., despite typically assembling at systematically different redshifts. Jing & Suto (2000) have also performed high resolution simulations of several galaxy, group and cluster size halos. Although they find that galaxy halos have a steeper inner profile than larger mass halos, they also find that the scatter in the properties of objects belonging to the same class is as large as the systematic differences between classes. We expect therefore, that our scaled simulation will represent reasonably well the formation process of a galactic halo, except that its assembly occurs at later times than the majority of such halos. From now on we will use this scaling in the paper, and assume, unless otherwise stated, that the simulation represents a galaxy halo.

## 3 The phase-space evolution of debris

Here we study the phase-space evolution of debris from the disrupted halos that end up forming the “Galactic” dark halo at the present time. With this goal in mind, we proceed by identifying halos at high-redshift which are directly accreted onto the main progenitor of the “Milky Way’s” halo.

We identify halos at each output time using a Friends-of-Friends (FOF) algorithm, which links particles separated by less than 20% of the mean interparticle separation. In this way we can construct a detailed merger history of the galaxy. As we step back in redshift we identify at each output the most massive halo which is part of the galaxy’s main progenitor at the subsequent time. We say that a halo identified at redshift $`z`$ will be directly accreted onto the main progenitor at $`z^{}`$ (the redshift of the next simulation output) if at least half of its particles and the most-bound particle have become part of the main progenitor at $`z^{}`$. Here the most-bound particle refers to the particle with the minimum potential energy in the (to be) accreted halo.

The centre of mass position $`𝐩_{CM}`$ and velocity $`𝐯_{CM}`$ of the main progenitor are computed as follows. We first determine the location of the most bound particle $`𝐩_{\mathrm{mb}}`$. We define concentric spheres of successively smaller radii (down to 10 kpc radius) around this particle, and compute the centre of mass $`𝐩_{}^{}{}_{CM}{}^{}`$ from the particles within these spheres. We stop this iterative procedure when $`|𝐩_{\mathrm{mb}}𝐩_{}^{}{}_{CM}{}^{}|<ϵ`$, where $`ϵ=0.3`$ kpc. This is then defined as the location of centre of mass of the main progenitor. The velocity of the centre of mass is then $`𝐯_{}^{}{}_{CM}{}^{}`$, measured by the velocities of the particles within the largest sphere for which the above condition is satisfied.

This is like following the (thickest) trunk of the merger tree, which would correspond to the galaxy’s main progenitor, and studying what happens to halos which join from other tree branches as time progresses. The idea is illustrated in Figure 2.

For our high resolution simulation we have identified 752 halos with at least 1000 particles (which corresponds to a minimum mass of $`8.66\times 10^7\mathrm{M}_{}`$), which fall onto the main progenitor between redshift $`z=2.4`$ and the present day. The accreted halos have a large spread in mass as shown in the bottom panel of Figure 3. From the top panel of this Figure we note that in some cases the satellite-to-primary mass ratio is close to unity, corresponding to a major merger. Such mergers, although few in number, contribute a substantial fraction of the total mass growth.

### 3.1 Spatial and kinematic evolution of halo debris

We now study in detail how the debris streams from one of the above 752 halos evolve in time. We follow a halo that merges with the main progenitor at redshift 1.8, and has an initial mass of $`4.3\times 10^{10}\mathrm{M}_{}`$. This halo is denoted with a diamond in Figure 3. Its mass at the time of accretion was about 25% that of the main progenitor. This satellite decays through dynamical friction to the centre of the main progenitor where it is fully disrupted. The mass stripped off at earlier times mainly populates the outer regions of the final object, while that lost at late times ends up closer to its centre.

#### 3.1.1 Evolution of debris in the outer galaxy

We identify material from the accreted halo which is part of a tidal stream in the outer galaxy at the present time. We select a reference particle in this structure, which can be traced back to $`z=1.8`$, and then followed forwards in time. In Figure 4 we show a time sequence of the spatial distribution of particles lost at $`z=1.8`$. We say that a particle has been lost by its progenitor halo if its binding energy has become positive. The binding energy $`ϵ`$ of a particle located at $`𝐫`$ and with velocity $`𝐯`$ with respect to the centre of mass of the satellite111The centre of mass is here defined by the 0.1% most bound particles of the system, which in this case corresponds to 500 particles. is defined as $`ϵ=\mathrm{\Phi }(𝐫)+0.5|𝐯|^2`$, where $`\mathrm{\Phi }(𝐫)`$ is the potential energy at $`𝐫`$ due to all particles in the satellite. From Fig. 4 we see that the initial distribution of particles is relatively compact, and that as time passes by the material is strung out in a characteristic stream-like structure (e.g at $`t=4.31`$ Gyr) over several tens to hundreds of kiloparsecs. At the final time the material appears to be more smoothly distributed over the whole box, which is 400 kpc (in physical units) on a side. Our simulations are not well-suited to address the effects of other dark-matter lumps on tidal tails (cf Johnston et al. 2001), which could be related to some of the transient structures observed in this Figure. Our main limitation is the large time interval between stored outputs which prevents us from determining the effect of close encounters and their relation to the features seen in Figure 4.

In many of the snapshots tight small substructures can be observed. These objects were already present as subhalos within the satellite halo before it was accreted by the galaxy. After accretion, they were released from their parent satellite, becoming subhalos of the main galaxy. It is worth mentioning that these subhalos only constitute a small fraction of the debris material lost by the satellite (less than 10%, e.g. Ghigna et al. 2000; Springel et al. 2001).

In Figure 5 we show the velocities of debris particles that are relatively close (inside a cube of side 20 kpc) to the reference particle, at three different times. For comparison, note that the virial radius of the satellite at the time of infall was approximately 37 kpc. The solid grey circles correspond to particles which

1. were always neighbours of the reference particle: to be in this set, particles need to have $`x`$, $`y`$ and $`z`$ coordinates within 10 kpc of those of the reference particle in all previous outputs;
2. are within 5 kpc of the reference particle in each coordinate in the current output.

We note that the initial velocity distribution of the halo is broad, and relatively clumpy. As discussed before, this clumpiness reflects the internal structure of the object that fell in. As time goes by, the motions of the neighbours of the reference particle become more similar, and the velocity-box is empty except for velocities close to that of the reference particle (middle row). At later times (bottom row), other moving groups are visible, showing that the system has now produced multiple intersecting streams even in the outer galaxy.

The evolution of the velocity dispersion $`\sigma `$ in the neighbourhood of the reference particle is shown in the top panel of Figure 6. Here $`\sigma `$ is defined as $`\sqrt{\sigma _x^2+\sigma _y^2+\sigma _z^2}/\sqrt{3}`$, and is therefore independent of the choice of the coordinate system. It is measured for the set of particles shown as grey solid circles in Figure 5, that is for particles that satisfy both conditions (i) and (ii) for being neighbours of the reference particle. We also measure $`\sigma `$ for a subset of even closer neighbouring particles (whose $`x`$, $`y`$, and $`z`$ coordinates in all outputs before the one under study were within 5 kpc of those of the reference particle, and in the current output are within 2.5 kpc). The different values of the velocity dispersion in the stream obtained for these two cases are due to velocity gradients, which can be as large as the measured dispersions themselves. Note as well the decrease by roughly a factor 10 in the velocity dispersions in only 2 Gyr.

The stream’s density is shown in the bottom panel of the same figure. It is measured at each output by the number of neighbours in the stream in boxes of a given size (either of 5 or 10 kpc on a side) around the reference particle. For the larger volume, this number evolves from being slightly larger than five thousand to only thirteen after less than 2 Gyr of evolution. This implies that after this time, the number counts are dominated by Poisson noise.

#### 3.1.2 Evolution of debris closer to the centre of the galaxy

To follow the evolution of streams closer to the centre of the galaxy we focus on material lost a few Gyr after the satellite halo shown as a diamond in Fig.3 was accreted. We here focus on debris lost from this satellite halo at $`t=7.16`$ Gyr, or 3.5 Gyr after infall. This is deposited at an intermediate distance from the centre of the galaxy. We decided against tracking material lost even later because it turns out to be very difficult to follow streams from such material with our numerical resolution and number of stored outputs. The material lost at later times mixes on shorter timescales. Streams more rapidly reach very small densities, and we are unable to determine their properties reliably.

Figure 7 shows the evolution of the spatial distribution of the particles lost at $`t=7.16`$ Gyr. Note that the particles are more smoothly distributed than in the case of the outer halo material shown in Figure 4. This is because the inner regions of (satellite) halos are in general much less lumpy than the outskirts due to the shorter dynamical timescales.

As in the previous section we choose a reference particle (amongst all those lost at $`t=7.16`$ Gyr) that orbits through the intermediate regions of the galaxy. Proceeding as before, we follow the evolution of the velocities of particles located close to the reference particle at all times. Figure 8 shows the velocity distribution in the neighbourhood of the reference particle. By comparing to the analogue for the outer halo stream (Fig. 5) we note that, not just the spatial distribution, but also the velocity distribution is much smoother initially. This is also true at late times, and very rapidly a regime is reached where multiple streams can be observed. Figure 9 corresponds to the evolution of the velocity dispersion and density of a stream. Again we are not able to follow the evolution of the velocity dispersion of the stream at late times, because we run out of neighbouring particles quite rapidly, particularly for the smaller box. After $`t=10.5`$ Gyr, or 2 Gyr of evolution, the number of neighbours is dominated by Poisson noise. After this time, the values of the 1-d velocity dispersion, when measured, are based on the velocities of only two or three particles.

The evolution of the intermediate and outer halo streams is characteristic of all streams originating in directly accreted halos. Naturally the properties of streams and their location in the galaxy halo will depend on their progenitor, in particular on its initial mass. In this respect, streams originating in smaller halos are narrower, more clearly defined, and, typically, they phase-mix on longer timescales.

### 3.2 Mixing in phase-space: Comparison to analytic models of stream evolution

In the previous section we found a very rapid decline in the density and velocity dispersion of streams, especially in the first few Gyr of evolution. After this initial period, the number of particles in a stream is so small in our simulation that it is dominated by Poisson noise. In this regime, it is difficult to quantify the properties of the stream, and it is therefore hard to determine whether the lack of variation of the velocity dispersion at late times is due to numerical limitations or is a real effect.

To gain insight into these issues, we will analyse the expected behaviour of streams evolving in a smooth time-independent potential, which should resemble that of the galaxy halo. Using the approach developed by HW, we can follow the evolution of streams produced in a spherical and static NFW potential. The basic idea here is to map the initial system onto action-angle space, then follow the much simpler evolution in this space, and finally transform back locally onto observable coordinates (all these being linear transformations; for details see HW). This method, which uses action-angle variables, is limited to applications in which the potential is separable (cf Goldstein 1953; Binney & Tremaine 1987). This includes all spherically symmetric potentials but only few axisymmetric and triaxial ones, such as the general class of Stäckel potentials (e.g. Lynden-Bell 1962; de Zeeuw 1985; Dejonghe & de Zeeuw 1988).

We also approximate the phase-space density around the reference particle by a multivariate Gaussian distribution. This is possible because the multivariate Gaussian is determined from the properties of particles in a volume much smaller than the size of the halo.

As discussed in Section 2.2, the (galaxy) halo can be fit by an NFW profile:

$$\rho (r)=\rho _0\frac{\delta _cr_s^3}{r(r+r_s)^2},$$
(3)
where $`\rho _0=3H^2(z)/(8\pi G)`$ and $`\delta _c`$ is a function of the concentration of the halo $`c_{\mathrm{NFW}}=r_{200}/r_s`$:

$$\delta _c=\frac{200}{3}\frac{c_{\mathrm{NFW}}^3}{\mathrm{log}(1+c_{\mathrm{NFW}})c_{\mathrm{NFW}}/(1+c_{\mathrm{NFW}})}.$$
(4)
Recall that $`H(z)=H_0\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_0(1+z)^3}`$. The potential associated with this density can be obtained by integrating Poisson’s equation, and is found to be:

$$\mathrm{\Phi }(r)=\mathrm{\Phi }_0(z)\frac{r_s}{r}\mathrm{log}\left(1+\frac{r}{r_s}\right).$$
(5)
Here $`\mathrm{\Phi }_0(z)=3/2H^2(z)\delta _cr_s`$. The values of the parameters $`r_s`$ and $`r_{200}`$ given in Section 2.2 should be multiplied by the factor $`\gamma `$ for the scaled halo.

In this (spherically averaged) potential we integrate the orbit of the reference particle in time. This integration is done in two complementary ways:

1. starting from the position and velocity of the reference particle at the “time of formation of the stream” (this is the time when the particles become unbound from their parent halo). In this case, the integration is performed forwards in time, and the potential used is given by Eq. (5), where $`z`$ is the redshift of formation of the stream.
2. starting from the position and velocity of the reference particle at the present time. The orbit is then integrated backwards in time, until the “time of formation of the stream”, in the present-day potential ($`z=0`$).

In both cases, the orbits are integrated in the reference frame of the centre of mass of the main progenitor of the galaxy.

In the top panel of Figure 10 we plot the radial oscillations of the orbit for the outer halo reference particle. The solid curve corresponds to the NFW potential at redshift $`z_{\mathrm{form}}=1.8`$, while the dashed curve to that at $`z=0`$. We note that for the reference particle in the outer halo stream, the best agreement is obtained when the integration is performed backwards in time. We also note that the potential seems to have fluctuated dramatically until $`t6`$ Gyr, inducing strong changes in the radial oscillations of the reference particle. As a whole, the orbit has evolved to a more bound state, due to the aggregation of mass during mergers.

The analogous plots for the reference particle orbiting the intermediate halo are shown in the bottom panel. Recall that here $`z_{\mathrm{form}}=0.73`$. In this case, none of the proposed orbits fit the actual orbit very well; it may be considered to lie in between these two cases. Thus, the characteristics and the evolution of the debris streams predicted for the two proposed orbits may perhaps encompass the actual behaviour of the streams in our simulations. In this “intermediate halo” case, we do not find a clear indication of evolution in the orbit. This is likely due to the mild increase in the mass within the orbit since $`z_{\mathrm{form}}`$ (see Sec. 4)

In Figures 11 and 12 we show the predicted evolution of the 1-d velocity dispersion and the behaviour of the density in the spherical static NFW potential for the outer and intermediate halo streams, respectively. The actual behaviour of the streams in the simulation is also shown for direct comparison. In both cases, we find that the decrease observed in the simulated streams is close to what is predicted for the evolution of streams moving in a static potential. This is the case at least for the first few Gyrs of evolution. Although after this period we expect some differences, due to the fact that the orbit evolves in shape, the agreement still appears to be quite good. We predict a rapidly varying velocity dispersion (and density) on top of a secular evolution. The spikes take place when particles in a stream go through a caustic surface defined by their orbital turning points. To observe this behaviour in the simulation would require output times spaced by 1/2 of the radial period at most, or roughly 0.25 Gyr for the stream orbiting the intermediate regions of the halo. Our outputs, on the other hand, are logarithmically spaced in the cosmological expansion factor $`a(t)`$, so that at $`z=1`$ the time elapsed between two outputs is $`\mathrm{\Delta }t=0.52`$ Gyr, and by $`z=0`$, $`\mathrm{\Delta }t=0.86`$ Gyr.

## 4 Mass growth history of the galaxy

Here we focus on determining the mass growth history of the halo as a function of distance from the halo centre. This is relevant in two different ways. First, debris from satellites accreted at late times will be generally less mixed, and could thus produce more massive streams. Determining where to expect these streams will enable us to understand the properties of the dark halo of our Galaxy. Secondly such satellites will have had more time to form stars in them (prior to their merging), thereby providing the galaxy with younger stars. The time of merging could thus be used as an indicator of the expected age distribution of stars in different regions of our Galaxy.

We select all satellites that merged with the galaxy since redshift $`z=2.4`$, and determine when these mergers took place and what is their final debris distribution. We proceed by dividing the halo in six spherical shells around the galaxy centre. These shells are located at: $`r<10`$ kpc, $`10r<25`$ kpc, $`25r<50`$ kpc, $`50r<75`$ kpc, $`75r<100`$ kpc and $`100r<200`$ kpc. For each particle in a given shell, we determine when it was accreted by the main progenitor of the galaxy. In Figure 13 we show the fraction of mass accreted normalised to the present mass for each shell as a function of redshift (and time). We note that the formation time of the inner galaxy is strongly biased towards high redshifts, with more than 60% of the mass already present at $`z=2.4`$ or 11 Gyr ago. Conversely, we find that mergers in the last 3 Gyr, contributed a relatively small amount of mass, less than 0.1%, to this region of the galaxy halo. For the outer regions of the galaxy, we note that the growth is much more gradual in time, with accretion being almost equally important at all times.

Of course, the detailed shape of the mass growth history depends on the detailed merger history of the halo, since the peaks observed in the different histograms correspond to individual mergers taking place at those times. Nonetheless, several authors have found that there is an almost “universal” form of the mass accretion histories of galaxies and galaxy clusters (e.g. van den Bosch 2002; Wechsler et al. 2002), which would also represent reasonably well the growth of our simulated halo. On the other hand, Zhao et al. (2002) have shown using high-resolution simulations that the mass in the inner regions of dark-matter halos (which they define as the mass within a scale radius, which in our case corresponds to approximately 30 kpc), is generally in place by $`z2`$. These results suggest that the formation history of our simulated halo is rather typical.

Another way of understanding the age structure of the galaxy halo is by focusing on what fraction of the mass was in place by different redshifts as a function of distance from the galaxy centre. To some extent this is the cumulative distribution of the plots shown in Fig. 13 for each shell. We focus on the fraction of the present-day mass in the shell that was in place at four different redshifts $`z=2.4`$, 1.5, 0.84, and 0.35. The results are shown in Figure 14. Again we notice that 60% of the mass in the inner galaxy was in place by $`z=2.4`$, and more than 90% by $`z=1.5`$, i.e. 9.5 Gyr ago. Thus any accreted stellar populations in the inner regions of the galaxy are predicted to be old. On the other hand, only about half of the particles in the intermediate regions of the galaxy were present 9.5 Gyr ago (90% by $`z=0.84`$). The formation of the outer galaxy is more biased towards late times, with half of its particles coming into place in the last 7 Gyr.

As a cautionary remark, let us recall that although we discuss here the expected “stellar” populations, we are actually describing the growth in mass of the dark-matter halo, rather than of the stellar halo. To describe properly the age distribution of the stellar halo we would need to take into account that the stars in the different satellites that merged to build up the halo probably were located primarily in their innermost regions. This would mean, as shown in Section 3, that they should end up closer to the galactic centre than most of the dark-matter particles of the same satellite (see also White & Springel 2000). Moreover, we are also ignoring the fact that satellites orbiting the intermediate and outer regions of the halo can survive until the present day. In those cases it is possible that they continue to form stars, even while orbiting inside the dark halo of the galaxy (like the Magellanic Clouds). This means that the material that is stripped off from these systems could also contain younger stars.

## 5 The number of streams

The results obtained in Section 3 suggest that the number of streams in the inner galaxy, and in particular in the vicinity of the Sun should be quite large. In this section we will estimate this quantity, as well as other characteristics of the streams, such as internal velocity dispersion and density. We will compute the number of streams inside cubes located in the inner $`([60,60])^3`$ kpc3 in the scaled “Milky Way” halo. We partition this space into boxes of 2 kpc on a side.

### 5.1 General definitions

Let there be $`L`$ separate true streams in a given box, each having mass $`\nu _k`$ ($`k=1,L`$). Let us assume that the observed number of particles in a stream follows a Poissonian distribution, and that the actual count is $`N_k`$ particles ($`k=1,L`$). Then the expectation value of $`N_k`$ is $`N_k=\nu _k`$, and $`N_k^2=\nu _k+\nu _k^2`$. (Here we measure mass in units of the particle mass of the simulation.)

The mass-weighted mean mass per stream $`\mu `$ is

$$\mu =\frac{_{k=1,L}\nu _k^2}{_{k=1,L}\nu _k}.$$
(6)
Therefore our estimate for $`\mu `$ corrected for Poisson noise will be

$$\widehat{\mu }=\frac{_{k=1,L}\left(N_k^2N_k\right)}{_{k=1,L}N_k}.$$
(7)
In the limit of very massive streams, the Poisson correction will be negligible since $`N_k^2N_k`$. In the limit of small number of particles per stream, the correction will be of the same order as the quantity we measure. Note as well, that streams which in this realization do not have any particle or just have one, do not contribute to the numerator of Eq. (7). However single particle streams do contribute to the total number of particles in the box. If in a box we find mostly one-particle streams, then $`\widehat{\mu }`$ can (correctly) become very small.

We define the mass-weighted number of streams $`F`$ in a box as the ratio of the total mass in the box to the mass-weighted mean mass $`\mu `$ per stream in the box. Therefore $`F`$ is

$$F=\frac{_{k=1,L}\nu _k}{\mu },$$
(8)
and our Poisson corrected estimate of $`F`$ is

$$\widehat{F}=\frac{_{k=1,L}N_k}{\widehat{\mu }},$$

or

$$\widehat{F}=\frac{\left(_{k=1,L}N_k\right)^2}{_{k=1,L}\left(N_k^2N_k\right)}.$$
(9)
For example, if in our realization all of the streams have only two particles $`N_k=2`$, our estimate of the mass-weighted filling factor becomes $`\widehat{F}N_T`$, where $`N_T`$ is the total number of particles in the box. In the regime where one massive stream dominates the distribution, $`\widehat{F}N_T/(N_T1)`$ and will thus be close to unity. Also note that $`\widehat{F}`$ can be larger than $`N_T`$ which will happen when most of the streams contain only 1 particle.

### 5.2 The minimum number of streams: the number of halos

We would like to obtain an estimate of the number of disrupted halos contributing to the density at each location in the galaxy. This is a lower limit to the total number of streams present, since a halo can (and usually does) give rise to multiple and spatially overlapping structures, as shown in Sec. 3. This lower limit will be particularly unrepresentative of the true number of streams in the inner halo both because of the very short dynamical timescales there and because more than 60% of the mass in this region comes from just one object identified at $`z=2.4`$, as shown in Sec.4.

To obtain an estimate of the number of halos contributing to any given position in the Galaxy, we determine which of the halos identified at $`z=2.4`$ contribute to each given box and with how many particles. In Figure 16 we plot the (Poisson corrected) mass-weighted mean number of particles per halo $`\widehat{\mu }_{\mathrm{halo}}`$ as a function of distance from the galaxy centre for each of the boxes considered. The thick black line corresponds to the mean number of particles in a box, averaged over all boxes at the same location. This Figure shows that most of the particles in the inner galaxy come only from a handful of disrupted halos. This can also be seen from Figure 16, where we plot the Poisson corrected estimate of the mass-weighted number of halos $`\widehat{F}_{\mathrm{halo}}`$ as a function of distance from the galaxy centre and for each one of the boxes considered. We note that the mean number of halos per box in the outer galaxy is large, but that each contributes only a handful of particles. However, the inner galaxy is dominated by just a few halos making up most of the mass. This trend (of increasing number of halos with distance) is due to a form of mass segregation: heavy halos can sink by dynamical friction in short timescales to the centre of the newly formed system, whereas lighter halos, unable to decay quickly, can only deposit their mass in the outskirts.

### 5.3 The observed number of streams in the simulations

In general, a stream may be defined by those particles having the same orbital phase and coming from the same halo at some initial (high) redshift. The orbital phase of a particle may be determined, in principle, by counting how many times it has crossed the planes $`x=0,y=0,z=0`$. However, because of the limited number of outputs which are logarithmically spaced in the expansion factor, it is difficult to carry this through correctly. As in section 3, we therefore define a stream as a set of particles which

1. have been neighbours in all previous outputs;
2. are relatively close in space at the present time.

In practice, we make a coarse partition of the 3-dimensional space, whose elements are boxes of 15 kpc on a side. At each output, we check in which box of the partition any given particle is located. We tag the particle by this box ID, and by the IDs of those nearest to it. So for example, a given particle will generally not be located right at the centre of a box, but will be closer to one of its edges. Thus, this particle is assigned six different numbers (corresponding to the IDs of eight neighbouring boxes) for every output, as shown in Figure 17 for the two-dimensional analogue. We repeat this procedure for each output of our simulation. Particles which are on the same stream should have the same box tags or have been in neighbouring boxes at all previous outputs. This is almost equivalent to defining a box of 15 kpc on a side around each particle at each output, and finding which particles fall within that region. Our procedure is however, much more efficient computationally.

For the output corresponding to the present day, and for each particle in one of the 2 kpc boxes of the partition of the dark-matter halo, we find which other particles in that same box, satisfy condition (i) above (condition (ii) is automatically fulfilled). Thus at the end we obtain a link list, where particles have the same tags if they belong to the same stream.

In Figure 19 we plot the (Poisson corrected) mass-weighted mean mass per stream $`\widehat{\mu }`$ as a function of distance from the galaxy centre for each one of the boxes considered. This Figure shows that $`\widehat{\mu }`$ is generally smaller than unity. This happens when inside a box there are one or two streams with several particles, and the rest are one-particle structures. Clearly this is the case for the vast majority of the boxes. However, at larger distances from the Galactic centre, massive streams can be found.

In Figure 19 we plot the Poisson corrected estimate of the mass-weighted number of streams $`\widehat{F}`$ as a function of distance from the galaxy centre and for each one of the boxes considered. The thick grey line shows the median value of $`\widehat{F}`$, where we consider all boxes located at the same distance from the centre. For the “Solar neighbourhood” the estimate of the mass-weighted number of streams is roughly $`2\times 10^5`$ streams, showing that the local dark halo of the galaxy is extremely well-mixed. The variations in the number of streams at fixed distance from the galaxy centre are due to variations in the number of particles themselves. These result from the flattened shape of the galaxy halo.

### 5.4 An analytic estimate for the number of streams

It may at first seem counterintuitive to find a number of streams larger than the number of particles actually observed in a box. This effect is due to the small number of particles populating each stream, which boosts up the Poisson correction. It is fair to say that the number of streams we are determining is based on the actual presence of a few hundred streams, each detected with only two or three particles, in the inner 30 kpc region of the halo.

To ensure that the determination is meaningful we estimate how many streams we expect to find using the analytic prescription for the evolution of streams developed by HW. In Sec.3 we saw that this prescription does represent a reasonable approximation to the evolution of debris streams, even during the hierarchical build up of a CDM halo, such as that studied here.

To obtain an analytic estimate for the number of streams expected in a given box we proceed as follows:

1. We select 20 representative boxes at different distances from the halo centre.
2. We identify the progenitor halos of the individual particles in each of these boxes.
3. Among the particles belonging to the same progenitor halo, we select, when possible, 10 particles (for boxes out to $`r30`$ kpc). In all other cases, we choose 2 particles, or just 1, if this is all the halo contributes.
4. For each of these particles we integrate the orbit:

* forwards in time: starting from the position and velocity of the particle at $`z=2.4`$, using the NFW spherical potential of Eq.5 with $`\mathrm{\Phi }_0(z=2.4)`$;
    * backwards in time, until $`z=2.4`$, starting from the current position and velocity of the particle. In this case we use the NFW spherical potential approximation to the present-day halo.
5. Once the orbit is known we can compute the evolution (always forwards in time) of the density in the stream where the particle is located. This is done using the action-angle formalism of HW.
6. From the median (in time) density of a stream222As discussed in Sec.3.2 the density of a stream oscillates strongly in time. We choose to take the median density over all its possible values in the time interval since $`z=2.4`$, to avoid being dominated by a peak value just at the present time. we determine the number of streams from a given progenitor halo. The number of streams is given by the ratio of the mean (coarse-grained) density of the debris to the actual density of the stream where the particle $`i`$ is located, i.e.

$$N_{\mathrm{stream},i}^{\mathrm{halo}}=\frac{M^{\mathrm{halo}}}{V_{\mathrm{orb}}}\frac{1}{\rho _{\mathrm{stream},i}(t)}.$$
(10)
Here $`M^{\mathrm{halo}}=4\pi /3r_{200}^3\rho ^{\mathrm{halo}}`$, $`V_{\mathrm{orb}}4\pi /3r_{\mathrm{apo},i}^3`$, and $`\rho _{\mathrm{stream},i}(t)=\rho _0f_i(t)`$, where we assume $`\rho _0\rho ^{\mathrm{halo}}`$. Therefore we can express the number of streams as

$$N_{\mathrm{stream},i}^{\mathrm{halo}}=\left(\frac{r_{200}^{\mathrm{halo}}}{r_{\mathrm{apo},i}}\right)^3\frac{1}{f_i(t)}.$$
(11)
7. Having obtained the number of streams for each particle, we derive the average number of streams for each halo from all $`n_p^{\mathrm{halo}}`$ particles that fall within the same box and belonged to this same halo:

$$N_{\mathrm{stream}}^{\mathrm{halo}}=\frac{1}{n_\mathrm{p}^{\mathrm{halo}}}\underset{i=1}{\overset{n_\mathrm{p}^{\mathrm{halo}}}{}}N_{\mathrm{stream},i}^{\mathrm{halo}}.$$
(12)
8. The total number-weighted number of streams in a box is now obtained by adding over all halos contributing to the box, where the weights are given by the number of particles from each halo in the box $`n_\mathrm{p}^{\mathrm{halo}}`$:

$$N_{\mathrm{stream}}=\underset{\mathrm{halos}\mathrm{in}\mathrm{box}}{}\left(n_\mathrm{p}^{\mathrm{halo}}N_{\mathrm{stream}}^{\mathrm{halo}}\right)/\underset{\mathrm{halos}\mathrm{in}\mathrm{box}}{}n_\mathrm{p}^{\mathrm{halo}}.$$
(13)

In Figure 20 we plot our estimate for the number of streams as a function of distance from the Galactic centre. The solid curve corresponds to the estimate obtained by integrating the orbit backwards in time, from the present day conditions. The dotted curve, on the other hand, corresponds to the orbit integration performed forwards in time from $`z=2.4`$. The geometric mean of these two estimates is the dashed curve in the same figure. The predicted number of streams given by the geometric mean is in reasonable agreement with what was shown in Figure 19.

### 5.5 Characteristics of the streams

Another interesting quantity which we can calculate for each box is the characteristic mass-weighted dispersion within a stream:

$$\sigma _{\mathrm{stream}}=\sqrt{\frac{1}{N_T}\underset{k=1,L}{}\frac{_{ij}|𝐯_i𝐯_j|^2}{N_k1}}$$
(14)
where the $`ij`$ sum runs over all particle pairs in the $`k`$-th stream in a given box, $`N_k`$ is the number of particles in this stream in this box, and as before, $`N_T`$ is the total number of particles in the box. This mass-weighted velocity dispersion is expected to be relatively small in view of our results of the phase-space evolution of halo debris. We plot this quantity in Figure 21 as a function of distance from the galaxy centre for each one of the boxes considered. For our boxes of 2 kpc on a side we find a typical velocity dispersion of 1 km s-1 at the position of the “Sun”. The values observed are probably upper limits to the typical velocity dispersion in a stream, since they could only be measured for the densest streams (with two or more particles), and so are biased towards high velocity dispersions.

## 6 Discussion

We have studied the phase-space evolution of debris from the progenitors that merge to build up a dark-matter halo in a $`\mathrm{\Lambda }`$CDM cosmology. Our analysis has shown that the debris streams originating in progenitors of different sizes and orbital characteristics all behave in a similar way: with velocity dispersions and local space density decreasing in time. The evolution of the debris streams that we were able to follow until the present time is consistent with phase mixing. Even for halos that we could not follow for a very long time – because of their smaller initial number of particles or their shorter orbital timescales– we find the debris to show similar behaviour. On the scales of the fine-grained distribution function, mixing is apparently not strongly chaotic. On the contrary, the phase-space evolution appears to be quite organised and simple, very similar to the mixing of streams orbiting in an idealised static potential.

In principle, a dark-matter halo formed in a $`\mathrm{\Lambda }`$CDM cosmology is not a smooth entity. Not only do dark-matter halos contain a large number of dark satellites, they also have large amounts of substructure in the form of streams. We predict, however, that dark matter in the Solar neighbourhood should be clumped in a few hundred thousand streams, producing a velocity ellipsoid which is close to a multivariate Gaussian (Helmi, White & Springel 2002). These streams have their origin in the different halos that merged to form the dark halo of the Galaxy. Most of these halos give rise to a large number of intersecting streams in the inner Galaxy.

Determining the characteristics of these hundred thousand streams is difficult even with the high-resolution simulation used here. We are mostly limited by the number of particles. Although the simulation as a whole has 66 million particles, inside a 4 kpc box centred on the “Sun”, we find only a couple of thousand particles. In such boxes we find on average two hundred streams with more than one particle, and typically each has only two particles! The internal stream velocity dispersion that we measure at the present-day is extremely small, of the order of only 1 km s-1.

It is encouraging that we find reasonable agreement between the behaviour of debris streams in static potentials and that observed in this high resolution simulation. This implies that our earlier estimates of the number of stellar streams in the vicinity of the Sun (Helmi & White 1999) may indeed apply. Since the initial stellar phase-space distribution in the progenitor objects was probably of lower dimensionality than we here assume – stars tend to form in disks so that the distribution in at least two of the six phase-space coordinates collapses – star streams may be colder than the dark matter streams we have analysed, and so may be more easily distinguishable. It is also interesting to note that because the material that ends up populating the inner galaxy was already in place 10 Gyr ago, the oldest stars are predicted to be near the galactic centre (White & Springel 2000). Because this material comes from only a few objects, one might expect the stellar populations to be quite homogeneous, although this of course depends on whether the stars themselves formed in these few massive objects, or whether they were accreted into these objects in the first place.

A crude estimate of the stellar content of a stream in the Solar neighbourhood can be obtained as follows. Let us first estimate the “mass-to-light ” ratio per particle $`[M_{\mathrm{DM}}/L_{}]`$ as the ratio of dark-matter mass to stellar halo light enclosed in a shell of thickness 2 kpc at the location of the Sun. We use the NFW profile of Eq.(3) for the dark-matter distribution. For the stellar distribution we assume a density profile $`\rho _{}(r)=\rho _{}r_{}^2/(r^3+r_c^3)`$, where $`\rho _{}3\times 10^4\mathrm{M}_{}`$ kpc-3, and where $`r_c`$ is the core radius, and is much smaller than $`r_{}`$. We obtain for the Solar Neighbourhood $`[M_{\mathrm{DM}}/L_{}]625`$, for $`\mathrm{{\rm Y}}_{}=2.5\mathrm{{\rm Y}}_{}`$. For one of our most massive streams (with 3 particles), this would imply an average $`0.54\mathrm{M}_{}`$ in stars in a sphere of 100 pc radius centred on the Sun. Assuming a Salpeter initial mass function, and down to an absolute magnitude of $`M_V=7^m`$ (Bergbusch & VandenBerg 1992), i.e. $`V=12^m`$ at a distance of 100 pc, this corresponds to approximately 2 stars per stream. On the other hand, streams originate from very localised regions of phase-space in the progenitor objects, so it seems more likely that only about one stream in 250 has any stars at all, that these streams have $`\mathrm{{\rm Y}}_{}=2.5`$ characteristic of the “stellar” regions of their progenitors, and that all the other streams are dark. In this case a massive, “luminous” stream might contain $`135\mathrm{M}_{}`$ in stars in a sphere of 100 pc radius centred on the Sun, implying approximately $`450`$ stars down to $`M_V=7^m`$ for a Salpeter IMF.

More reliable estimates of the stellar content and number of streams, as well as considerably more insight into the properties of the Galactic stellar halo would be obtained by combining semi-analytic techniques (e.g Kauffmann et al. 1993) with high-resolution simulations (see, for example, Springel et al. 2001). This would enable one to predict trends in the chemical composition, age, spatial distribution and kinematics of halo stars as a function of position throughout the Galaxy.

## Acknowledgments

We wish to thank the referee for the positive remarks which helped us improve this manuscript.