text:
# Do Chaotic Trajectories Care About Self-Similarity?

\[

## Abstract

We investigate the relation between the chaotic dynamics and the hierarchical phase-space structure of generic Hamiltonian systems. We demonstrate that even in ideal situations when the phase space is dominated by an exactly self-similar structure, the long-time dynamics is not dominated by this structure. This has consequences for the power-law decay of correlations and Poincaré recurrences.

\]

Generic Hamiltonian systems are neither integrable nor chaotic , but rather exhibit a mixed phase space, where regular and chaotic regions coexist. Each island of regular motion is surrounded by infinitely many chains of smaller islands. As the same holds for any of these smaller islands a very complex hierarchical phase-space structure is found for generic Hamiltonian systems, which is well understood and nowadays appears in textbooks on classical mechanics.

The dynamical properties, however, are still unclarified. The most fundamental statistical quantity for characterizing dynamics is the decay of correlations in time. It determines transport properties and is directly related to the distribution of Poincaré recurrences $`P(t)`$, which is the probability to return to a given region in phase space with a recurrence time larger than $`t`$. This probability decays on average like a power-law

$$P(t)t^\gamma ,$$
(1)
due to the trapping of chaotic trajectories in the hierarchically structured vicinity of islands of regular motion. The power-law decay is a universal property of Hamiltonian systems. It has dramatic consequences for transport (anomalous diffusion) and quantum mechanics (conductance fluctuations and eigenfunctions), which sensitively depend on the value of $`\gamma `$. The exponent $`\gamma `$, as determined from finite time numerical experiments, seems to be non-universal, varies with system and parameter, and typically ranges between $`1`$ and $`2.5`$. It is a fundamental question of Hamiltonian chaos, how the exponent $`\gamma `$ of the dynamics is related to the structure of the hierarchical phase space.

Recently, it was argued by Chirikov and Shepelyansky that for asymptotically large times the exponent is independent of the specific system and parameter and is given by the universal exponent $`\gamma =3`$. Their arguments are based on the universal presence of critical tori in phase space and are supported by a numerical investigation of the kicked rotor at kicking strength $`K=K_c=0.97163540631`$. At this parameter value the golden torus is critical, i.e., it can be destroyed by an arbitrarily small perturbation. The self-similar vicinity of the critical golden torus \[see Fig. 1(a)\] has been studied using renormalization methods and the asymptotic value $`\gamma =3`$ for the power-law decay of $`P(t)`$ was predicted long time ago . The fact that it has never been observed led to the speculation that the universal decay should appear for larger times. In Ref.  a numerical approach allowed to estimate the onset of this decay in agreement with the presented data for $`P(t)`$.

In addition to the sticking of trajectories in the vicinity of critical tori, the trapping of trajectories in island-around-island structures has been studied . Zaslavsky et al.  showed that for the kicked rotor at $`K=K^{}=6.908745`$ the phase space possesses an island-around-island structure of sequence $`3888\mathrm{}`$ \[see Fig. 1(c)\]. They used this self-similarity to derive the trapping exponent $`\gamma =2.25`$ by renormalization arguments, which was recently supported by the numerics.

In fact, these renormalization approaches for single self-similar phase-space structures can be considered as special cases of the more general binary tree model by Meiss and Ott. In this model a chaotic trajectory can at any stage of the tree either go to a boundary circle (level scaling) or to the island-around-island structure (class scaling). The universal coexistence of the two routes of renormalization at any stage led to the exponent $`\gamma =1.96`$. In contrast, the recent findings claim that just one of these scalings is relevant for the trapping mechanism: While in Ref.  it is argued that universally (and in particular for $`K=K_c`$) the level scaling should dominate, in Ref.  it is claimed that for $`K=K^{}`$ the class scaling describes the trapping mechanism.

In order to clarify these contradictions, we numerically investigate $`P(t)`$ for the kicked rotor at $`K_c`$ and $`K^{}`$ for times larger than in previous studies. We find strong deviations from the predictions of the renormalization theories that only consider a single self-similar phase-space structure. In addition, our numerical approach allows to analyze where chaotic trajectories are trapped in phase space. For large times the majority of trajectories is not trapped in those phase-space regions, that are described by the simple renormalization theories. We thereby reveal the mistaken assumption of these theories. In particular for $`K=K_c`$, the self-similar vicinity of the critical golden torus does not dominate the trapping mechanism for large times and thus even in this ideal situation the proposed universal exponent $`\gamma =3`$ is not found. For $`K=K^{}`$, the majority of long trapped trajectories does not follow the self-similar island-around-island structure, which leads to a smaller exponent than predicted. Although the phase space in both cases is dominated by exactly self-similar structures, our analysis shows that they do not dominate the dynamics.

We use the standard map (kicked rotor) defined by

$$q_{n+1}=q_n+p_n\mathrm{mod}\mathrm{\hspace{0.17em}2}\pi p_{n+1}=p_n+K\mathrm{sin}q_{n+1},$$
(2)
which has a $`2\pi `$-periodic phase space in $`p`$ direction. We concentrate on two parameters: (i) The dynamics for $`K=K_c`$ is bounded in $`p`$ direction by the golden torus, which is critical . The route towards the critical golden torus is determined by the principal resonances given by the approximants of the golden mean $`\sigma =(\sqrt{5}1)/2`$ and the scaling has been analyzed in detail . The dynamics along this route was described by a Markov chain leading to $`\gamma =3.05`$  and alternatively via the scaling of the local diffusion rate leading to $`\gamma =3`$ . (ii) The phase space for $`K=K^{}`$ consists of two small accelerator modes embedded in an otherwise completely chaotic phase space. Each mode shows an island-around-island structure of sequence $`3888\mathrm{}`$ . This exact scaling relation was used to predict the exponent $`\gamma =2.25`$.

In order to check the predictions for $`P(t)`$, in case (i) we start several long trajectories initially located near the unstable fixed point $`(q,p)=(0,0)`$. We measure the times $`\tau `$ for which an orbit stays close to the critical torus by monitoring successive crossings of the line $`p=0`$, as was also done in Ref. . In case (ii) we start many trajectories at four different regions in the chaotic part of phase space. Whenever a chaotic trajectory is trapped to one of the island structures, it follows the dynamics of the accelerator mode and jumps to the neighboring unit cell in $`p`$ direction. We measure the time $`\tau `$ it continuously jumps one unit cell per iteration in the same direction. From the set of trapping times $`\tau `$ one determines the fraction $`\stackrel{~}{P}(t)`$ of orbits with $`\tau t`$. This quantity decays with the same power-law exponent as the Poincaré recurrences $`P(t)`$ and was chosen for numerical convenience. The total computer time corresponds to (i) $`1510^{12}`$ and (ii) $`810^{12}`$ iterations of the standard map. We have checked if our statistical data for large times are sensitive to the unavoidable finite numerical precision, by comparing data for double ($``$16 significant digits) and quadruple ($``$32 digits) precision. We found no difference and present a combination of both data sets in Fig. 2.

In Fig. 2 we compare our numerical findings for $`P(t)`$ for case (i) with the prediction $`P(t>10^7)=3.910^{12}t^3`$ extracted from Ref. . This power law is not compatible with our data, even though we are in the time regime, where it should be observable according to Ref. . For $`10^5t10^9`$ we rather see an exponent $`\gamma 1.9`$. For case (ii) with $`K=K^{}`$, we find a power-law decay of $`P(t)`$ with $`\gamma =1.85`$ (Fig. 3) for various starting conditions  contradicting the renormalization prediction $`\gamma =2.25`$ .

Our numerical results raise the question why the renormalization theories predict the wrong exponents. This is particularly surprising since the structure of the phase space for the specific parameters $`K_c`$ and $`K^{}`$ is dominated by one exactly self-similar hierarchy. All these approaches rely on the assumption that this single self-similar phase-space structure also dominates the long-time trapping of chaotic trajectories. We will now check this assumption.

To this end we calculate the density in phase space for trajectories, that are trapped for long times. For case (i) we show two examples in Fig. 1(a,b). Figure 1(a) shows a trajectory of length $`t510^7`$ that follows the route to the critical golden torus up to the principal resonance with winding number $`55/144`$. This is consistent with the renormalization theory according to the data presented in Ref. . In contrast, the trajectory shown in Fig. 1(b) approaches the critical torus only up to the resonance $`3/8`$ and is predominantly trapped around a non-principal resonance. This trajectory is not captured by the renormalization theory, since it has the same length as the trajectory in Fig. 1(a) and therefore should be trapped around the $`55/144`$ resonance or one of its neighbors. In order to quantify this observation we introduce the fraction $`f(t)`$ of trajectories with trapping time $`t`$ that follow the route of renormalization. Numerically, we determine $`f(t)`$ by considering trajectories with trapping times in an interval around $`t`$. We classify these trajectories, as was done in Fig. 1, according to their phase-space densities. For $`t>210^7`$ this fraction decreases to zero (Fig. 2, inset). This clearly demonstrates that for large times the majority of trajectories is not trapped in the self-similar phase-space regions approaching the critical golden torus. The assumption of the renormalization theories leading to $`\gamma 3`$ is violated and thus they are not applicable for predicting $`P(t)`$. From the ratio of the predicted $`P(t)t^3`$ for the trajectories trapped in the self-similar phase-space structure and the observed $`P(t)t^{1.9}`$ we expect their fraction to decay as $`f(t)t^{1.1}`$ for large times. This is confirmed by our numerical data in the inset of Fig. 2. We find that the majority of trajectories is trapped around non-principal resonances, which is in agreement with the binary tree model . This is in strong contrast to the conclusions of Ref.  that are based on the computation of exit times from the vicinity of unstable fixed points of principal resonances. As for the investigation of local diffusion rates also the analysis of the mean exit time yields information only about trajectories trapped in the region studied. One cannot conclude, however, that these trajectories dominate the global trapping mechanism, as is convincingly shown by our numerics.

We carry out the same analysis for case (ii), $`K=K^{}`$. In Figure 1(c,d) we show the phase-space densities of two long trajectories. Although both trajectories have the same length, only the trajectory in Fig. 1(c) follows the self-similar island-around-island structure, while the trajectory shown in Fig. 1(d) is trapped around other islands. The fraction $`f(t)`$ of trajectories following the route of renormalization decays to zero for large times (Fig. 3, inset). This decay is well described by the estimate $`f(t)t^{2.25}/t^{1.85}t^{0.4}`$, i.e., the ratio of the predicted $`P(t)t^{2.25}`$ and the observed $`P(t)t^{1.85}`$. This shows that the renormalization theory for the island-around-island structure is not capable of explaining $`P(t)`$. It should be noted that this difference is not caused by the effect, that the finite precision of $`K^{}`$ eventually leads to a breakdown of the self-similarity on very small scales.

In conclusion, our analysis shows that even in the presence of an exactly self-similar phase-space structure it is not sufficient to describe the trapping mechanism of chaotic trajectories by only this structure, as was recently claimed in the literature. We find that additional island structures may dominate the trapping mechanism for large times and thus affect the power-law decay of $`P(t)`$. Our analysis supports qualitatively the tree model by Meiss and Ott, which allows for the coexistence of two routes of renormalization at any stage. Quantitatively, we find slightly smaller exponents. It remains an open question if there exists a universal asymptotic exponent for the trapping of chaotic trajectories in Hamiltonian systems.

We thank R. Fleischmann for helpful discussions. M.W. acknowledges financial support by an EMBO fellowship.


format: text

text:
# 1 Introduction

## 1 Introduction

The energy loss of heavy quarks in QCD matter has excited considerable interest. In particular, collisional energy loss in quark–gluon plasma has been evaluated  by the methods of finite–temperature QCD (see  for a comprehensive review). However, the energy loss of fast partons in medium is dominated by radiation of gluons  and, to the best of our knowledge, a consistent treatment of radiative energy loss of heavy quarks in QCD matter is still lacking. Practical importance of this issue has become clear from a number of studies  which showed a strong dependence of charmed hadron and lepton spectra on the assumed heavy quark energy loss.

In this letter we follow the recent analysis of the quenching of inclusive particle spectra in QCD medium  to compare the yields of hadrons from light and heavy quarks produced in hard interactions in heavy ion collisions. We find that the specific pattern of gluon bremsstrahlung off heavy quarks (the dead cone effect) results in a much smaller heavy quark quenching, leading to a relative enhancement of heavy particle production in heavy ion collisions.

## 2 Medium induced gluon radiation – a brief summary

We first recall the basic features of gluon radiation caused by propagation of a fast parton (quark) through QCD medium.

As was pointed out in , the accompanying radiation is determined by multiple rescattering of the radiated gluon in the medium. The gluon, during its formation time

$$t_{form}\frac{\omega }{k_{}^2},$$
(1)
accumulates a typical transverse momentum

$$k_{}^2\mu ^2\frac{t_{form}}{\lambda },$$
(2)
with $`\lambda `$ the mean free path and $`\mu ^2`$ the characteristic momentum transfer squared in a single scattering. This is the random walk pattern with an average number of scatterings given by the ratio $`t_{form}/\lambda `$.

Combining (2) and (1) we obtain

$$N_{coh}=\frac{t_{form}}{\lambda }=\sqrt{\frac{\omega }{\mu ^2\lambda }}$$
(3)
describing the number of scattering centres which participate, coherently, in the emission of the gluon with a given energy $`\omega `$. For sufficiently large gluon energies, $`\omega >\mu ^2\lambda `$, when the coherent length exceeds the mean free path, $`N_{coh}>1`$. In this situation the standard Bethe-Heitler energy spectrum per unit length describing independent emission of gluons at each centre gets suppressed:

$$\frac{dW}{d\omega dz}=\left(\frac{dW}{d\omega dz}\right)^{\text{BH}}\frac{1}{N_{coh}}=\frac{\alpha _\text{s}C_R}{\pi \omega \lambda }\sqrt{\frac{\mu ^2\lambda }{\omega }}=\frac{\alpha _\text{s}C_R}{\pi \omega }\sqrt{\frac{\widehat{q}}{\omega }}.$$
(4)
Here $`C_R`$ is the “colour charge” of the parton projectile ($`C_R=C_F=\frac{N_c^21}{2N_c}=4/3`$ for the quark case we are interested in).

In (4) we have substituted the characteristic ratio $`\mu ^2/\lambda `$ by the so-called gluon transport coefficient

$$\widehat{q}\rho \frac{d\sigma }{dq^2}q^2𝑑q^2,$$
(5)
which is proportional to the density $`\rho `$ of the scattering centres in the medium and describes the typical momentum transfer in the gluon scattering off these centres.

The transport coefficient for cold nuclear matter was expressed in  as

$$\widehat{q}\frac{4\pi ^2\alpha _sN_c}{N_c^21}\rho [xG(x,Q^2)],$$
(6)
with $`\rho 0.16\mathrm{fm}^3`$ the average nuclear density and $`[xG(x)]`$ the gluon density in a nucleon. The latter should be evaluated at the characteristic scale $`Q^2\widehat{q}L`$, which makes (6), strictly speaking, an equation for $`\widehat{q}`$. Given a slow logarithmic $`Q^2`$ dependence of the gluon density, this equation can be solved iteratively. Choosing for the sake of estimate $`L5`$ fm results in $`Q0.5\mathrm{GeV}`$. Taking at this scale $`\alpha _s0.5`$ and $`[xG(x)]1`$ (at $`x<0.1`$), yields

$$\widehat{q}_{\text{cold}}0.01\mathrm{GeV}^38\rho .$$
(7)
This value is an agreement with the result of the analysis of the gluon $`p_{}`$ broadening from the experimental data on $`J/\psi `$ transverse momentum distributions , which in the present notation yielded

$$\widehat{q}=(9.4\pm 0.7)\rho .$$
(8)
An estimate  for a hot medium based on perturbative treatment of gluon scattering in quark–gluon plasma with $`T250`$ MeV resulted in the value of the gluon transport coefficient of about factor twenty larger than (7):

$$\widehat{q}_{\text{hot}}0.2\mathrm{GeV}^3\mathrm{\hspace{0.25em}20}\widehat{\mathrm{q}}_{\text{cold}}.$$
(9)
Multiplying (4) by the length $`L`$ of the medium traversed,222For the sake of simplicity we assume here that the medium is static and uniform. we arrive at the following expression for the inclusive energy distribution of gluons radiated by a quark:

$$\frac{dW}{d\omega }\frac{\alpha _\text{s}C_F}{\pi \omega }\sqrt{\frac{\omega _1}{\omega }},\omega <\omega _1\widehat{q}L^2.$$
(10)
A more accurate treatment shows that the semi-quantitative estimate (4) is valid for $`\omega \begin{array}{c}<\hfill \\ \hfill \end{array}\omega _1/2`$, while at very large gluon energies $`\omega \omega _1`$ the gluon yield is very small and decreases as $`dI/d\omega \omega ^3`$ (for details, see ).

The fact that the medium induced radiation vanishes for $`\omega >\omega _1`$ has a simple physical explanation, as according to (3) the formation time of such gluons starts to exceed the length of the medium:

$$t_{form}=\lambda \sqrt{\frac{\omega }{\mu ^2\lambda }}=\sqrt{\frac{\omega }{\widehat{q}}}=L\sqrt{\frac{\omega }{\omega _1}}>L.$$

Another important feature of medium induced radiation is the relation between the transverse momentum and the energy of the emitted gluon. Indeed, from (1) and (2) (see also (5)) we derive

$$k_{}^2\sqrt{\widehat{q}\omega }.$$
(11)
This means that the angular distribution of gluons with a given energy $`\omega `$ is concentrated at a characteristic energy- (and medium-) dependent emission angle

$$\theta \frac{k_{}}{\omega }\left(\frac{\widehat{q}}{\omega ^3}\right)^{1/4}.$$
(12)
## 3 Radiation off heavy quarks: dead cone

Gluon bremsstrahlung off a heavy quark differs from the case of a massless parton (produced in a process with the same hardness scale) in one respect: gluon radiation is suppressed at angles smaller than the ratio of the quark mass $`M`$ to its energy $`E`$. Indeed, the distribution of soft gluons radiated by a heavy quark is given by

$$dP=\frac{\alpha _\text{s}C_F}{\pi }\frac{d\omega }{\omega }\frac{k_{}^2dk_{}^2}{(k_{}^2+\omega ^2\theta _0^2)^2},\theta _0\frac{M}{E},$$
(13)
where the strong coupling constant $`\alpha _s`$ should be evaluated at the scale determined by the denominator of (13). Equating, in the small-angle approximation, $`k_{}`$ with $`\omega \theta `$ we conclude that the formula (13) differs from the standard bremsstrahlung spectrum

$$dP_0\frac{\alpha _\text{s}C_F}{\pi }\frac{d\omega }{\omega }\frac{dk_{}^2}{k_{}^2}=\frac{\alpha _\text{s}C_F}{\pi }\frac{d\omega }{\omega }\frac{d\theta ^2}{\theta ^2}$$
(14)
by the factor

$$dP_{\text{HQ}}=dP_0\left(1+\frac{\theta _0^2}{\theta ^2}\right)^2$$
(15)
This effect is known as the “dead cone” phenomenon. Suppression of small-angle radiation has a number of interesting implications, such as perturbative calculability of (and non-perturbative $`\mathrm{\Lambda }/M`$ corrections to) heavy quark fragmentation functions , multiplicity and energy spectra of light particles accompanying hard production of a heavy quark .

In the present context we should compare the angular distribution of gluons induced by the quark propagation in the medium with the size of the dead cone. To this end, for the sake of a semi-quantitative estimate, we substitute the characteristic angle (12) into the dead cone suppression factor (15) and combine it with the radiation spectrum (4) to arrive at

$$I(\omega )=\omega \frac{dW}{d\omega }=\frac{\alpha _sC_F}{\pi }\sqrt{\frac{\omega _1}{\omega }}\frac{1}{(1+(\mathrm{}\omega )^{3/2})^2},$$
(16)
where

$$\mathrm{}\widehat{q}^{1/3}\left(\frac{M}{E}\right)^{4/3}.$$
(17)
The suppression factor (16) is displayed in Fig. 1 for two quark energies ($`x=\omega /p_{}`$).

To see whether the finite quark mass essentially affects the medium induced gluon yield, we need to estimate the product $`\mathrm{}\omega `$ for the maximal gluon energy $`\omega \omega _1`$ to which the original distribution (4) extends:

$$\mathrm{}\omega _1=\widehat{q}^{1/3}\left(\frac{M}{E}\right)^{4/3}\widehat{q}L^2=\left(\frac{E_{\text{HQ}}}{E}\right)^{4/3},E_{\text{HQ}}M\sqrt{\widehat{q}L^3}.$$
(18)
This shows that the quark mass becomes irrelevant when the quark energy exceeds the characteristic value $`E_{\text{HQ}}`$ which depends on the size of the medium and on its “scattering power” embodied into the value of the transport coefficient.

Which regime is realized in the experiments on heavy quark production in nuclear collisions? Taking $`M=1.5\mathrm{GeV}`$ for charm quarks and using the values (7) and (9) we estimate

$`E_{\text{HQ}}`$ $`=`$ $`\sqrt{\widehat{q}_{\text{cold}}}L^{3/2}M20\mathrm{Ge}\mathrm{V}\left({\displaystyle \frac{L}{5\mathrm{fm}}}\right)^{3/2},`$ (19)
$`E_{\text{HQ}}`$ $`=`$ $`\sqrt{\widehat{q}_{\text{hot}}}L^{3/2}M92\mathrm{Ge}\mathrm{V}\left({\displaystyle \frac{L}{5\mathrm{fm}}}\right)^{3/2},`$ (20)
for the cold and hot matter, respectively. We observe that for the transverse momentum (energy) distributions of heavy mesons the inequality $`EE_{\text{HQ}}`$ always holds in practice, especially for the hot medium. We thus conclude that the pattern of medium induced gluon radiation appears to be qualitatively different for heavy and light quarks in the kinematical region of practical interest.

## 4 Quenching

The issue of in-medium quenching of inclusive particle spectra was recently addressed in . The $`p_{}`$ spectrum is given by the convolution of the transverse momentum distribution in an elementary hadron–hadron collision, evaluated at a shifted value $`p_{}+ϵ`$, with the distribution $`D(ϵ)`$ in the energy $`ϵ`$ lost by the quark to the medium-induced gluon radiation:

$$\frac{d\sigma ^{med}}{dp_{}^2}=𝑑ϵD(ϵ)\frac{d\sigma ^{vac}}{dp_{}^2}(p_{}+ϵ)\frac{d\sigma ^{vac}}{dp_{}^2}(p_{})Q(p_{}),$$
(21)
with $`Q(p_{})`$ the medium dependent quenching factor. The two facts, namely that in the essential region $`ϵp_{}`$ and that the vacuum cross section is a steeply falling function, allow one to simplify the calculation of the quenching factor $`Q`$ by adopting the exponential approximation for the $`ϵ`$-integral in (21):

$$Q(p_{})𝑑ϵD(ϵ)\mathrm{exp}\left\{\frac{ϵ}{p_{}}\right\},\frac{d}{d\mathrm{ln}p_{}}\mathrm{ln}\left[\frac{d\sigma ^{vac}}{dp_{}^2}(p_{})\right].$$
(22)
This integral results in the Mellin moment of the quark distribution,

$$Q(p_{})=\stackrel{~}{D}(\nu )=\mathrm{exp}\left[\nu _0^{\mathrm{}}𝑑\omega N(\omega )e^{\nu \omega }\right],\nu =\frac{}{p_{}},$$
(23)
where $`N(\omega )`$ is the integrated gluon multiplicity defined according to (see  for details)

$$N(\omega )_\omega ^{\mathrm{}}𝑑\omega ^{}\frac{dW(\omega ^{})}{d\omega ^{}}.$$
(24)
Here we present an approximate evaluation based on a simplified energy spectrum (16). Integration of (16) leads to the following expression for the gluon multiplicity (24):

$$N(\omega )=\frac{2\alpha _sC_F}{\pi }\sqrt{\mathrm{}\omega _1}K(\mathrm{}\omega ),$$
(25)
where

| $`K(x)`$ | $`{\displaystyle \frac{1}{\sqrt{x}}}{\displaystyle \frac{2\pi }{3\sqrt{3}}}+{\displaystyle \frac{x}{3\left(1+x^{3/2}\right)}}+{\displaystyle \frac{4}{3\sqrt{3}}}\mathrm{arctan}\left({\displaystyle \frac{1+2\sqrt{x}}{\sqrt{3}}}\right)`$ |
| --- | --- |
|  | $`+{\displaystyle \frac{2}{9}}\mathrm{ln}\left[{\displaystyle \frac{1\sqrt{x}+x}{(1+\sqrt{x})^2}}\right].`$ |

(26)
The function $`K(x)`$ has essential support in the region $`\mathrm{}\omega =x1`$, that is where the gluon yield is not suppressed by the quark mass (dead cone) effect (cf. (16)). To find the quenching factor we need to evaluate its Mellin transform according to (23), with $`\nu `$ given in (35). The essential energies $`\omega `$ in the multiplicity function in (23) are then restricted to $`\omega \begin{array}{c}<\hfill \\ \hfill \end{array}1/\nu `$. It is easy to check that, due to the bias effect (large value of $`10`$), this restriction is stronger than the quark mass energy cutoff:

$$\frac{\nu }{\mathrm{}}=\frac{}{p_{}}\widehat{q}^{1/3}\left(\frac{E}{M}\right)^{4/3}\left(\frac{p_{}\widehat{q}}{M^4}\right)^{1/3}>\mathrm{\hspace{0.25em}2}$$
(27)
for $`p_{}>5`$ GeV, even for cold matter. Therefore, in this practically interesting domain, we can approximate

$$N(\omega )\frac{2\alpha _sC_F}{\pi }\sqrt{\omega _1}\left(\frac{1}{\sqrt{\omega }}\frac{8\pi \sqrt{\mathrm{}}}{9\sqrt{3}}\right)\left[\mathrm{\hspace{0.17em}1}+𝒪\left(\mathrm{}\omega \right)\right],$$
(28)
and get

$$\stackrel{~}{D}(\nu )\mathrm{exp}\left[\frac{2\alpha _sC_F}{\pi }\sqrt{\widehat{q}}L\left(\sqrt{\pi \nu }\frac{8\pi \sqrt{\mathrm{}}}{9\sqrt{3}}\right)\right].$$
(29)
The use of (23) furnishes our final result:

$$Q_H(p_{})\mathrm{exp}\left[\frac{2\alpha _sC_F}{\sqrt{\pi }}L\sqrt{\widehat{q}\frac{_H}{p_{}}}+\frac{16\alpha _sC_F}{9\sqrt{3}}L\left(\frac{\widehat{q}M^2}{M^2+p_{}^2}\right)^{1/3}\right].$$
(30)
The first term in the exponent in (30) represents the quenching of the transverse momentum spectrum which is universal for the light and heavy quarks,

$$Q_L(p_{})\mathrm{exp}\left[\frac{2\alpha _sC_F}{\sqrt{\pi }}L\sqrt{\widehat{q}\frac{_L}{p_{}}}\right]$$

(modulo the difference of the $``$ parameters determined by the $`p_{}`$ distributions in the vacuum). The second term is specific for heavy quarks. It has a positive sign, which means that the suppression of the heavy hadron $`p_{}`$ distributions is always smaller than that for the light hadrons. This is a straightforward consequence of the fact that the heavy quark mass suppresses gluon radiation. At very high transverse momenta, both terms vanish – this is in accord with the QCD factorization theorem, stating that the effects of the medium should disappear as $`p_{}\mathrm{}`$. How fast this regime is approached depends, however, on the properties of the medium encoded in the value of the transport coefficient $`\widehat{q}`$ and in the medium size $`L`$.

Constructing the ratio of the quenching functions, we estimate the heavy-to-light enhancement factor as

$$\frac{Q_H(p_{})}{Q_L(p_{})}\mathrm{exp}\left[\frac{16\alpha _sC_F}{9\sqrt{3}}L\left(\frac{\widehat{q}M^2}{M^2+p_{}^2}\right)^{1/3}\right].$$
(31)
This simple expression provides a reasonably good approximation to more accurate quantitative results presented below.

## 5 Quantitative estimates

To perform a quantitative calculation of the transverse momentum distributions, our considerations have to be combined with a realistic treatment of nuclear geometry, including its centrality dependence, gluon multiple interactions in the initial state (the “Cronin effect”), and a model for the time evolution of QCD matter in the final state. Such an analysis is beyond the scope of the present paper. Nevertheless, we provide some semi-quantitative illustrations of the expected consequences of our results, based on a simplified model of a static uniform medium and a fixed quark path length $`L`$.

In what follows we use the full inclusive soft gluon emission spectrum off a massless quark from , taking into account the first non-soft correction $`𝒪\left(x\right)`$:

$$I(x)=(1x)\frac{C_F\alpha _\text{s}}{\pi }\mathrm{ln}\left(\mathrm{cosh}^2\left(\sqrt{\frac{\omega _1}{4\omega }}\right)\mathrm{sin}^2\left(\sqrt{\frac{\omega _1}{4\omega }}\right)\right),x\frac{\omega }{E},$$
(32)
which, for the heavy quark case, is supplied with the dead-cone suppression factor (15). To illustrate the “dead cone” effect, in Fig. 2 a comparison of the inclusive one-gluon emission spectra off light and charm quarks is shown, for hot medium with $`L=5\mathrm{fm}`$.

As we have discussed above, the quenching of $`p_{}`$ distributions of heavy hadrons caused by QCD matter is much weaker than for pions. This is because the gluon cloud around the heavy quark is “truncated” by the large quark mass, in a way which depends on the properties of the medium. This interesting effect can be illustrated by the transverse momentum dependence of the ratio of hadrons originating from the fragmentation of heavy and light quarks (for example, $`D/\pi `$ ratio) in heavy ion collisions.

To estimate this ratio we need to know the behaviour of the vacuum spectra. For light hadrons (which we assume to originate from light quarks333we remark that the quenching of hadrons produced in the fragmentation of primary gluons contains an additional factor $`9/42`$ in the exponent) the parameterization of the $`p_{}`$ distribution which describes the first RHIC hadroproduction data for $`p_{}`$ up to $`6`$ GeV :

$$\frac{d\sigma }{dp_{}^2}=A\left(\frac{1}{p_0+p_{}}\right)^m,$$
(33)
with $`p_0=1.71\mathrm{Ge}\mathrm{V}`$ and $`m=12.42`$.

The $`p_{}`$ distributions of $`D`$ mesons produced in hadron collisions were experimentally found  to be well described by the following simple parameterization:

$$\frac{d\sigma }{dp_{}^2}=C\left(\frac{1}{bm_c^2+p_{}^2}\right)^{\frac{n}{2}},$$
(34)
with $`b=1.4\pm 0.3`$, $`n=10.0\pm 1.2`$, and $`m_c=1.5`$ GeV. This parameterization also provides a good fit to the theoretical calculations  in perturbative QCD. Using this parameterization in (22) gives

$$\nu =\frac{np_{}}{a^2+p_{}^2},a^2bm_c^2.$$
(35)
for the $`\nu `$ parameter in (23).

Fig. 3 shows the ratio of quenching factors (30) for heavy and light quarks in cold nuclear matter ($`L=5\mathrm{fm}`$), which is relevant for large-$`p_{}`$ particle production in $`pA`$ collisions. A small value of the transport coefficient, $`\widehat{q}_{\text{cold}}0.01\mathrm{Ge}\mathrm{V}^3`$, translates into a $``$ 15% enhancement of the heavy-to-light ratio (about 1% for $`L=2\mathrm{fm}`$).

According to (31), we expect a larger quenching ratio for a hot medium. Indeed, as Fig. 4 demonstrates, the $`D/\pi `$ ratio should become significantly enhanced as compared to $`pp`$ collisions: assuming a fixed length $`L=5\mathrm{fm}`$ of the hot medium traversed by the quarks, we find a factor of $`2`$ enhancement at $`p_{}5÷10\mathrm{Ge}\mathrm{V}`$.

We only present the ratio of quenching factors for the following reason. As noticed in , in the region of transverse momenta we are considering, the absolute magnitude of quenching turns out to be extremely sensitive to gluon radiation in the few-hundred-$`\mathrm{Me}\mathrm{V}`$ energy range and therefore cannot be quantitatively predicted without detailed understanding of the spectral properties of the medium.

The heavy-to-light ratio, however, proves to be much less sensitive to the infrared region, since gluon radiation off heavy and light quarks is universal in the $`x0`$ limit, see (16). To illustrate this point, in Fig. 4 we show the ratio of quenching factors calculated with the energy restriction upon gluon energies, $`\omega >500\mathrm{Me}\mathrm{V}`$. We see that in the $`510\mathrm{Ge}\mathrm{V}`$ range of $`p_{}`$, this modifies the ratio by $`2030`$ %.444The quenching factors themselves change (increase) by an order of magnitude when the radiation of gluons with energies smaller than $`500\mathrm{Me}\mathrm{V}`$ is vetoed.

The predicted ratio should not be taken at its face value: the enhancement factor originates mainly from a very strong quenching of light quark jets, and in reality such a dramatic suppression will be washed out by light-particle production at the periphery of the collision zone. In a central $`AA`$ collision, the number of collisions in this peripheral zone, depending on the strength of absorption, scales as $`A^{1/3}÷A^{2/3}`$, whereas the total number of collisions which contribute to the production of $`D`$ mesons scales as $`A^{4/3}`$. Thus, even with account of the collision geometry, we may still expect an order of $`A^{2/3}÷A`$ enhancement of the $`D/\pi `$ ratio, which for a $`AuAu`$ collision may translate in a strong effect.

Moreover, given a large nucleus and moderately large $`p_{}`$, both light and (especially) heavy quarks will be tempted to turn into hadrons prior to leaving the interaction arena. In the pion case, this is likely to add to quenching via pion absorption in the medium , thus enhancing further the $`D/\pi `$ ratio.

Clearly, detailed calculations have to be performed before a defendable number for the magnitude of the $`D/\pi `$ enhancement can be presented. Nevertheless, the $`D/\pi `$ ratio appears to be extremely sensitive to the density of colour charges in QCD matter, and we eagerly await the results of experimental studies of this quantity in relativistic heavy ion collisions.

## Acknowledgments

We are grateful to Rolf Baier, Al Mueller and Dominique Schiff for valuable discussions. D.K. wishes to thank the groups of LPTHE, Université Pierre et Marie Curie and LPT, Université Paris-Sud, for their hospitality during the period when this work was done. The research of D.K. is supported by the U.S. Department of Energy under contract No. DE-AC02-98CH10886.


text:
# HBT and Fluctuations: Recent Results

## 1 Introduction

The study of small relative momentum correlations, a technique also known as HBT  interferometry, is one of the most powerful tools at our disposal to study complicated space-time dynamics of heavy ion collisions . It provides crucial information which helps to improve our understanding of the reaction mechanisms and to constrain theoretical models of the heavy ion collisions. It is also considered to be a promising signature  of the Quark Gluon Plasma (QGP). Interpretation of the extracted HBT parameters in terms of source sizes and lifetime is more or less straightforward for the case of chaotic static sources. In the case of expanding sources with strong space-momentum correlations (due to flow, etc.) the situation is more difficult, but the concept of length of homogeneity  provides a useful framework for the interpretation of data.

## 2 HBT Results

New data from NA45  (CERES Collaboration) and NA49  experiments at CERN shed new light on the evolution of the reaction dynamics at SPS energies. Both collaborations showed results of the systematic studies of the HBT parameters in Pratt-Bertsch parameterization  for Pb+Pb(Au) collisions at 40A GeV. Taking into account preliminary status of both analyses the agreement between them is fairly good. The NA49 Collaboration analysis showed a striking similarity (including detailed $`K_T`$ dependence, where $`K_T`$ is the average transverse momentum of a pair) between the radius parameters at 40 and 160A GeV. This, coupled with results of the transverse radial flow analysis , suggests that the reaction dynamics, despite the difference in the initial energy density and final state multiplicities, are rather similar.

The STAR Collaboration has shown  preliminary results of the extensive analysis of the pion (both $`\pi ^+`$ and $`\pi ^{}`$) correlation functions measured at RHIC in Au+Au collisions at $`\sqrt{S_{NN}}`$=130 GeV. While an increase in radius parameters with charged particle multiplicity was more or less expected, the observed $`K_T`$ dependence, shown in Fig. 1, is somewhat of a surprise. The general trend of the HBT radius parameters with $`K_T`$ is consistent with strong space-momentum correlations due to transverse flow, however, the $`K_T`$ dependence of the ratio of $`R_O`$/$`R_S`$ contradicts the model predictions . Model calculations predict the ratio to be greater than unity due to system lifetime effects which cause $`R_O`$ to be larger than $`R_S`$. They also predict that the $`R_O`$/$`R_S`$ ratio increases with $`K_T`$. Such an increase seems to be a generic feature of the models based on the Bjorken-type , boost-invariant expansion scenario. Hence, it was surprising to see that the experimentally observed ratio is less than unity and is decreasing as a function of $`K_T`$. Currently, it is far from clear what kind of scenario can lead to such a puzzling $`K_T`$ dependence. In principle, the observed behaviour is consistent with the creation of an opaque source  in the collisions. However, to address this particular possibility seriously one needs to perform a detailed analysis of the pion correlation function in Yano-Koonin-Podgoretski parameterization .

Preliminary results of the HBT analysis presented  by the PHENIX Collaboration agree with the STAR results within error bars but do not provide a sufficient cross check due to the limited statistics of the first analysis. Clearly, further cross checks of these results by other RHIC experiments are very desirable.

New results from RHIC and lower SPS energy significantly extended HBT systematics. Figure 2 shows evolution with energy of the Pratt-Bertsch parameters of the two-particle correlation functions measured  in central heavy-ion collisions. It is evident from the figure that HBT parameters change most rapidly at the AGS energies ($`\sqrt{S_{NN}}`$2-5 GeV). This pronounced decrease with energy is driven by an increase in transverse flow. Growing space-momentum correlations lead to a shrinking apparent source size. Saturation in flow velocity  (which is reached around the highest available beam energy at the AGS) leads to a different evolution of the HBT parameters at the SPS energy range ($`\sqrt{S_{NN}}`$10 - 20 GeV). One can see that transverse HBT radii change very little at SPS. The observed behaviour strongly suggests that the HBT correlation functions directly reflect properties of the collective transverse flow generated in the collision.
A consistent picture seems to emerged from the energy dependence systematics of the pion phase space densities reported during the Conference. Preliminary results shown  by the STAR Collaboration seem to confirm the hypothesis of a universal freeze-out phase space density  at RHIC energies. The $`P_T`$ dependence of the phase space density at midrapidity at RHIC clearly rules out a static thermal source and is consistent with a Bose-Einstein distribution modified to take into account radial flow. The flow velocity extracted from the analysis of the $`P_T`$ dependence is quite high ($`0.58c`$) with a freeze-out temperature of about 94 MeV.
Study of the freeze-out phase space density by the NA45 Collaboration  also confirmed a universal phase space density at 40A GeV. They also reported that the phase space density seems to be independent of centrality.
The E895 Collaboration showed  that at AGS energies (2-8A GeV) the phase space density distribution is lower than the “universal” one and gradually approaches it as energy increases.
It was evident during the Conference that studies of the correlations of neutral ($`\mathrm{\Lambda }`$, $`K^0`$) strange particles attracted a lot of attention. Even though the details of the final state interaction of neutral stange particles are not very well known, which complicates interpretation of the measurements, the physics which can be accessed via this type of correlations is clearly attractive. In combination with measurements of strange particles yields and spectra they can provide valuable new information about strangeness dynamics in the collision.
From the experimental point of view such correlations have certain advantages. Being neutral they pose no Coulomb repulsion problems, and two-track resolution problems are largely avoided as well. They are also less influenced by the resonance decays. On the other hand such studies pose significant experimental challenges since they usually require sizable datasets and large acceptance detectors due to the relatively low yields and reconstruction efficiencies of neutral strange particles. Such data were generally unavailable until now in heavy ion collisions.
The NA49 Collaboration presented  preliminary analysis of the $`\mathrm{\Lambda }`$-$`\mathrm{\Lambda }`$ correlations from Pb+Pb collisions at 158A GeV. The shape of the measured correlation function shows a pronounced dip at low relative momentum which is consistent with Fermi-Dirac statistics of $`\mathrm{\Lambda }`$ hyperons. Also there seems to be evidence of weak repulsive $`\mathrm{\Lambda }`$-$`\mathrm{\Lambda }`$ potential which may have some consequences for H-dibaryon searches. This interesting topic certainly requires further investigation.
The STAR Collaboration showed preliminary results of the analysis of $`K^0`$-$`K^0`$ correlations measured in Au+Au collisions at $`\sqrt{S_{NN}}`$=130 GeV. The extracted Gaussian radius $`R_{inv}6fm`$ is consistent with pion HBT measurements.
The E895 Collaboration studied another exotic system - $`\mathrm{\Lambda }`$-p - in Au+Au reactions at lower AGS energies. Preliminary results  are consistent with an attractive strong interaction between proton and $`\mathrm{\Lambda }`$ hyperon but cannot be described by the simple source modes with an Urbana type potential .
The E895 Collaboration also showed  results of the HBT analysis in the impact-parameter-fixed coordinate frame. In this approach event-by-event information about orientation of the reaction plane is combined with information from two-particle correlations. Such analysis is not only sensitive to the length of homogeneity and time, but also to the tilt of the source in the reaction plane. Experimental access to this level of geometrical detail of the freeze-out distribution is unprecedented, and represents a new opportunity to study dynamics of heavy-ion collisions.

## 3 Fluctuations

In recent years analyses of the event-by-event fluctuations assumed one of the central positions in the experimental studies of the relativistic heavy ion collisions. That was evident during the Conference.
The NA45 Collaboration showed  analysis of the mean transverse momentum fluctuations in Pb+Au central collisions at 40, 80 and 158A GeV. They reported observation of small ($`<3\%`$), but statistically significant dynamical fluctuations at all three energies. The observed $`\mathrm{\Phi }_{P_T}`$ value  is $``$8 MeV at 158A GeV. This interesting result contradicts the small $`\mathrm{\Phi }_{P_T}`$ 0.6 MeV/c value reported earlier  by the NA49 Collaboration. Clearly, this discrepancy requires further study.
The NA49 Collaboration reported  preliminary results of the analysis of the event-by-event charge fluctuations in central Pb+Pb collisions at 40 and 158A GeV. They studied fluctuations of ratios of the number of positive and negative hadrons as a function of the width of the pseudorapidity window and applied acceptance corrections following the procedure suggested by Bleicher, Jeon and Koch . Fully corrected results are somewhat surprising: fluctuations of charge particle ratios at both energies are similar and consistent with fluctuations produced in a pion gas. Naively one would expect at SPS energies somewhat smaller fluctuations consistent with a resonance gas or perhaps even with QGP.
The WA98 Collaboration presented results of the search for localized fluctuations in the multiplicity of charged particles and photons produced in Pb+Pb reactions at 158A GeV. They reported individual fluctuations in photons and charged particle multiplicities but no correlated event-by-event fluctuations of $`N_\gamma N_{ch}`$.
The NA44 Collaboration reported  preliminary results of their search for local fluctuations in the number of charged particles using a novel event texture analysis. No evidence of critical fluctuations was found.
For RHIC, the STAR Collaboration reported  preliminary results of their extensive event-by-event fluctuations analysis. They presented an investigation of the charge-dependent and charge-independent $`<P_T>`$ fluctuations individually, as well as the overall $`<P_T>`$ fluctuations. A significant excess in fluctuations of $`\mathrm{\Delta }\sigma _{P_T}35`$ MeV/c beyond the statistical expectation was observed in the charge-independent analysis. Also, a somewhat smaller but still statistically significant excess was observed in the charge-dependent analysis. They also studied the centrality dependence of these excess fluctuations where hints of a decrease with centrality in the charge-dependent part was observed. The results of multiplicity fluctuations analysis show no evidence for predicted fluctuation suppression associated with rapid hadronization of the QGP.
Event-by-event fluctuations also drew a lot of attention from theorists. Most of the theoretical studies reported at this Conference were devoted to the fluctuations of so-called “conserved quantities”: charge, baryon number, strangeness, etc.
Koch discussed the theoretical situation with conserved quantities fluctuations in general and charge ratios in particular. He pointed out effects of limited detector acceptance and the importance of related corrections. It was apparent during the Conference that current theoretical interest is aimed at the studies of the “survivability” of the proposed observables in the real world situation. Most of the results were devoted to the studies of the influence of the system expansion and limited experimental acceptance. These were discussed in talks by Asakawa , Stephanov , Gavin .
Pratt  presented a new observable for determination of the nature of hadronization in heavy ion collision - the balance functions. The balance function describes the conditional probability that a particle in the momentum bin $`p_1`$ will be accompanied by a particle of opposite charge in the bin $`p_2`$ . Balance functions probe hadronization by quantifying charge correlations in momentum space and potentially allow us to distinguish scenarios with QGP formation from purely hadronic ones. First experimental attempts to apply this observable to heavy ion collisions are underway.

## 4 Conclusions

Many new and interesting HBT results have appeared since last Quark Matter. Combined results from the AGS, SPS and RHIC allow us to have a systematical view of the evolution of the HBT parameters with beam energy. Presentations of the first results from RHIC experiments were some of the highlights of the conference. The most intriguing feature of the preliminary HBT results from RHIC is $`K_T`$ dependence of the $`R_O/R_S`$ ratio reported by the STAR Collaboration. It is not clear yet what physics can lead to such a dependence.
It would be interesting to perform measurements between AGS and SPS energy domains, which seems to be a transition region in terms of HBT systematics. It would be also useful to probe the energy range between SPS and RHIC. Such measurements will significantly extend our understanding of spatio-temporal evolution of heavy ion collisions.
The situation with fluctuation measurements at SPS is far from clear and further experimental and theoretical efforts are needed to resolve existing ambiguities.
Preliminary fluctuations results from RHIC are tantalizing and warrant further investigation.


text:
# Halo Substructure and Disk Heating in a ΛCDM universe

## 1. Introduction

Over the past few years cosmological N-body simulations have achieved unprecedented resolution. One recent highlight of these studies has been the discovery that during the merger events that characterize the assembly and growth of dark matter halos the central regions of accreted halos may survive for several orbital times as self-bound entities in the parent halo (Klypin et al 1999, Moore et al 1999, Ghigna et al 2000).

This population of surviving halo cores, or “subhalos”, typically contributes less than $`10\%`$ of the total mass of the system, with the bulk of the mass in a smooth monolithic structure, as envisioned by White & Rees (1978). Despite the small fraction of the total mass they make up, at any given time a large number of subhalos are expected within the virialized region of a cold dark matter halo. In a galaxy cluster, it is tempting to identify these subhalos as dynamical tracers of the individual galaxy population (Moore et al 1999, hereafter M99), but in the case of a galaxy-sized halo such association would lead to the expectation of several hundred satellites in the vicinity of the Milky Way. For example, M99 find that up to $`500`$ satellites with circular velocities exceeding $`10`$ km/s may have survived within $`300`$ kpc from the center of the Galaxy. Comparing this large number with the dozen or so known Milky Way satellites with comparable velocity dispersion clearly implies that most subhalos must have failed to form a significant number of stars if these models are to match observations.

The lack of a significant luminous component associated with most subhalos was suggested several years ago by semianalytic models of galaxy formation (see, e.g., Kauffmann, White & Guiderdoni 1993), and is analogous to the need to prevent a large number of dwarf galaxies from forming in low-mass halos first noted by White & Rees (1978) and confirmed later by semianalytic models of galaxy formation (White & Frenk 1991, Cole et al 1994, Somerville & Primack 1999). Energetic mechanisms that operate more efficiently on low-mass systems, such as feedback from evolving stars (White & Rees 1978), as well as heating by an ionizing UV background (Efstathiou 1992, Bullock, Kravtsov & Weinberg 2000), are usually invoked to explain this decoupling of luminous and dark components on small mass scales and are a staple of hierarchical galaxy formation models.

Even if luminous galaxies fail to “turn on” in most subhalos a potential difficulty with the substructure expected in CDM halos has been cited by M99: the fluctuating gravitational potential induced by the clumpy structure of the halo may act to heat and thicken fragile stellar disks beyond observational constraints (see, e.g., Toth & Ostriker 1992). M99 estimate that the energy pumped into a disk like that of the Milky Way by subhalos may add up to a significant fraction of the binding energy of the disk over a Hubble time. Estimating analytically the heating rate due to the cumulative effect of disk/subhalo collisions is, however, quite uncertain, and in practice simple calculations that neglect the self-gravitating response of the disk often overestimate the effects of heating by infalling satellites (Huang & Carlberg 1997, Velázquez & White 1999).

In this Letter we present the results of an attempt to gauge the dynamical effects of substructure on the dynamical evolution of thin stellar disks embedded in dark matter halos. Our approach uses the mass spectrum, orbit distribution, and internal structure of subhalos identified in cosmological simulations of galaxy-sized CDM halos.

## 2. Substructure in CDM halos

The mass function of substructure halos, their internal structure, and the parameters of their orbits are the main properties of the subhalo population that determine the tidal effects of substructure on stellar disks. Subhalos discussed here and in M99 have been identified using the groupfinder SKID111http://www-hpcc.astro.washington.edu/tools/skid.html. SKID’s main input parameter is a characteristic length scale, $`l_S`$, which determines the characteristic size of the smallest groups and, indirectly, the minimum separation between clumps.

Fig. 1. Cumulative circular velocity function of substructure halos. Velocities are normalized to the circular velocity of the parent halo measured at the virial radius. Only subhalos within the virial radius of the main halo are included. Short-dashed curves correspond to two galaxy-sized halos formed in a SCDM universe, derived from data kindly provided by M99. Long-dashed curves correspond to a halo formed in the $`\mathrm{\Lambda }`$CDM cosmogony. The two sets of curves, labeled $`V_{\mathrm{peak}}`$ and $`V_{\mathrm{outer}}`$, indicate the result of adopting the maximum circular speed or the circular velocity at the outermost bound radius, respectively. The dot-dashed curve illustrates the results reported by Klypin et al (1999). Filled squares correspond to Milky Way satellites assuming that stars in dwarf spheroidals are on isotropic orbits in isothermal potentials. The dotted line corresponds to the same data but where circular velocities have been increased by a factor of three to reflect uncertainties in the structure of dark halos surrounding dwarf spheroidals (see, e.g., White 2000).

The short-dashed curves in Figure 1 show the subhalo velocity function corresponding to two galaxy-sized dark matter halos formed in the “Local Group” run reported by M99. The “Local Group” simulation assumed the former “standard” CDM cosmogony (SCDM: $`\mathrm{\Omega }_M=1,\mathrm{\Omega }_\mathrm{\Lambda }=0,h=0.5,\sigma _8=0.7`$). The two halos shown here have $`V_{200}193`$ and $`172`$ km s-1, respectively, and each has roughly $`10^6`$ particles within $`r_{200}`$, the “virial” radius within which the mean density is $`200`$ times the critical density for closure.

Circular velocity (instead of mass) is used to characterize subhalos because of its weaker dependence on the exact way in which substructure is identified. Still, subhalo circular velocities change with radius, and there is no unique way of defining subhalo circular velocities. SKID reports $`V_{\mathrm{peak}}`$ and $`V_{\mathrm{outer}}`$, which correspond to the maximum circular speed within the subhalo, and to its value at the outermost bound radius, respectively. These two measures may differ in some cases by up to $`50\%`$, implying rather different substructure velocity functions, as can be seen in Figure 1. The substructure function presented in M99’s Figure 2 is consistent with the set of curves labeled $`V_{\mathrm{peak}}`$ in Figure 1222Note that the caption to M99’s Figure 2 suggests, apparently in error, that $`V_{\mathrm{outer}}`$ is used there rather than $`V_{\mathrm{peak}}`$..

Fig. 2. Mass vs pericentric radii for subhalos identified in two CDM halos, after scaling to the virial radius and circular velocity of our Milky Way disk model.

As pointed out by M99 and Klypin et al (1999), the substructure velocity function, scaled to the circular speed of the parent halo, is roughly independent of the mass of the parent halo. The long dashed curves in Figure 1 indicate that this function is also rather insensitive to the cosmological model adopted. These curves correspond to a $`205`$ km s-1 halo formed in the $`\mathrm{\Lambda }`$CDM cosmogony ($`\mathrm{\Omega }_M=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,h=0.65,\sigma _8=0.9`$). This system was selected at random from a large simulation of a $`32.5h^1`$ Mpc periodic box, and resimulated at high resolution using PKDGRAV. The procedure for resimulation is analogous to that described in Eke, Navarro & Steinmetz (2001), and is similar (but not identical) to that used for the SCDM “Local Group” halos. The mass and force resolution are similar in the SCDM and $`\mathrm{\Lambda }`$CDM runs, and they have both been carried out with the same N-body code (PKDGRAV).

The insensitivity to cosmogony of the substructure velocity function was first reported by Klypin et al (1999). Their results are shown by the dot-dashed line in Figure 1 and are in reasonable agreement with our determination, provided that the circular velocities computed by their substructure-finding algorithm fall, as it seems plausible, between $`V_{\mathrm{peak}}`$ and $`V_{\mathrm{outer}}`$.

How does the subhalo velocity function compare with the observed number of satellites in the vicinity of the Milky Way? The solid squares in Figure 1 (joined by a solid line) illustrate the cumulative number of known Milky Way satellites as a function of the circular velocity of their halos, as plotted by M99. Here, circular velocities are derived for the halos of dwarf spheroidals assuming that stars in these systems are on isotropic orbits in isothermal potentials. This is a plausible, but nevertheless questionable, assumption. Dark halos differ significantly from simple isothermal potentials, and numerical simulations indicate that circular velocities decrease substantially near the center. If stars populate the innermost regions of subhalos their velocity dispersions may substantially underestimate the subhalo peak circular velocities.

This has been discussed by White (2000), who finds, using the mass model proposed by Navarro, Frenk & White (1996, 1997, NFW), that dwarf spheroidals may plausibly inhabit potential wells with circular velocities up to a factor of $`3`$ times larger than inferred under the isothermal assumption. Such correction may reconcile, at the high mass end, the Milky Way satellite velocity function with the subhalo function, as shown by the dotted line in Figure 1. Thus, it is possible that the number of massive satellites expected in the CDM scenario might not be in gross conflict with observation.

On the other hand, this does not erase the large number of low mass “dark” satellites that should inhabit the main halo. Can thin stellar disks survive unscathed in this clumpy environment? The tidal heating rate by substructure can be shown to scale as $`dE/dtn(m_s)m_s^2𝑑m_s`$, where $`m_s`$ is the subhalo mass (White 2000). Since, according to numerical simulations, $`n(m_s)m_s^{1.8}`$ (Ghigna et al 2000, Springel et al 2000), the tidal effects of substructure are dominated by the few most massive subhalos. Stochastic effects are therefore expected to dominate the heating rate and a statistically significant sample of systems should be simulated in order to obtain conclusive results. The results reported here (which are based on two realizations) should then be viewed as preliminary.

## 3. The Evolution of Stellar Disks in a $`\mathrm{\Lambda }`$CDM halo

### 3.1.  The Galaxy Model

We have carried out a series of numerical simulations that follow the evolution of a disk of particles within a dark matter halo with substructure similar to that of the high-resolution $`\mathrm{\Lambda }`$CDM halo alluded to in Figure 1. Since we are interested in isolating the effects of substructure on the disk, we choose to minimize the influence of other effects (such as the triaxiality of the halo, spiral structure, disk warping) by carefully constructing an equilibrium disk/bulge/halo model following the procedure outlined by Hernquist (1993). We choose parameters so that the model reproduces many of the observational characteristics of the Milky Way galaxy, as in Velázquez & White (1999). The disk is assumed to be exponential, with a total mass of $`5.6\times 10^{10}M_{}`$, exponential scale length of $`3.5`$ kpc, and vertical scale height of $`700`$ pc. The bulge is modeled as a Hernquist (1993) model with total mass $`1.87\times 10^{10}M_{}`$ and $`525`$ pc scale radius. Finally, we adopt a non-singular isothermal model for the dark halo, with a total mass of $`3.09\times 10^{12}M_{}`$, a core radius $`r_c=3.5`$ kpc, and exponentially truncated at a distance of $`300`$ kpc from the center. The scaled circular velocity profile of the galaxy model closely resembles that of the $`\mathrm{\Lambda }`$CDM halo out to $`600`$ kpc, about twice its virial radius.

The galaxy model is first evolved without substructure halos in an attempt to assess the deviations from equilibrium induced by noise in the particle distribution. N-body disks are notoriously unstable, and substantial numbers of particles are needed to ensure stability over several rotation periods. The simulations we report here use $`40,000`$ particles in the disk, $`13,333`$ in the bulge, and $`2.2\times 10^6`$ in the halo, all of equal mass, $`m_p=1.4\times 10^6M_{}`$. All simulations were run with PKDGRAV using a spline softening length of $`175`$ pc and a hierarchy of timesteps, the smallest of which is typically less than $`10^5`$ yrs.

### 3.2. Including Substructure

Substructure halos, despite their large numbers, make up a small fraction of the total mass of the system, implying that it is possible to add subhalos to our galaxy model without seriously compromising the global equilibrium of the system. We include subhalos using the following procedure. We first run SKID on the CDM halos discussed in §2, using the characteristic linking-length scales that maximize the number of subhalos found within the virial radius. This yields $`685`$ and $`737`$ satellites with masses exceeding $`10^8M_{}`$ within $`2r_{200}`$ for the $`\mathrm{\Lambda }`$CDM and SCDM halo, respectively.

We choose the peak circular velocity, $`V_{\mathrm{peak}}`$, and the radius at which it occurs, $`r_{\mathrm{peak}}`$, to characterize each subhalo. The values of $`r_{\mathrm{peak}}`$ and $`V_{\mathrm{peak}}`$ are then used to generate, for each of the satellites, an equilibrium NFW model with particle masses equal to those used in the disk/halo galaxy model333Note that these two parameters fully specify an NFW halo model.. In order to ensure a finite mass we truncate each NFW model at $`5r_{\mathrm{peak}}`$. The positions and velocities of each satellite are scaled so that they correspond to the same fractions of $`r_{200}`$ and $`V_{200}`$, respectively, in both the parent CDM halos and in the galaxy model. This procedure ensures that the subhalos have approximately the same masses, densities and orbits as those in the CDM halos. Figure 2 shows the masses and pericentric distances of the subhalos within $`2r_{200}`$ included in the simulations. We include satellites within twice the virial radius because many of them have time to reach the pericenter of their orbits (which lie within $`r_{200}`$) during the timespan of our simulations. Note that there are no satellites with masses exceeding $`10^{11}M_{}`$ and that only a few satellites are on orbits that take them within the solar circle ($`8.5`$ kpc). All of those satellites have masses below $`10^9M_{}`$. Figure 3 shows an edge-on view of the disk/bulge/halo/subhalo model.

### 3.3. Disk Heating

We evolve the disk model first for about $`3.5`$ Gyrs without including substructure in order to quantify the heating rate due to noise in the particle distribution. The ($`R`$, $`z`$, $`\varphi `$) velocity dispersions of disk particles at the solar circle (i.e., $`8.5`$ kpc from the center) grow from ($`31`$, $`27`$, $`26`$) km s-1 to ($`43`$, $`30`$, $`33`$) km s-1 over the same period (see solid lines in Figure 4). Quantifying this heating rate by the usual expression $`\sigma _{tot}^2=\sigma _0^2+Dt`$, we find $`D200`$ km2 s-2 Gyr-1, about a factor of two less than inferred for stars in the solar neighborhood from the age-velocity dispersion (Wielen 1977, Lacey & Ostriker 1985, Fuchs et al 2000).

The two-body heating rate in our equilibrium stellar disk model is thus low compared with the actual heating experienced by stars in the disk of the Galaxy. This implies that our disk model is stable enough to verify numerically whether substructure in the halo leads to heating rates inconsistent with observational constraints. This is shown by the dashed lines in Figure 4, which show the evolution of the disk velocity dispersion when the substructure halos are added to the system. Clearly, the heating rate is approximately the same with and without substructure, a result that may be traced to the fact that there are no satellites more massive than $`10^{11}M_{}`$ and that their orbits seldom take them near the disk. If substructure in the halos we consider is representative of galaxy-sized CDM halos (and we have no reason to suspect it is not), this would imply that tidal heating rates of thin stellar disks by substructure halos may be consistent with the observational evidence.

Fig. 4. Stellar disk velocity dispersion as a function of time in our Milky Way model compared with the age-velocity dispersion relation in the solar neighborhood, as compiled recently by Fuchs et al (2000, dot-dashed curve). Solid (dashed) lines indicate the heating rate of our disk/bulge/halo model evolved for about $`3.5`$ Gyrs excluding (including) substructure. “Starting” times for simulation data are arbitrary and have been chosen to lie on the age-$`\sigma _{tot}`$ relation for the solar neighborhood.

## 4. Discussion and Conclusions

Our results suggest that concerns regarding excessive tidal heating of thin stellar disks by substructure in CDM halos may be less serious than previously thought. Once subhalos with the mass, structure, and orbital distribution expected in a CDM universe at $`z0`$ are considered, direct numerical simulation shows that thin stellar disks at the center of such halos may survive undisturbed by collisions with subhalos. This is because subhalo masses in galaxy-sized halos are typically below $`10^{10}`$-$`10^{11}M_{}`$ and have orbits that usually do not take them close to the disk. We conclude that the substructure observed in virialized CDM halos is not clearly inconsistent with the existence of thin stellar disks such as that of the Milky Way.

These conclusions are subject to a number of caveats. The most obvious one is that our study explores only two numerical realizations of a disk galaxy within clumpy dark matter halos, and it is always hazardous to extrapolate from such a small number of cases. Our study does show, however, that it is possible at least in some cases to maintain a stellar disk in spite of substructure. A second caveat is that we have explored a model motivated by the present-day structure of the Milky Way and by the $`z=0`$ substructure of a CDM halo. Models that take into account the ongoing formation of the disk and a more realistic treatment of the evolution of substructure are clearly desirable in order to refine the conclusions presented here. Finally, for simplicity we have adopted a non-singular isothermal model for the dark halo in the simulations reported here, rather than the ‘cuspy’ density profiles found in cosmological N-body simulations (NFW). Although we do not expect this assumption to undermine the conclusions reported here, we cannot rule out subtle effects that may arise from the interaction between disk and halo in a cosmological context.

Finally, as discussed by Navarro & Steinmetz (1997, 2000), it is quite difficult to account simultaneously for the masses, luminosities, rotation speeds, and angular momenta of galaxy disks in cosmogonies such as CDM, where much of the mass of a virialized halo tends to be assembled through a sequence of mergers. Until these issues are fully resolved it would be premature to extend a clean bill of health to the CDM paradigm regarding the formation and evolution of spiral galaxies like our own Milky Way.

We would like to thank B. Moore, F. Governato, S. Ghigna and G. Lake, for sharing with us an electronic version of the data showed in Figure 1. Lars Hernquist kindly made available the software used to set up the equilibrium disk/bulge/halo models. We thank the referee for a useful report that led to a number of clarifications and improvements. The Natural Sciences & Engineering Research Council of Canada (NSERC) and the Canadian Foundation for Innovation have supported this research through various grants to JFN.