text:
# Ionization history of the cosmic plasma in the light of the recent CBI and future PLANCK data.

## 1 Introduction

Recent observations of the CMB anisotropy such as the BOOMERANG (de Bernardis et al. 2000), the MAXIMA-1 (Hanany et al. 2000), the DASI (Halverson et al. 2001), the VSA (Watson et al. 2002), the Archeops (Benoit et al. 2002) and especially the new CBI data (Mason et al. 2002) and the DASI polarization measurements (Kovac et al. 2002) provide the good base for the progress in understanding the most general properties of the early Universe. In these experiments, major attention is on the statistical properties of the signal and the determination of the power spectrum of the CMB anisotropy and polarization. The obtained cosmological parameters for adiabatic cold dark matter model (CDM) are well consistent with the Big Bang nucleosynthesis predictions, SN Ia observational data and results obtained from investigations of galaxy surveys. Further progress will be achieved with more sensitive experiments such as the launched MAP and the upcoming PLANCK missions. They will be able to measure the CMB anisotropy and polarization power spectra up to the multipole range $`l3410^3`$ where the lensing effect, Sunyaev-Zel’dovich effect in clusters, secondary anisotropies and some peculiarities of the ionization history of the Universe dominate.

The high precision of the future CMB measurements allows to reveal some distortions with respect to the standard model of the hydrogen recombination process occurred at redshifts $`z10^3`$. Here we show that polarization power spectrum contains important information about the kinetics of hydrogen recombination and allows to trace the ionization history of the cosmic plasma in wide range of redshifts $`10^3z610`$.

For the baryon–dominated Universe the classical theory of the hydrogen recombination was developed in Peebles (1968) and Zel’dovich, Kurt and Sunyaev (1968). For the dark matter dominated Universe, it was generalized in Zabotin and Naselsky (1982), Jones and Wyse (1985), Seager, Sasselov and Scott (1999) (see detailed review in White, Scott and Silk, 1994). Many distortions of the standard model of recombination have also been investigated. The delay of recombination due to evaporation of primordial black holes has been discussed in Naselsky (1978) and Naselsky and Polnarev (1987). Avelino et al. (2000), Battye et al. (2001), Landau et al. (2001) have noted that possible time variations of the fundamental constants could also influence the ionization history of the cosmic plasma. Sarkar and Cooper (1983), Scott et al. (1991), Ellis et al. (1992), Adams et al. (1998), Peebles, Seager and Hu (2000), Doroshkevich and Naselsky (2002) discussed distortions of the kinetics of hydrogen recombination caused by decays of hypothetical particles. It is worth noting that for such models the energy injection delays the recombination at $`z10^3`$ and distorts the ionization history of the Universe down to the period of galaxies formation at $`z510`$.

Recently Naselsky and Novikov (2002) have proposed a model of accelerated recombination with non–uniform spatial distribution of baryonic matter. In this model the recombination proceeds faster within the clumps and slowly in the intercloud medium. In such models at redshifts $`z10^3`$ the mean ionization fraction decreases faster than that in the standard model with the same cosmological parameters which mimics the acceleration of recombination. Here we show that, for suitable parameters of the clouds, such models can improve the best–fit of the observed power spectrum of the CMB anisotropy.

The non-standard models of hydrogen recombination are characterized by a few additional parameters which must be included in the list of the best-fit cosmological parameter determination using the CMB anisotropy and polarization data. We call these parameters as the “missing parameters” of the theory and show how important these parameters may be for the CMB physics.

In this paper we compare the observed CMB anisotropy power spectrum for the reference model of the standard recombination, with models of delayed and accelerated recombination. For the reference model we accept the cosmological parameters from the best-fit of the CBI observational data (Mason et al 2002). We show that the models considered in this paper can provide better fits for the already observed and the future MAP and PLANCK power spectra and can change usual estimates of the cosmological parameters. Moreover, the measurements of the polarization allow us to discriminate between some of such complex models of the Universe. These results can be important for the interpretation of the MAP and upcoming PLANCK data.

The paper is organized as follows. In section 2 we discuss the generic models of distorted recombination. In section 3 we demonstrate distortions of the CMB anisotropy power spectrum due to delay or acceleration of recombination. In section 4 we introduce the models of accelerated and delayed recombination and discuss the redshift variations of the hydrogen ionization fraction in these models. In section 5 we compare the observed CMB anisotropies and expected polarization power spectra with ones for the models under discussion. Section 6 is devoted to discussion of late reionization and of the polarization power spectrum for the delayed recombination models. In Conclusion we make some predictions for the MAP and the upcoming PLANCK experiments.

## 2 Deviations from the standard model of hydrogen recombination

The deviations from the standard recombination process caused by the injection of additional Ly-$`\alpha `$ and Ly-c photons at the recombination epoch or by strong small scale clustering of the baryonic component can be suitably described in terms of an additional source of Ly-$`\alpha `$ and Ly-c photons which was proposed in Peebles, Seager and Hu (2000). In the generalized recombination model the rates of production of resonance, $`n_\alpha `$, and ionized, $`n_i`$, photons are described by functions $`\epsilon _\alpha (z)`$ and $`\epsilon _i(z)`$, defined as follows:

$$\frac{dn_\alpha }{dt}=\epsilon _\alpha (z)H(z)n_b(z),$$
(1)
$$\frac{dn_i}{dt}=\epsilon _i(z)H(z)n_b(z),$$
(2)
where $`H(z)`$ and $`n_b(z)`$ are the Hubble parameter and the mean baryonic density.

### 2.1 Delay of recombination

Peebles, Seager and Hu (2000) considered the models with $`\epsilon _\alpha (z)=const.`$ and $`\epsilon _i(z)=const.`$ However, in general case both parameters $`\epsilon _\alpha (z)`$ and $`\epsilon _i(z)`$ could depend on redshift (see, e.g., Doroshkevich & Naselsky 2002; Naselsky & Novikov 2002). As was shown in these papers, for the generation of ionized photons from the decay of super heavy dark matter (SHDM) particles we get

$$\epsilon _\alpha (z)\epsilon _i(z)(1+z)^1.$$
(3)
For models of particle or strings decays the particle number density, $`n_x`$, is a decreasing function of time defined by the life time of the particles, $`\tau _x(n_x,t)`$, as follows:

$$\frac{dn_x}{dt}+3H(t)n_x=\frac{n_x}{\tau _x(n_x,t)}$$
(4)
For the simplest models with $`\tau _x=const.`$ we have

$$\frac{d[n_x(1+z)^3]}{dt}=\epsilon _x(t)H(t)n_b,$$
(5)
$$\epsilon _x=\frac{1}{H(t)\tau _x}\mathrm{exp}\left(\frac{tt_u}{\tau _x}\right)\frac{n_x(t_u)}{n_b(t)}$$

where $`t_u`$ is the age of the Universe.

Obviously, $`\epsilon _i=\kappa _i\epsilon _x,\epsilon _\alpha =\kappa _\alpha \epsilon _x`$ where factors $`\kappa _i`$ and $`\kappa _\alpha `$ characterize the efficiency of resonance and ionized photons production. The decays of $`x`$-particles cannot significantly distort the thermodynamics of the Universe at high redshifts $`z10^3`$ but this energy injection changes the kinetic of recombination at $`z10^3`$ that leads to the observable distortions of the power spectra of CMB anisotropy and polarization. To evaluate these distortions we firstly have to follow the transformation of high energy injected particles to UV photons influenced directly on the recombination process.

The electromagnetic cascades initiated by the ultra–high energy cosmic rays (UHECR) are composed of photons, protons, electrons, positrons and neutrinos. At high redshifts, the cascades develop very rapidly via interaction with the CMB photons and pair creation ($`\gamma _{UHECR}+\gamma _{CMB}e^++e^{}`$), proton-photon pair production ($`p_{UHECR}+\gamma _{CMB}p^{^{}}+\gamma {}_{}{}^{}+e^++e^{}`$), inverse Compton scattering ($`e_{UHECR}^{}+\gamma _{CMB}e^{^{}}+\gamma ^{^{}}`$), pair creation ($`e_{UHECR}^{}+\gamma _{CMB}e^{^{}}+e^{}+e^++\gamma ^{^{}}`$), and, for neutrino, electron- positron pair creation through the Z-resonance. As was shown by Berezinsky et al. (1990) and Protheroe et al. (1995), these processes result in the universal normalized spectrum of a cascade with a primary energy $`E_\gamma `$ which can be written as follows:

$$N_\gamma (E,E_\gamma )=F(E,E_\gamma )\{\begin{array}{cc}\sqrt{\frac{E}{E_a}}& EE_a\\ 1& E_aEE_c\\ 0& E_cE\end{array}$$
(6)
$$F(E,E_\gamma )=\frac{E_\gamma E^2}{2+\mathrm{ln}(E_c/E_a)},_0^{E_\gamma }EN_\gamma 𝑑E=E_\gamma $$

where $`E_c4.610^4(1+z)^1`$GeV, $`E_a=1.810^3(1+z)^1`$GeV. At the period of recombination $`z10^3`$ and for lower redshifts both energies, $`E_a`$ and $`E_c`$, are larger than the limit of the electron-positron pair production $`E_{e^+,e^{}}=2m_e=1`$ MeV and the spectrum (6) describes both the energy distribution at $`EE_{e^+,e^{}}`$ and the injection of UV photons with $`EE_{e^+,e^{}}`$. However, the spectrum of these UV photons is distorted due to the interaction of photons with the hydrogen - helium plasma.

Let us briefly discuss the relationship between the initial spectra of the SHDM decay and the parameters $`\epsilon _\alpha `$ and $`\epsilon _i`$ of hydrogen excitation and ionization. Evidently, the energy density of UV photons is proportional to the rate of energy injection, $`\dot{Q}(t)`$,

$$N(E)\dot{Q}(t)(E_a/E)^{3/2}=\omega (t)\dot{Q}(t_u)(E_a/E)^{3/2},$$
(7)
$$\dot{Q}(t)=\omega (t)\dot{Q}(t_u)=_0^{\mathrm{}}𝑑EEN(E),$$

where $`\dot{Q}(t_u)`$ is the rate of energy injection at $`t=t_u`$. For the net of ionizing UV photons with $`EI`$ 13.6 eV we have

$$\dot{n}_i\frac{\omega (t)\dot{Q}(t_u)}{I}\left(\frac{I}{E_a(z)}\right)^{\frac{1}{2}}$$
(8)
Using normalization $`\dot{Q}(t_u)H^1(t_u)ϵ_{EGRET}`$ to the EGRET energy density $`ϵ_{EGRET}410^7eVcm^3`$ (Berton, Sigl & Silk 2001) we can estimate the functions $`\epsilon _i`$ and $`\epsilon _\alpha `$. Namely, for

$$\omega (t)t^{4p}(1+z)^{\frac{3}{2}(4p)}$$
(9)
we get

$$\epsilon _i\epsilon _\alpha \xi (z)\left(\frac{ϵ_{EGRET}}{In_b(z=0)}\right)(1+z)^{\frac{3}{2}(1p)}(1+z)^{2(13p/4)}$$
(10)
$$\xi (z)=\left(\frac{I}{E_a(z)}\right)^{\frac{1}{2}}310^6\sqrt{1+z}$$

The power index $`p=1`$ corresponds to the model involving the release of $`x`$-particles from topological defects, such as ordinary cosmic strings, necklaces and magnetic monopoles (see, e.g., Sigl 2001). The model with $`p=2`$ corresponds to the exponential decay of SHDM particles described by Eq. (3) with $`\tau _xt_U`$. The model discussed in Peebles, Seager and Hu (2000) with $`\epsilon _i=const.`$ and $`\epsilon _\alpha (z)=const.`$ corresponds to $`p=4/3`$.

For some models such as those of evaporation of primordial black holes or decays of SHDM relics with shorter life time, $`\tau _X=t(z_x)t_u`$, more complex normalization is required. Thus, for $`\tau _X0.1t_u`$ that corresponds to redshift $`z_x67`$, these decays occur before or just at the epoch of galaxy formation and now their number is small. For such models and for $`zz_x`$ we get, instead of Eq. (10),

$$\xi _{mod}(z)=\xi (z)\mathrm{\Theta }(z_x)=\xi (z)\mathrm{exp}\left((1+z_x)^{\frac{3}{2}}\right).$$
(11)
For some models with different $`p`$ and $`\mathrm{\Theta }(z_x)`$, functions $`\epsilon _\alpha (z)`$ are plotted in Fig. 1. As one can see, for some of the models we should have $`\epsilon _\alpha 10^310^2`$ at $`z<1500`$ and distortions of the kinetics of hydrogen recombination can be quite strong. For other models these distortions are small and cannot be observed.

Other approach to restrictions of admissible energy injection is discussed in Sec. 6 .

### 2.2 Acceleration of recombination

The basic idea of accelerated recombination is very transparent (Naselsky and Novikov 2002). Let us consider the model in which a fraction of baryons is concentrated in high density low mass clouds. The hydrogen recombination inside the clouds proceeds faster than in the low density intercloud medium, and the mean ionization fraction is

$$x_e(z)=x_{e,in}Z_m+x_{e,out}(1Z_m)$$
(12)
where $`Z_m`$ is the mass fraction accumulated by clouds, $`x_{e,in}`$ and $`x_{e,out}`$ are the ionization fractions within and between clouds. This expression clearly shows that in such models the redshift variations of $`x_e(z)`$ differ from those in the standard models with homogeneous spatial distribution of baryons and at the same values of cosmological parameters.

Cosmological models with such clumpy distribution of baryons appear when we consider the mixture of adiabatic and isocurvature (or isotemperature) perturbations (Dolgov & Silk 1993). In suitable models isocurvature perturbations provide the concentration of some fraction of baryons within low massive clouds with $`M10^210^5M_{}`$ while adiabatic perturbations dominate on larger scales and are responsible for galaxy formation and other effects.

The ionization history of hydrogen in such models can also be described by Eqs. (1, 2) with $`\epsilon _\alpha 0`$ and $`\epsilon _i0`$ (Naselsky and Novikov 2002). Of course, this choice does not mean that Ly-$`\alpha `$ photons disappear or are really transformed to Ly-c ones. But the suitable choice $`\epsilon _\alpha 0`$ allows to describe correctly the faster hydrogen recombination within the clouds while the choice $`\epsilon _i0`$ allows for the excess of Ly-$`\alpha `$ photons and the delay of hydrogen recombination in the low density intercloud medium.

### 2.3 Three scenarios of recombination distortions

For wide set of models, the functions $`\epsilon _\alpha `$ and $`\epsilon _i`$ can be written as follows:

$$\epsilon _{i,\alpha }(z)=\underset{j}{}\mathrm{exp}\left[\zeta ^{m_j}\right]c_j^{i,\alpha }\zeta ^{n_j},\zeta =(1+z)/(1+z_d),$$
(13)
where parameters $`c_j^{i,\alpha },z_d,m_j,`$ and $`n_j`$ characterize the $`j^{th}`$ source of distortion. It is natural to expect that $`\epsilon _i(z)`$ 0 for all redshifts $`z10^3`$. This additional ionization suppresses partly the CMB anisotropy because of the growth of the optical depth for the Compton scattering at $`z10^3`$.

Three scenarios can be discussed depending upon the shape of the function $`\epsilon _\alpha (z)`$. Namely, the models with $`\epsilon _\alpha (z)`$ 0 and $`\epsilon _\alpha (z)`$ 0 describe the delay or acceleration of recombination, respectively. More complex models with the acceleration of the recombination at $`zz_{ad}`$ and the delay of recombination at $`zz_{ad}`$ are described by complex functions with $`\epsilon _\alpha (z)`$ 0 at $`zz_{ad}`$ and $`\epsilon _\alpha (z)`$ 0 for $`zz_{ad}`$.

All these scenarios can be based on the realistic physical models. For example, for the decay of the long lived SHDM particles we get

$$\epsilon _\alpha \epsilon _i(1+z)^1,$$
(14)
which corresponds to the first scenario. The same scenario is realized for the decay of topological defects and evaporation of primordial black holes discussed above. In all the cases, the Ly-c and Ly-$`\alpha `$ photons are formed via the electromagnetic cascades and $`\epsilon _\alpha \epsilon _i`$. Third scenario is realized for the clumpy baryonic model discussed above. The list of the physical models can be essentially extended.

## 3 Distortion of the CMB anisotropy power spectrum

As was mentioned in Peebles, Seager and Hu (2001), the distortions of recombination shift the redshift of the last scattering surface, $`z_r`$, and change its thickness, $`\mathrm{\Delta }_z`$. In turn, these variations shift positions of the peaks in the CMB anisotropy and polarization power spectrum (Hu and Sugiyama 1995) and change their amplitudes. Here we give the rough analytical estimates of these influences. More accurate numerical results are presented below.

Both parameters of the last scattering surface, $`z_r`$ and $`\mathrm{\Delta }_z`$, are roughly expressed through the Thomson optical depth, $`\tau _T`$, and its derivations, $`\tau _T^{^{}}=d\tau _T/dz`$ :

$$\tau _T=_0^{z_r}\tau _T^{^{}}(z)𝑑z,\tau _T^{^{}}(z)=\frac{c\sigma _Tx_e(z)n_b(z)}{H(z)(1+z)}=\frac{x_e\mathrm{\Phi }(z)}{1+z}.$$
(15)
Here $`\sigma _T`$ is the Thomson cross-section and $`c`$ is the speed of light.

The redshift of the last scattering surface, $`z_r`$, is defined by the position of the maximum of the so–called visibility function, $`g(\tau _T)=\tau _T^{^{}}\mathrm{exp}(\tau _T)`$,

$$\frac{dg(\tau _T)}{dz}|_{z=z_r}=0,\tau _T^{^{\prime \prime }}(z_r)=\left(\tau _T^{^{}}(z_r)\right)^2,$$
(16)
and its width is

$$\mathrm{\Delta }_z\left(2\left(\tau _T^{^{}}(z_r)\right)^2\frac{\tau _T^{^{\prime \prime \prime }}(z_r)}{\tau _T^{^{}}(z_r)}\right)^{1/2}$$
(17)
(see, e.g., Naselsky & Polnarev 1987). Taking into account that recombination is a fast process and $`\mathrm{\Delta }_z0.1z_r`$, we can further consider only redshift variations of the ionization fraction, $`x_e`$, and take the function $`\mathrm{\Phi }`$ in (15) at $`z=z_r`$. This means that the parameters $`z_{}`$ and $`\mathrm{\Delta }_z(z_{})`$ can be found from the following equations:

$$x_e^{^{}}(z_{})x_e^2(z_{})\mathrm{\Phi }(z_r)/z_r,\mathrm{\Delta }_zx_e(z_{})/x_e^{^{}}(z_{})x_e^1(z_{})z_r/\mathrm{\Phi }(z_r).$$
(18)
Here $`z_r`$ and $`z_{}`$ are the redshifts of the last scattering surface in the standard and the distorted models. For small distortions of ionization fraction we have

$$x_ex_{e,st}+x_{e,st}^{^{}}z_r(1z_{}/z_r),x_e^{^{}}x_{e,st}^{^{}}n_b^1dn_\alpha /dz=x_{e,st}^{^{}}\epsilon _\alpha (z_r)/z_r,$$

$$z_{}z_r(1\epsilon _\alpha (z_r)\nu ),\nu =0.5\mathrm{\Phi }^2(z_r)x_{e,st}^3,$$
(19)
$$\mathrm{\Delta }_z(z_{})\mathrm{\Delta }_{z,st}[1+0.5\epsilon _\alpha (z_r)\mathrm{\Phi }^1(z_r)x_{e,st}^2].$$

Here $`x_{e,st}`$ and $`\mathrm{\Delta }_{z,st}`$ are the ionization fraction and the thickness of the last scattering surface for the standard model.

The leading order dependence of the first peak location is (Peebles, Seager and Hu 2001)

$$l_1(\epsilon _\alpha )l_{1,st}\sqrt{z_{}/z_r}l_{1,st}[10.5\nu \epsilon _\alpha (z_r)].$$
(20)
For the second peak of the CMB power spectrum the shift of the position is practically the same as for the first one, but the ratio between their amplitudes is a linear function of ratio $`z_{}/z_r`$. For the next few peaks growth of the thickness of the last scattering surface, $`\mathrm{\Delta }_z`$, increases the damping and decreases the amplitudes of the peaks (White 2001). In contrast, for models with accelerated recombination and $`\epsilon _\alpha 0`$, peaks are shifted to larger $`l`$ and the damping of perturbations decreases together with $`\mathrm{\Delta }_z`$ as compared with the standard model.

## 4 Five models of recombination and distortions of the ionization history.

### 4.1 Models of recombination.

To illustrate the possible manifestation of the external sources of resonance and ionized photons and non–uniform distribution of baryonic matter in the ionization history of the Universe we compare four models with different ionization history with the basic spatially flat $`\mathrm{\Lambda }`$CDM model, $`\mathrm{\Omega }_K=0`$, the standard ionization history, $`\epsilon _i=\epsilon _\alpha =0`$, and

$$\mathrm{\Omega }_\mathrm{\Lambda }=0.7,\mathrm{\Omega }_c+\mathrm{\Omega }_b=0.3,h=0.7,n_s=1,\mathrm{\Omega }_bh^2=0.022,\mathrm{\Omega }_ch^2=0.125.$$
(21)
This baryonic matter density is close to the most probable value found from the comparison of the Big-Bang nucleosynthesis prediction with the observable abundance of the light chemical elements. The basic model includes also the optical depth for the Compton scattering due to reionization of the Universe by young galaxies and quasars at redshift $`z_{rei}`$

$$\tau _{rei}=63x_e\left(\frac{1Y_p}{0.76}\right)\left(\frac{\mathrm{\Omega }_bh^2}{0.02}\right)\left(\frac{\mathrm{\Omega }_ch^2}{0.122}\right)^{\frac{1}{2}}\left(\frac{1+z_{rei}}{10^3}\right)^{\frac{3}{2}},$$
(22)
where $`x_e1`$ is the fractional ionization at the redshift $`z_{rei}`$ and $`Y_p`$ is the mass fraction of helium. For parameters (21), $`Y_p=0.24`$, $`z_{rei}14`$ and $`\tau _{rei}0.1`$, this model provides the best–fit for the recent CBI data. We will refer to the basic model as to Model 1 .

Recent observations of the farthest quasars (Fan et al. 2001) give the low limit for the epoch of reionization, $`z_{rei}6`$, which does not contradict the adopted value $`z_{rei}14`$. For $`z_{rei}6`$, the Thompson optical depth is only $`\tau _{rei}0.04`$. We will refer to the model with parameters (21), $`z_{rei}6`$ and $`\tau _{rei}0.04`$ as to the Model 1a and will discuss it in Section 6.

To demonstrate the influence of distortions of the ionization history, we consider three models with various distortions. In Model 2 we change the baryonic fraction and accept $`\mathrm{\Omega }_bh^2=0.03`$ preserving the total density, $`\mathrm{\Omega }_b+\mathrm{\Omega }_c=0.3`$. The Model 3 is the cloudy baryonic model with $`\mathrm{\Omega }_bh^2=0.022`$ and reionization of hydrogen at redshift $`z_{rei}14`$ as for the Model 1. We assume also that the density contrast between clouds, $`\rho _{in}`$, and interclouds medium, $`\rho _{out}`$, is $`\eta =\rho _{in}/\rho _{out}=11`$, clouds occupy the volume fraction $`f_v=0.1`$ and accumulate the mass fraction $`Z_m=\eta f_v\left[\eta f_v+(1f_v)\right]^1=0.55`$. For simplicity, we consider all clouds with the same mass in range $`10^2M_{}M_{cl}10^5M_{}`$ (Naselsky and Novikov, 2002).

The Model 4 differs from the Model 1 by action of external sources of ionization with

$$\epsilon _\alpha =\alpha \xi ^{\frac{3}{2}}\mathrm{exp}\left(\xi ^{\frac{3}{2}}\right),\epsilon _i=\beta \xi ^{\frac{3}{2}}\mathrm{exp}\left(\xi ^{\frac{3}{2}}\right),$$

$$\xi =(1+z)/(1+z_d)10^3(1+z).$$

The choice $`\alpha 0.3`$, $`\beta 0.13`$ coincides with the model of primordial black hole evaporation at the moment of recombination (Naselsky 1978; Naselsky & Polnarev 1987). In this Model we exclude the impact of late reionization which mimics roughly the model with $`z_{rei}6`$. For some applications we consider also the Model 5 which differ from the Model 4 by the baryonic density, $`\mathrm{\Omega }_bh^2=0.03`$ preserving the total density, $`\mathrm{\Omega }_b+\mathrm{\Omega }_c=0.3`$. In this respect, the Model 5 is identical to the Model 2.

Main parameters of the Models are listed in Table 1.

For all the models the anisotropy and polarization power spectra are found with modified RECFAST and CMBFAST codes. To characterize the differences between the models and to compare them with the expected sensitivity of the PLANCK experiment we describe the deviations of the anisotropy and polarization power spectrum in terms of the functions

$$D_{i,j}^a(l)=2\left[C_i^a(l)C_j^a(l)\right]/\left[C_i^a(l)+C_j^a(l)\right],$$

$$D_{i,j}^p(l)=2\left[C_i^p(l)C_j^p(l)\right]/\left[C_i^p(l)+C_j^p(l)\right],$$
(23)
where $`i,j`$ relate to the Models and $`C^a`$ and $`C^p`$ are the power spectrum of anisotropy and $`E`$ component of polarization, respectively. The differences between the ionization fractions in the models are characterized by the ratios

$$\mathrm{\Delta }x_{i,j}=2\frac{x_ix_j}{x_i+x_j}.$$
(24)
### 4.2 Distortions of the ionization history

The redshift variations of the relative difference between the ionization fractions in the basic Model 1 and three other models, $`\mathrm{\Delta }x_{1,j},j=2,3,4`$ are plotted in Fig. 2 together with the redshift variations of the ionization fraction in the basic Model 1.

For the Model 2, this difference does not exceed $`20\%`$ at $`z10^3`$ and increases up to $`30\%`$ at lower redshifts which demonstrates a relatively weak influence of the variations of baryonic fraction. For the Model 3 with accelerated recombination, the relative difference $`\mathrm{\Delta }x_{1,3}<0`$ at redshifts $`z700`$ and the ionization fraction decreases with time faster than in the basic model. For $`z700`$ the ionization fraction is larger than that in the basic model, $`\mathrm{\Delta }x_{1,3}>0`$, but the excess is less than $`10\%`$ and the distortions of the CMB anisotropy and polarization power spectra remain small due to the small optical depth.

The conspicuous delay of recombination and distortions of the ionization fraction appear for the Model 4. They achieve $`2030\%`$ already at $`z10^3`$ and are maximal at $`z500`$, which provides the essential growth of the optical depth for the Thomson scattering and suppresses the amplitude of the CMB anisotropy.

These results indicate that we can expect noticeable variations of the CMB anisotropy and polarization power spectra for the Model 3 and the Model 4 while for the Model 2 they can be quite moderate.

## 5 Anisotropy and polarization power spectra as a test of the ionization history of the Universe.

To obtain the power spectra of the CMB anisotropy and polarization in our models we have to use the modification of the CMBFAST code (Seljak and Zaldarriaga, 1996) taking into account more complicated ionization history of the Universe discussed above. For the reference Model 1 and 2, we use the standard CMBFAST code for the optical depth $`\tau _r`$ and the ionization fraction $`x_e=1`$ achieved at the redshift of reionization $`z_{rei}`$.

In Fig.3 we plot the CMB power spectrum for the models 1, 3 and 4 in comparison with the observational data from the BOOMERANG, MAXIMA-1, and CBI (Mason et al, 2002) experiments at the multipole range $`l`$ 2 000 where the possible influence of the Sunyaev-Zel’dovich clusters is not yet crucial.

As one can see from Fig. 3, for all models the power spectra are very similar to each other and the Model 1, 3 and 4 are consistent with the available measurements of $`C(l)`$. This means that the more complicated ionization history of hydrogen violates the standard fitting procedure. In particular, the comparison of results listed in Table 2 for the Model 1, 2, and 5 shows that the determination of baryonic density from the CMB measurements is ambiguous. This fact is important for the interpretation of the MAP and PLANCK measurements.

To characterize the differences between the Models and to compare them with the expected sensitivity of the PLANCK experiment we plot in Fig. 4 the functions $`D_{13}(l)`$ and $`D_{14}(l)`$ defined by (23) for the multipole range $`2l`$ 2000. These functions can be directly compared with the expected error bars of $`C(l)`$ for the PLANCK mission.

As is clear from Fig. 4 for the multipole range $`500l2000`$, $`0.05D_{13}(l)0.2`$ and $`0.05D_{14}(l)0.3`$, and there is a quite regular tilt at high multipole range. Obviously, the tilt is related to the damping of the $`C(l)`$ power spectrum in the basic model because of the late reionization with $`\tau _{rei}=0.1`$. In order to clarify the manifestations of the more complex ionization regimes in the Model 3 and 4 we have to compare the peak–to–peak variations of amplitudes, $`D_{13}(l)`$ and $`D_{14}(l)`$, with the expected errors of measurements for the PLANCK experiment. If the systematic effects will be successfully removed, the error should be close to

$$\frac{\mathrm{\Delta }C(l)}{C(l}\frac{1}{\sqrt{f_{sky}(l+\frac{1}{2})}}\left[1+w^1C^1(l)W_l^2\right]$$
(25)
$$w=(\sigma _p\theta _{FWHM})^2,W_l\mathrm{exp}\left[\frac{l(l+1)}{2l_s^2}\right],f_{sky}0.65.$$

Here $`f_{sky}`$ is the sky coverage during the first year of observations, $`\sigma _p`$ is the sensitivity per resolution element $`\theta _{FWHM}\times \theta _{FWHM}`$ and $`l_s=\sqrt{8\mathrm{ln}2}\theta _{FWHM}^1`$.

For example, for the HFI 353 GHz channel the $`\theta _{\mathrm{FWHM}}`$ is 5 arcmin and the error at $`l2000`$ are defined by the first term in Eq (25) . As one can see from Fig. 4, at $`l10001500`$ the peak–to–peak amplitudes are $`1020\%`$, while $`\mathrm{\Delta }C(l)/C(l)57\%`$. This means that such peculiarities of the anisotropy power spectrum can be detected by the PLANCK mission, if the systematic effects would be less than $`\mathrm{\Delta }C(l)/C(l)[f_{sky}(l+1/2)]^{1/2}`$.

The polarization power spectra are even more sensitive to distortions of the ionization history of the Universe. Recently the DASI experiment (Kovac et al. 2002) has detected the CMB polarization without significant manifestation of polarized foregrounds. It indicates that the launched MAP and the PLANCK missions will be able to measure the polarization of the CMB with unprecendent accuracy.

The polarization power spectrum $`\mathrm{\Delta }T_p^2(l)=l(l+1)C_p(l)/2\pi `$ are plotted in Fig.5 for the Model 1, 3 and 4 and their relative variations, $`D_{1,j}^p(l)`$, with respect to the basic model are plotted in Fig. 6 . As one can see the differences between spectra are quite small especially at $`l200`$. But the bump at $`l510`$ which is caused by the late reionization in the basic model can be essentially suppressed for the Model 4 wherein the required optical depth is gained at redshifts $`z100500`$.

For the PLANCK mission, even if all the systematic and foreground contaminations are removed from the polarization map (see Franco, Fosalba and Tauber 2002), one of the most important sources of uncertainty at $`l200`$ is the cosmic variance. For $`f_{sky}0.65`$ and $`l2000`$ the cosmic variance limit is close to $`\mathrm{\Delta }C(l)/C(l)3\%`$ for $`1\sigma `$ confidence level and all peculiarities of the polarization power spectrum should be observable. As it follows from Fig. 6, all polarization features which corresponds to the Model 1, 3 and 4 can be observed by the PLANCK mission, and, for the low multipole part $`l200`$, can be probably observed already by the MAP mission. This means that future polarization experiments can be crucial for investigation of the ionization history of the cosmic plasma as near the redshift of the last scattering surface as at the redshift of reionization.

## 6 Low multipole peculiarities of the polarization power spectrum in models with delayed recombination

The signature of late reionization is the bump at small multipoles, $`l10`$, in the polarization power spectrum. Its amplitude and shape depend upon the ionization history characterized by $`z_{rei}`$ and $`\tau _{rei}`$. This bump is clearly seen in Fig. 5 for the basic model with $`z_{rei}=14`$ and $`\tau _{rei}=0.1`$ but it is strongly suppressed for the later reionization such as, for example, for $`z_{rei}6`$ and $`\tau _{rei}=0.04`$ widely discussed in connection with the recent observations of the farthest quasars (Fan et al. 2001). For the non–standard ionization history discussed in this paper, the shape of the polarization power spectrum is more complicated but its observations allow to restore some features of the ionization hystory.

As is seen from Fig. 6, for the Model 3 with accelerated recombination this bump is just the same as in the basic Model 1. However, in the Model 4 with delayed recombination and without any sudden reionization at low redshifts significant variations of $`D_{1,j}^p(l)`$ are seen in the Fig. 6 at $`l500`$. They can be observed by the PLANCK experiment even for relatively small $`\epsilon _i(z)`$, which would clarify the ionization history for such models.

To illustrate these abilities, we plot in Fig. 7 the relative distortions of the polarization power spectrum at $`l1000`$ for several models with delayed recombination and different $`\epsilon _i(z)`$. First of all, to illustrate the influence of the optical depth, $`\tau _{rei}`$, we compare the Model 1 with the Model 1a with the same parameters and $`z_{rei}6`$ and $`\tau _{rei}=0.04`$. Large values $`D_{1,j}^p(l)`$ for these models at $`l20`$ allows to easily discriminate between these models.

In the same Fig. 7 we compare also the Model 1a with several models of delayed recombination with the same cosmological parameters (21) and with

$$\epsilon _i(z)=\epsilon _i(z_r)\left(\frac{1+z}{1+z_r}\right)^{2(1\frac{3}{4}p)}$$
(26)
where $`z_r10^3`$ is the redshift of the last scattering surface. We see that for the models with $`p=1,\&\epsilon _i(z_r)=10^2`$ and with $`p=2,\&\epsilon _i(z_r)=10^3`$ the relative distortions, $`D_{1,j}^p(l)`$, are quite similar and can be observed at $`l`$ 100 . However, for the model with $`p=1,\&\epsilon _i(z_r)=10^3`$ distortions become small and are under observational limits. These comparisons demonstrate the abilities of feature polarization measurements with respect to restrictions of the ionization history of the Universe at $`z10^3`$.

In wide class of such models the injection of ionized photons at $`zz_r`$ increases the ionization fraction and the Thompson optical depth with respect to the standard models. In such case, the already available limitations of the depth $`\tau (z)\tau _{thr}0.1`$ restrict the admissible $`\epsilon _i(z)`$ even at high redshifts where estimates (10) are not effective. To estimate the optical depth for such models we have to trace numerically the ionization fraction at all redshifts. However, even rough analytical estimates are informative and allow to reject some models.

For such rough estimates we use the equation of ionization balance

$$\alpha _cn_bx_e^2\epsilon _i(z)H(z)$$
(27)
where $`\alpha _cT_m^{0.6}`$ is the coefficient of hydrogen recombination and $`T_m`$ is the temperature of the matter. As is well known, $`T_m1+z`$ for $`zz_{eq}10^2`$ and $`T_m(1+z)^2`$ for lower redshifts and, so, $`\alpha _cT_m^{0.6}(1+z)^{2\omega }`$, $`\omega 0.3`$ for $`zz_{eq}`$ and $`\omega 0.6`$ for $`zz_{eq}`$ . Thus, for $`\epsilon _i(z)`$ given by (26) we get for the ionization fraction and for the optical depth

$$x_e(z)\epsilon _i^{\frac{1}{2}}(z_r)(1+z)^\mu ,\mu =0.75+\omega +(13p/4),$$

$$\tau _r36\nu ^1\epsilon _i^{1/2}(z_r)\left(\frac{1+z_\tau }{1+z_r}\right)^\nu \left(\frac{1Y_p}{0.76}\right)\left(\frac{\mathrm{\Omega }_bh^2}{0.02}\right)^{1/2}\left(\frac{\mathrm{\Omega }_ch^2}{0.125}\right)^{1/4}$$
(28)
where $`Y_p`$ is the mass fraction of helium, $`\nu =1.5+\mu =0.75+\omega +(13p/4)`$, and $`z_\tau `$ is the maximal redshift where the ionization balance is described by (27). This relation is directly applied for $`z_\tau z_{eq}`$, $`\omega =0.6`$ and can be easily extended for the case $`z_\tau z_{eq}`$.

The observations of farthest quasars show that the reionization occurs at $`z_{rei}6`$ when $`\tau 0.04`$ and, so, for cosmological parameters (21), $`Y_p=0.24`$ and $`\tau _{thr}0.1`$ we get

$$\epsilon _i(z_r)10^6\nu ^2\left(\frac{1+z_\tau }{1+z_r}\right)^{2\nu }.$$
(29)
As compared with the normalization (10), this expression more effectively restricts the possible energy injection at high redshifts.

## 7 Discussion and conclusion

The already achieved precision of measurements of the CMB anisotropy and especially expected one for the MAP and PLANCK missions allows to discriminate wide set of cosmological models and to reveal noticeable distortions of the standard model of hydrogen recombination process at redshifts $`z10^3`$ and of the ionization history. Many physical models can be considered as a basis for such discussions and can be tested and restricted via such measurements.

In this paper we compare available observations of the CMB anisotropy with those expected in several models with different kinetics of recombination. The differences can be caused by possible external sources of resonance and ionized radiation, possible small scale clustering of baryonic component and other factors. We show that the models with more complicate ionization history provide better description of the available observational data. This means also that the accepted estimates of cosmological parameters can be changed when such models are incorporated in the fitting procedure. It is especially crucial for the estimates of baryonic density, $`\mathrm{\Omega }_bh^2`$, which is closely related to the Big-Bang nucleosynthesis predictions.

In this paper we do not consider the impact of Sunyaev–Zel’dovich (SZ) effect and the weak lensing by the large scale structure elements. However, the SZ effect is frequency dependent that allows to reduce its contribution for multifrequency measurements. As was shown in Schmalzing, Takada and Futamase (2000), the influence of lensing effect can really distort the measured characteristics of the CMB only at larger $`l10^3`$ while the discussed deviations of ionizing history distort these characteristics at all $`l`$. This means that the measurements of the CMB spectra in wide range of $`l`$ allows to discriminate between these distortions. However, this problem requires more detailed discussions for various models of both lensing and recombination.

Near the first Doppler peak, the cosmic variance effect restricts an accuracy of estimates of $`C^a(l)`$ as $`\mathrm{\Delta }C^a(l)/C^a(l)10\%`$ for all experiments with the sky coverage $`f_{sky}65\%`$. For the recent balloon-borne and ground–based experiments the achieved accuracy in this range of $`l`$ is even worse than this limit. With such precision, we cannot discriminate the basic model and models with delayed recombination and another $`\mathrm{\Omega }_bh^2`$. To illustrate this statement we compare quantitatively the observational data with the power spectra for Model 1, 2 & 3 parameters which are given in Sec. 4.1 . Instead of the Model 4 we here present the results for Model 5 which differs from Model 4 by the values of parameters $`\mathrm{\Omega }_bh^2=0.03`$ and $`\mathrm{\Omega }_m=0.256`$ which are the same as in the Model 2.

The quality of models in comparison with the recent BOOMERANG, MAXIMA-1, CBI and VSA observational data can be characterized quantitatively by the $`\chi ^2`$-parameter listed in Table 2 for all models. Here the reference Model 1 provides the best– fit for the CBI data for the standard ionization history while $`\chi _2^2\chi _1^2`$ for the Model 2 with larger $`\mathrm{\Omega }_bh^2`$ . However, the Model 3 demonstrates an excellent agreement with the CBI data, $`\chi _3^2=0.5\chi _1^2`$, and $`\chi _3^2\chi _1^2`$ for other observations. For the Model 5 with $`\mathrm{\Omega }_bh^2=0.03`$ we get also $`\chi _5^2\chi _1^2`$ for the CBI data and only for the BOOMERANG data $`\chi _5^2\chi _1^2`$.

Obviously, such extension of the cosmological models increases the number of parameters using to fit observed power spectra of the anisotropy and, in perspective, of the polarization. We show that the influence of these ”missing” parameters, namely $`\epsilon _i(z)`$ and $`\epsilon _\alpha (z)`$ , can improve the fits used and, in particular, to obtain better agreement between the observed and expected positions and amplitudes of peaks for higher multipoles. At the same time, including of the ”missing” parameters changes the measured values of cosmological parameters usually discussed.

We show also that the expected sensitivity of the MAP and the PLANK missions in respect of the measurements of the polarization will allow to discriminate between main families of such models and, in particular, between models with small and large optical depth at redshifts $`10^3z`$ 20-50. We would like to note that the realistic values of the cosmological parameters can not be obtained from the CMB data without the PLANCK observations of the polarization.

## 8 Acknowledgments

This paper was supported by Denmark’s Grundforskningsfond through its support for an establishment of Theoretical Astrophysics Center and, in part, by the Danish Natural Science Research Council through grant No. 9701841, and grant RFBR 17625. We should like to thank the anonymous referee for useful comments.


format: text

text:
# DIFFRACTION AT RHIC

## 1 Introduction

Diffractive events occur via the exchange of color singlet objects, Pomerons, with the same quantum numbers as the vacuum ($`J^{PC}=0^{++}`$). $`^\mathrm{?}`$ They are characterized by final states which include rapidity gaps, regions of phase space containing no final state particles. Examples are $`pp`$ elastic scattering, photoproduction of vector mesons, and diffractive $`W^\pm `$ production. Typical hadronic interactions, in contrast, involve the exchange of colored particles, where the probability of finding a particle-free region with a width of $`\mathrm{\Delta }y`$ units of rapidity is exponentially suppressed as $`\mathrm{exp}(\mathrm{\Delta }ydN/dy)`$, where $`dN/dy`$ is the mean number of particles per unit rapidity. Meson exchange can also produce rapidity gaps. However, the meson exchange contribution drops with increasing collision energy, while the Pomeron exchange cross section increases slightly. So, at high (RHIC) energies, processes with rapidity gaps are expected to be dominated by Pomeron exchange.

Despite 40 years of studying diffraction, many questions remain. Soft (low energy) diffraction is usually characterized in terms of the optical model as represented by the absorptive part of the cross section, and via a Regge trajectory. The optical model can describe phenomena such as elastic scattering and vector meson photoproduction, but does very poorly with hard reactions like $`W^\pm `$ and jet production. These higher momentum-transfer reactions may be explainable within perturbative QCD as via Pomeron exchange, $`^\mathrm{?}`$ where the Pomeron is represented as a gluon ladder (two gluons connected by additional gluon and quark loop ‘rungs’). These models are moderately successful in explaining many hard diffraction phenomena.

Despite these successes, many questions and controversies remain. In $`pp`$ collisions, it is still unclear whether diffractive jet production can be explained as due to the collisions of two Pomerons. Can the Pomeron flux from a relativistic proton be described in terms of a Pomeron distribution, akin to the usual parton distributions? Another example involves photoproduction. At HERA, the $`J/\psi `$ and $`\psi ^{}`$ cross sections rise sharply with energy, $`^\mathrm{?}`$ in contrast to the expected slow energy dependence of the Pomeron exchange model. This has been argued to be a threshold effect, or due to a hard component of the Pomeron.

The nature of the Pomeron is also debated. In the simplest models, it is almost entirely gluonic. However, the observation of diffractive deep inelastic scattering $`^\mathrm{?}`$ and diffractive $`W^\pm `$ production in $`p\overline{p}`$ collisions $`^\mathrm{?}`$ both require a significant quark component. The role of the 3-gluon triangle (Fig. 1b) is also debated. It is also difficult to explain hard and soft diffraction within a single theory; $`^\mathrm{?}`$ some papers argue that they are two distinct phenomena. $`^\mathrm{?}`$ For the soft Pomeron, this is little understanding of how diffraction (i.e. the absorptive part of the cross section) meshes with colored QCD reactions.

RHIC can bring significant new data to bear on many of these questions. One key feature of RHIC is its polarized proton beams, allowing for the first time studies of polarized hard diffraction up to center of mass energies of 500 GeV. As the first high energy $`pp`$ collider, RHIC can also look for differences between $`pp`$ and $`p\overline{p}`$.

These studies will complement studies of $`pp`$ elastic scattering by the $`pp2pp`$ collaboration. $`^\mathrm{?}`$ The $`pp2pp`$ collaboration will measure single and double spin asymmetries in elastic scattering, and so probe the spin structure of the Pomeron at very soft momentum scales.

Diffraction is also a tool to study meson spectroscopy. Double-diffractive (two rapidity gaps) central production of mesons has been studied extensively at the CERN SPS. $`^\mathrm{?}`$ The experiments found that invariant mass spectra changed drastically when the protons were ’kicked’ in opposite directions, compared to when the protons went in the same direction. This ’$`p_T`$ filter’ is still poorly understood, but is clearly important for both meson spectroscopy and for studies of diffraction. RHIC studies of the spin dependence of the $`p_T`$ filter could yield qualitatively new information.

Events with rapidity gaps can also be a place to search for new physics. $`^\mathrm{?}`$ The new physics can occur via either double-Pomeron or two-photon interactions. The two event classes look similar, except that in two-photon interactions, the outgoing protons have a smaller average $`p_T`$; the classes may be statistically separated on this basis. $`^\mathrm{?}`$

Diffractive $`pA`$ collisions can be used to study the nucleus. If the Pomeron is made up of two gluons, then the Pomeron should exhibit nuclear shadowing (the EMC effect), with a magnitude roughly the square of the size of gluon shadowing. $`^\mathrm{?}`$

Diffraction usually requires that the interacting protons remain intact. Dedicated forward detectors called Roman pots are commonly used to observe the outgoing protons. Roman pots detect protons travelling inside the beam pipe of the accelerator, making it possible to measure very small scattering angles. Since the Roman pots are usually placed past the accelerator magnets, they can act as spectrometers, measuring the kinematics of the scattered protons. At RHIC, diffraction could be studied by adding two Roman pot systems to a central detector. Because of the angular sensitivity of the $`p_T`$ filter and all spin measurements, pots with full azimuthal acceptance are very desirable.

This physics was discussed at a workshop on “Diffraction and Glueball Searches at RHIC,” at Brookhaven National Laboratory, May 17-18, 2002. $`^\mathrm{?}`$ The workshop found ample physics justification to install a Roman pot system around one of the existing detectors. This writeup will now present three core topics discussed at the workshop: meson spectroscopy, hard diffraction in $`pp`$ and diffraction in $`pA`$.

## 2 Meson Production in Soft $`pp`$ Diffraction

Mesons can be produced diffractively via Pomeron-Pomeron fusion, $`ppppPPppX`$, Figure 1. Pomeron-meson and meson-meson fusion are also possible. At lower energies, the meson contributions are important; at RHIC, they should be small. Because Pomerons are believed to have a largely gluonic content, double-Pomeron collisions may be a good channel to search for glueballs.

Double-Pomeron production of mesons has been studied extensively by experiments WA 76/91/102 at the CERN SPS. $`^\mathrm{?}`$ These experiments made a remarkable discovery: the character of the produced mesons depended strongly on the direction of the momentum transfer from the protons. When the transverse momentum transfers pointed in the same direction, one $`K^+K^{}`$ (for example) invariant mass spectrum was obtained, showing mostly conventional $`q\overline{q}`$ mesons. When the momentum transverse pointed in opposing directions, the $`K^+K^{}`$ spectrum was very different, with many suspected ’exotic’ (non $`q\overline{q}`$) states appearing clearly.

The difference can be seen by selecting events on the basis of $`dp_T`$, the $`p_T`$ difference between the two momentum transfers. Figure 2 shows two example of the spectral changes in different $`dP_T`$ bins. Conventional mesons are prominent at high $`dP_T`$, while at low $`dP_T`$, ’unconventional’ mesons appear clearly. Partial wave analyses were used to confirm the spin/parity assignments. Similar spectral changes have been observed for many other channels, and the statistical significance of the effect is unassailable.

Many theoretical interpretations have been given for this effect. It may be evidence that the Pomeron is a vector particle, or at least acts as a non-conserved vector current, $`^\mathrm{?}`$ in contrast to the expected spin 0 behavior. It has also been cited as evidence for a QCD scale anomaly, $`^\mathrm{?}`$ and for the existence of instantons. $`^\mathrm{?}`$ The $`p_T`$ filter could also be due to a large contribution from the gluon triangle diagram, Fig. 1b. Alternately, it may be that Reggeon-Reggeon (meson-meson) and meson-Pomeron interactions are still important at SPS energies.

RHIC should be able to shed much light on this effect. If the $`p_T`$ filter effect is due to meson exchange, then it should be much reduced at RHIC. If the Pomeron acts as a vector particle, then it is likely to display spin effects which should be visible in polarized collisions.

The cross section for double-diffractive interactions at RHIC is about 100-200 $`\mu `$b, $`^\mathrm{?}`$ corresponding to an interaction rate of 20 kHz, even at a a center of mass energyy of 200 GeV and a luminosity of $`10^{31}`$cm-2s-1. RHIC will produce 100’s of millions of diffractive final states in a one month run, and any rate limitation is likely to be in the detector or trigger. The production rates are very high even for rather specific final states.

## 3 Hard $`pp`$ Diffraction

Precision measurements of hard diffraction in $`ep`$ and $`p\overline{p}`$ scattering at the HERA and the Tevatron, respectively, raise many questions about factorization and the universality of the hard Pomeron. Measurements of the diffractive structure function $`F_2^{D(3)}`$ $`^\mathrm{?}`$ at HERA and associated studies of the diffractive-induced hadronic final state $`^\mathrm{?}`$ imply a universal Pomeron composed predominantly of gluons. However, efforts to apply the next-to-leading order QCD fits from $`ep`$ data to jet production in $`p\overline{p}`$ scattering have not been entirely successful.

The structure of the Pomeron can be measured in terms of $`\beta `$, the fraction of the total Pomeron momentum carried by the interacting parton. The $`\beta `$ dependence of the diffractive structure function $`F_{JJ}^D`$, as derived from jet measurements at CDF, $`^\mathrm{?}`$ is shown in figure 3 for a given region of $`\xi `$, the fraction of the proton momentum carried by the Pomeron. A comparison is made with parameterizations of the diffractive structure function measured at HERA. The $`\beta `$ dependence is similar although the overall normalization is inconsistent between the two scattering environments. Possible reasons for this include the need for a modified Pomeron flux at Tevatron energies owing to unitarity constraints or poorly understood gap-destroying multiple interactions present in hadronic collisions.

The study of $`pp`$ scattering at RHIC can extend this field. With center of mass energies up to several hundred GeV, RHIC covers a kinematic region of diffractive scattering which overlaps with HERA and the Tevatron. If Roman Pot detectors were installed around one of the RHIC experiments then proton tagging would allow the study of the $`t`$dependence of the diffractive exchange. A better understanding of proton dissociation would help to reduce one of the dominant systematic errors in diffraction measurements.

Further, measurements of single and double spin asymmetries will probe the spin structure of the hard Pomeron from a polarized proton. The cross sections for most diffractive processes in $`p\overline{p}`$ collisions are usually a small percentage (typically around 1%) of the cross section of the corresponding non-diffractive final state.

Nothing is known about the spin sensitivity of this production. Non-diffractive $`W^\pm `$, $`Z^0`$ and jet production comes largely from $`q\overline{q}`$ annihilation and is expected to show significant polarization asymmetries; $`^\mathrm{?}`$ diffraction might show similar asymmetries. However, if diffraction is mediated by spin-0 Pomerons, then, to lowest order, no asymmetry is expected. At lower energies, there is some evidence that diffractive reactions may be sensitive to spin. $`^\mathrm{?}`$ Either case would be very interesting.

In addition to measurements sensitive to the hard scales of the scattering such as jet production, the RHIC detectors can track particles down to low $`p_T`$ in high-multiplicity events; these abilities will allow for high-quality measurements of particle multiplicities. Studies of multiplicity fluctuations have been used in the past to study the gluonic nature of the Pomeron. $`^\mathrm{?}`$

The rates for hard diffractive events may be scaled from the rates for corresponding hard non-diffractive events. Here, we give the rates for a 1 month ($`10^6`$ s) $`pp`$ run at the maximum 500 GeV center of mass energy, with an assumed luminosity of $`10^{32}`$cm-2s-1. The D0 collaboration finds that the rate of diffractive to non-diffractive jet events is $`1.07\pm 0.10_{0.13}^{+0.25}\%`$$`^\mathrm{?}`$ Assuming that this ratio holds, RHIC will produce about 2 million diffractive events per year containing a jet with pseudorapidity $`|\eta |<1`$ and transverse energy $`E_T>30`$ GeV. $`^\mathrm{?}`$

The CDF collaboration finds that diffractive $`W^\pm `$ production is $`1.15\pm 0.55\%`$ of non-diffractive $`W^\pm `$ production. $`^\mathrm{?}`$ With this ratio, each year RHIC will also produce 700 diffractive $`W^+e^+\nu `$, 20 $`W^{}e^{}`$, and 4 $`Z^0e^+e^{}`$, with lepton rapidity $`|\eta |<1`$$`^\mathrm{?}`$ The jet rates are very high, allowing for precision studies of single and double spin asymmetries and the like. Even diffractive vector boson production should be clearly visible.

## 4 Diffraction with Heavy Ions

The $`pA`$ program at RHIC enables the first studies of high energy diffraction with large and varying target sizes to take place.

Diffraction in $`pA`$ collisions can be an important probe of the gluon density of heavy nuclei at low-$`x`$, where gluon densities are expected to be very high, and, even with shadowing, may saturate the available phase space. $`^\mathrm{?}`$ In this regime, it may be possible to describe the gluons as a classical field, and new effects, such as a colored glass condensate may be observed. $`^\mathrm{?}`$

Diffractive processes are sensitive to the total absorptive content of the nucleus, and, through that, to the gluon content of the medium. Several specific processes have been proposed to study the gluon density at small $`x`$. For example, one could compare the rate of forward jet production in $`pA`$ to that in $`pp`$ to study nuclear modifications of the gluon distribution.

The $`t`$dependence of the diffractive exchange depends critically on the target size. The $`t`$ and $`A`$ variation of diffractive cross sections could be used to infer the relative contributions of soft and hard diffraction and their variation with $`A`$.

The breakdown in universality of the picture of the partonic Pomeron developed in deep-inelastic $`ep`$ scattering could be probed as a function of $`A`$. Using $`pp`$ scattering as a baseline, $`A`$dependent deviations of multi-jet production and multiplicity fluctuations from the simple partonic picture could be studied. This should be sensitive to low-$`x`$ gluon shadowing and gap-destroying multiple interactions.

These studies complement the existing RHIC program for diffraction in $`AA`$, in the form of ultra-peripheral photon-Pomeron and two-photon interactions. $`^\mathrm{?}`$ The cross sections for photon-Pomeron fusion into vector mesons are very large, $`^\mathrm{?}`$ and $`\rho ^0`$ production is observed at the predicted level. $`^\mathrm{?}`$ Ultra-peripheral collisions may be used to study nuclear shadowing; production of both open heavy quarks $`^\mathrm{?}`$ and quarkonium $`^\mathrm{?}`$ is sensitive to the gluon distribution in heavy nuclei.

## 5 Conclusions

With its polarized beams, high luminosity, and high-quality detectors, RHIC is capable of making significant contributions to our understanding of diffraction, in the areas of meson spectroscopy, hard diffraction, and the low-$`x`$ gluon distributions in nuclei. Much of this physics depends on measurements of the outgoing scattered protons. The proton kinematics could be determined by adding Roman pot systems upstream and downstream of a central detector at RHIC. We thanks Andrew Kirk for permission to use Fig. 2. This work was supported by the U.S. Department of Energy under Contract No. DE-AC-03076SF00098.

## References


text:
# Substructure and the halo model of large-scale structure

## 1 Introduction

Sheth & Jain (1997) described how the halo model for clustering can allow one to model the distribution of matter in the highly nonlinear regime. The model falls within the broader framework described by Neyman & Scott (1954) and Scherrer & Bertschinger (1991). It combines results from Peebles (1974) and McClelland & Silk (1977) with the work of Press & Schechter (1974), and is able to provide a good description of nonlinear clustering seen in numerical simulations of hierarchical gravitational clustering. The halo model assumes that most of the mass in the Universe is bound up in virialized dark matter halos, and that statistical measures of clustering on small scales are dominated by the internal structure of the halos. The agreement with simulations shows that is possible to provide an accurate description of clustering in the small-scale nonlinear regime even if one has no knowledge of if and how the halos themselves are clustered. This halo–model of clustering has been the subject of much recent interest (e.g. Seljak 2000; Peacock & Smith 2000; Ma & Fry 2000; Scoccimarro et al. 2001; Cooray & Sheth 2002).

To date, almost all analytic work based on the halo–model approach assumes that halos are spherically symmetric, and that the density run around each halo center is smooth. Halos which form in numerical simulations of hierarchical clustering are neither spherically symmetric nor smooth (e.g., Navarro, Frenk & White 1997; Moore et al. 1999; Jing & Suto 2002). About ten percent of the mass of a halo is associated with subclumps (Tormen, Diaferio & Syer 1998; Ghigna et al. 1999). The main purpose of the present work is to derive a model which accounts for this substructure.

Section 2 shows how to compute two-point statistics, the correlation function and its Fourier transform, the power spectrum, when substructure is important. Section 3 shows that the model can be easily extended to estimate higher-order statistics. Section 4 shows how our formalism can incorporate a range of parent halo and subclump masses, and Section 5 provides a few explicit examples. This section includes a discussion of how to model the shape of the subclump mass function. Section 6 summarizes our results, and suggests various other applications of our formalism.

## 2 Two-point correlations

The halo model approach assumes that all mass is bound up in dark matter halos. All statistical measures of clustering are then decomposed into two distinct types of contributions: one comes from sets of particles which are in the same halo, and the other comes from particles which are in different halos. The assumption is that, on the smallest scales, it is the single-halo contribution which dominates the statistic. In the context of the present paper, this means that we expect that the contribution to, e.g., the power spectrum, which comes from substructure cannot be important on scales larger than a typical halo. Therefore, even if substructure changes the single-halo contribution substantially compared to the case of smooth halos, it should be accurate to ignore substructure when estimating the two- and higher-order halo terms which dominate on larger scales. The primary consequence of substructure, then, is that the single-halo term becomes more complicated. For example, for two-point statistics, the one-halo term will now have three types of pairs: both particles from the smooth component, both particles from the substructure component, and one particle from the smooth component with the other from the substructure. Our goal will be to derive expressions for these different contributions.

### 2.1 Preliminaries

The probability of finding a subclump in the volume element $`d^3𝒓`$ at distance $`𝒓`$ from the center of the parent halo of mass $`M`$ is

$$p(𝒓)d^3𝒓=\frac{n_c(𝒓)d^3𝒓}{d^3𝒓n_c(𝒓)}.$$
(1)
Here $`n_c`$ describes the number density of subclumps at $`𝒓`$. In what follows, it will prove convenient to multiply this number density by a mass and so define a mass density $`\rho _c`$. We choose this mass so that $`d^3𝒓\rho _c(𝒓)=M`$, the total mass in the halo. (For simplicity we will sometimes assume that the subclumps have the same density run as the smooth dark matter distribution, but our analysis is not confined to this case. In this special case, $`\rho _c`$ is the same as the mass density profile of the smooth component.)

The density at $`𝒓`$ from the halo center is the sum of the smoothly distributed mass plus the contribution from the subclumps. If $`F`$ and $`f_i`$ denote the density run in units of the mean density $`\overline{\rho }`$ around the parent halo and the $`i`$th subclump respectively, then

$$\frac{\rho (𝒓)}{\overline{\rho }}=F(𝒓)+\underset{i=1}{\overset{N}{}}f_i(𝒓𝒓_i).$$
(2)
The mass in the smooth and subclump components is defined by

$$M_s𝑑r\mathrm{\hspace{0.17em}4}\pi r^2\overline{\rho }F(r)\mathrm{and}m_i𝑑r\mathrm{\hspace{0.17em}4}\pi r^2\overline{\rho }f(r).$$
(3)
The total mass $`M`$ is the sum of the smooth and clumped components.

Our strategy will be to first derive expressions for, e.g., the correlation function for fixed values of $`M`$, $`m_i`$ and $`N`$, and to average over the distributions of these variables later. We will keep the discussion as general as possible; realistic choices for $`F`$, $`f`$, $`p`$, and the halo and subclump mass functions will be inserted into the formalism later.

### 2.2 The correlation function of the mass

The ensemble-averaged two-point correlation function $`\xi (r)`$ will be obtained by averaging over halos with distributions of sub-structure specified by $`p(r)`$. We will first obtain the two-point correlation function within a halo, denoted $`𝒞(𝒓)`$, then average over substrucutre to get $`\xi (r)`$, and finally in section 4 consider the averaging over varying numbers of subclumps and of parent halo masses. For a halo of total mass $`M`$ we define

$$𝒞(𝒓)\overline{n}d^3𝒔\delta (𝒔)\delta (𝒔+𝒓),\mathrm{where}\delta (𝒓)\frac{\rho (𝒓)}{\overline{\rho }}1,$$
(4)
and $`\overline{n}`$ denotes the number density of the parent halos. If the mass densities within the subclumps are all much larger than the mean density $`\overline{\rho }`$, then $`\delta (𝒓)\rho (𝒓)/\overline{\rho }`$. Writing out all the terms explicitly shows that

$`{\displaystyle \frac{𝒞(𝒓)}{\overline{n}}}`$ $`=`$ $`{\displaystyle d^3𝒔F(𝒔)F(𝒔+𝒓)}`$ (5)
$`+{\displaystyle \underset{i}{}}{\displaystyle d^3𝒔f_i(𝒔𝒓_i)f_i(𝒔+𝒓𝒓_i)}`$
$`+{\displaystyle \underset{j}{}}{\displaystyle \underset{ij}{}}{\displaystyle d^3𝒔f_i(𝒔𝒓_i)f_j(𝒔+𝒓𝒓_j)}`$
$`+{\displaystyle \underset{i}{}}{\displaystyle d^3𝒔F(𝒔+𝒓)f_i(𝒔𝒓_i)}`$
$`+{\displaystyle \underset{i}{}}{\displaystyle d^3𝒔F(𝒔)f_i(𝒔+𝒓𝒓_i)}.`$
This is the correlation function associated with a smooth component and the $`i`$ subclumps at the specified positions. The first term is the contribution from particles which are not in the subclumps, and is familiar from previous work; this is the contribution from the first pair-type. The second term is from pairs where both particles are in the same subclump, whereas the third term is from pairs where the two particles are in separate subclumps. These two terms represent the contribution from the second pair-type. The fourth and fifth terms are from pairs in which one particle is in a subclump and the other is in the smoother component.

We are less interested in the correlation function associated with a specific realization of the subclump distribution, than we are with what happens upon averaging over the various possible subclump distributions. That is, we are more interested in

$$\xi (𝒓)=d^3𝒓_id^3𝒓_jp(𝒓_i)p(𝒓_j)𝒞(𝒓).$$
(6)
The next step is to compute these averages.

Let $`\lambda _{ss}(𝒓)`$ denote the value of the first integral in expression (5). It is the same for all the possible subclump distributions, so the result of averaging it over the distributions is

$$\lambda _{ss}(𝒓)=d^3𝒔F(𝒔)F(𝒔+𝒓).$$
(7)
This term is the convolution of $`F`$ with itself, and is the only term which most halo–models use. The remaining terms are due to the substructure. Before we compute them, it will prove convenient to rewrite this first term as

$$\lambda _{ss}(𝒓)=\left(\frac{M_s}{M}\right)^2\lambda _{\mathrm{smooth}}(𝒓)=\left(1\underset{i}{}\frac{m_i}{M}\right)^2\lambda _{\mathrm{smooth}}(𝒓).$$

Here $`\lambda _{\mathrm{smooth}}(𝒓)`$ denotes the correlation function if the total mass were smoothly distributed around the center of the parent halo (i.e, if there were no subclumps). Thus, the factor of $`(M_s/M)^2`$ simply denotes the fact that now only a fraction of the mass is in the smooth component.

Let $`\lambda _{ii}(𝒓)`$ denote the result of averaging each of the second terms in equation (5) over $`p(𝒓_i)`$. Then

$`\lambda _{ii}(𝒓)`$ $`=`$ $`{\displaystyle d^3𝒓_ip(𝒓_i)d^3𝒔f_i(𝒔𝒓_i)f_i(𝒔+𝒓𝒓_i)}`$ (8)
$`=`$ $`{\displaystyle d^3𝒙f_i(𝒙)f_i(𝒙+𝒓)},`$
where we have set $`𝒙=𝒔𝒓_i`$. This term is the convolution of each subclump profile $`f_i`$ with itself.

Averaging the third term in equation (5) is more complicated so we will do it last. The result of averaging the fourth term is

$`\lambda _{si}(𝒓)`$ $`=`$ $`{\displaystyle d^3𝒓_ip(𝒓_i)d^3𝒔F(𝒔+𝒓)f_i(𝒔𝒓_i)}`$ (9)
$`=`$ $`{\displaystyle d^3𝒔d^3𝒓_if_i(𝒔𝒓_i)F(𝒔+𝒓)\frac{\rho _c(𝒓_i)}{d^3𝒓\rho _c(𝒓)}}`$
$`=`$ $`{\displaystyle d^3𝒙f_i(𝒙)d^3𝒓_iF(𝒙+𝒓_i+𝒓)\frac{\rho _c(𝒓_i)}{M}}`$
$`=`$ $`{\displaystyle \frac{d^3𝒙}{M/\overline{\rho }}f_i(𝒙)\lambda _{sc}(𝒙+𝒓)}`$
where $`\lambda _{sc}`$ denotes the convolution of $`F`$ with $`\rho _c/\overline{\rho }`$. (If the subclumps were distributed around the halo center similarly to the dark matter, i.e., $`\rho _c/\overline{\rho }F`$, then $`\lambda _{sc}\lambda _{ss}`$.) By symmetry, the average of the fifth term is the same.

And finally, the average of the third term is

$`\lambda _{ij}(𝒓)`$ $`=`$ $`{\displaystyle d^3𝒓_id^3𝒓_jp(𝒓_i)p(𝒓_j)}`$ (10)
$`\times {\displaystyle }d^3𝒔f_i(𝒔𝒓_i)f_j(𝒔+𝒓𝒓_j),`$
where the integral over $`𝒔`$ is the convolution of the two subclump profiles. Written this way, $`\lambda _{ij}`$ is the weighted sum of this convolution over all pairs of subclump positions. It is interesting to re-arrange the order of the integrals above:

$`\lambda _{ij}(𝒓)`$ $`=`$ $`{\displaystyle d^3𝒓_id^3𝒓_jp(𝒓_i)p(𝒓_j)}`$ (11)
$`\times {\displaystyle }d^3𝒔f_i(𝒔𝒓_i)f_j(𝒔+𝒓𝒓_j)`$
$`=`$ $`{\displaystyle d^3𝒔d^3𝒓_if_i(𝒔𝒓_i)}`$
$`\times {\displaystyle }d^3𝒓_jf_j(𝒔𝒓_i+𝒓+𝒓_{ij})p(𝒓_{ij}+𝒓_j)p(𝒓_j)`$
$`=`$ $`{\displaystyle d^3𝒙d^3𝒓_{ij}f_i(𝒙)f_j(𝒙+𝒓+𝒓_{ij})\frac{\lambda _{cc}(𝒓_{ij})}{(M/\overline{\rho })^2}}`$
$`=`$ $`{\displaystyle d^3𝒙d^3𝒚f_i(𝒙)f_j(𝒚)\frac{\lambda _{cc}(𝒚𝒙𝒓)}{(M/\overline{\rho })^2}}`$
where $`\lambda _{cc}`$ denotes the convolution of $`\rho _c/\overline{\rho }`$ with itself. The second line follows from setting $`𝒓_{ij}=𝒓_i𝒓_j`$, the third line from setting $`𝒙=𝒔𝒓_i`$ and the final expression from setting $`𝒚=𝒙+𝒓+𝒓_{ij}`$. When written in this way, $`\lambda _{ij}`$ appears very similar to the two-halo term when substructure is absent; $`\lambda _{cc}`$ plays the role of the ‘subclump–subclump correlation function’ (compare equation 2 in Sheth et al. 2001).

Using the above results, the correlation function $`\xi (r)`$ defined in equation (6) can be expressed as

$$\frac{\xi (𝒓)}{\overline{n}}=\lambda _{ss}(𝒓)+\underset{i}{}\lambda _{ii}(𝒓)+2\underset{i}{}\lambda _{si}(𝒓)+\underset{i}{}\underset{ji}{}\lambda _{ij}(𝒓).$$
(12)
Note that all the $`\lambda `$ terms on the right-hand side have dimensions of volume, so that $`\xi `$ is dimensionless.

### 2.3 Limiting cases for small-sized subclumps

Now consider some simple limiting cases. If the subclumps are much smaller than the parent halo, then we can treat them as point masses when performing the integrals which define $`\lambda _{si}`$ and $`\lambda _{ij}`$. The associated delta function profiles simplify the integrals, so that

$`{\displaystyle \frac{\xi (𝒓)}{\overline{n}}}`$ $``$ $`\lambda _{ss}(𝒓)+{\displaystyle \underset{i}{}}\lambda _{ii}(𝒓)+2{\displaystyle \underset{i}{}}{\displaystyle \frac{m_i}{M}}\lambda _{sc}(𝒓)`$ (13)
$`+{\displaystyle \underset{i}{}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{m_i}{M}}{\displaystyle \frac{m_j}{M}}\lambda _{cc}(𝒓).`$
If in addition, the distribution of subclumps around the halo center is similar to the dark matter, $`\rho _cF`$, then we can set $`\lambda _{sc}\lambda _{ss}`$ and $`\lambda _{cc}\lambda _{ss}`$. In this case, the previous expression becomes

$`{\displaystyle \frac{\xi (𝒓)}{\overline{n}}}`$ $``$ $`\lambda _{ss}(𝒓)\left[1+2{\displaystyle \underset{i}{}}{\displaystyle \frac{m_i}{M_s}}+{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}{\displaystyle \frac{m_im_j}{M_s^2}}\right]`$ (14)
$`+{\displaystyle \underset{i}{}}\lambda _{ii}(𝒓)`$
$`=`$ $`\lambda _{\mathrm{smooth}}(𝒓)\left[1{\displaystyle \underset{i}{}}\left({\displaystyle \frac{m_i}{M}}\right)^2\right]+{\displaystyle \underset{i}{}}\lambda _{ii}(𝒓).`$
(Recall that $`\lambda _{\mathrm{smooth}}`$ denotes the correlation function if all the mass was smoothly distributed around the halo center.) If the fraction of the total mass which is in subclumps is small, then $`\xi (𝒓)\overline{n}\lambda _{\mathrm{smooth}}(𝒓)+_i\overline{n}\lambda _{ii}(𝒓)`$. In this limit, the total correlation function is well approximated by taking what one would have got if the mass was smoothly distributed, and then adding the contribution from the individual subclump components. In the small mass and size limit, the contribution from the subclumps can only be important on scales smaller than the typical subclump, so that this additional contribution is only important on very small scales. On scales larger than the diameter of a typical subclump, the correlation function looks just as though the mass within halos is smoothly distributed. This simple assumption is probably sufficiently accurate for most applications.

We remarked that the term $`\lambda _{cc}(𝒓)`$ could be thought of as the correlation function of the subclumps. To see why, suppose there is no mass in the smooth component, and all the subclumps are infinitesimally small: i.e., we replace all factors of the subclump density profile $`f_i`$ with delta functions $`\delta _\mathrm{D}(𝒓_i)`$. Then, when $`𝒓0`$, only the third term in equation (5) contributes any pairs. Using the delta functions reduces this term to

$`\lambda _{ij}(𝒓)`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}{\displaystyle d^3𝒓_id^3𝒓_jp(𝒓_i)p(𝒓_j)}`$ (15)
$`\times {\displaystyle }d^3𝒔f_i(𝒔𝒓_i)f_j(𝒔+𝒓𝒓_j)`$
$`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}\left({\displaystyle \frac{m_im_j}{\overline{\rho }^2}}\right){\displaystyle d^3𝒓_ip(𝒓_i)p(𝒓_i+𝒓)},`$
which is proportional to the convolution of the subclump distribution $`\rho _c`$ with itself.

### 2.4 The cross–correlation between subclumps and mass

We can use a similar argument to compute the cross correlation between subclumps and mass. That is, we imagine we sit at the center of the $`i`$th subclump and we then compute the typical density at distance $`𝒓`$ from it. Because the first term in equation (5) is not centered on a subclump it no longer contributes. Since our constraint only requires one member of each pair of positions to be centered on a subclump we now need to replace one factor of $`f_i`$ with a delta function. This means

$`\overline{n}\lambda _{ii}(𝒓)`$ $``$ $`f_i(𝒓),\overline{n}\lambda _{si}{\displaystyle \frac{\lambda _{sc}(𝒓)}{(M/\overline{\rho })}},\mathrm{and}`$
$`\overline{n}\lambda _{ij}(𝒓)`$ $``$ $`{\displaystyle d^3𝒚f_j(𝒚)\frac{\lambda _{cc}(𝒚r)}{(M/\overline{\rho })^2}},`$
so that the subclump–mass cross correlation function is

$`\xi _\times (𝒓)`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{\lambda _{sc}(𝒓)}{M/\overline{\rho }}}+{\displaystyle \underset{i}{}}f_i(𝒓)`$ (16)
$`+{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}{\displaystyle d^3𝒔f_j(𝒔)\frac{\lambda _{cc}(𝒔𝒓)}{(M/\overline{\rho })^2}}.`$
In the limit in which the subclumps are much smaller in size than the parent halo the last term simplifies to give,

$$\xi _\times (𝒓)\underset{i}{}\frac{\lambda _{sc}(𝒓)}{M/\overline{\rho }}+\underset{i}{}f_i(𝒓)+\underset{i}{}\underset{ji}{}\frac{m_j}{M}\frac{\lambda _{cc}(𝒓)}{(M/\overline{\rho })}.$$

When $`\rho _c/\overline{\rho }F`$ this becomes

$$\xi _\times (𝒓)\underset{i}{}\frac{\lambda _{ss}(𝒓)}{M/\overline{\rho }}\left[1+\underset{ji}{}\frac{m_j}{M}\right]+\underset{i}{}f_i(𝒓),$$
(17)
and further, if the total mass in subclumps is small we obtain

$$\xi _\times (𝒓)\underset{i}{}\lambda _{ss}(𝒓)/(M/\overline{\rho })+\underset{i}{}f_i(𝒓).$$
(18)
It is easy to verify that the result for $`\xi _\times `$ can be obtained simply by averaging

$`{\displaystyle \underset{i}{}}\rho _{cm|i}(𝒓_i+𝒓)`$ $`=`$ $`{\displaystyle \underset{i}{}}F(𝒓_i+𝒓)+{\displaystyle \underset{i}{}}f_i(𝒓)`$ (19)
$`+{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}f_j(𝒓_i+𝒓𝒓_j)`$
over all realizations of subclump positions, $`p(𝒓_i)`$ and $`p(𝒓_j)`$.

### 2.5 Power spectra

The power spectrum is the Fourier Transform of the correlation function. Since the correlation function involves a number of convolution-type integrals, the power spectrum is given by simple multiplications of the various density profile factors. Let

$$U(k)=\frac{𝑑r\mathrm{\hspace{0.17em}4}\pi r^2F(r)\mathrm{sin}(kr)/kr}{(2\pi )^3𝑑r\mathrm{\hspace{0.17em}4}\pi r^2F(r)}$$
(20)
denote the Fourier transform of the density run of the smooth component $`F`$ normalized by the total mass contained in the profile, and define $`u_i(k)`$ and $`U_c(k)`$ similarly (we have assumed we are working in three-dimensions). Then

$`P(k)`$ $`=`$ $`{\displaystyle \frac{dr\mathrm{\hspace{0.17em}4}\pi r^2}{(2\pi )^3}\frac{\mathrm{sin}(kr)}{kr}\xi (r)}`$ (21)
$`=`$ $`\overline{n}\left({\displaystyle \frac{M}{\overline{\rho }}}\right)^2[\left({\displaystyle \frac{M_s}{M}}\right)^2U(k)^2`$
$`+2{\displaystyle \underset{i}{}}{\displaystyle \frac{m_iM_s}{M^2}}u_i(k)U(k)U_c(k)`$
$`+{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}{\displaystyle \frac{m_im_j}{M^2}}u_i(k)u_j(k)U_c(k)^2`$
$`+{\displaystyle \underset{i}{}}\left({\displaystyle \frac{m_i}{M}}\right)^2u_i(k)^2],`$
which we could write as

$$P(k)P_{ss}(k)+2\underset{i}{}P_{si}(k)+\underset{i}{}\underset{ji}{}P_{ij}(k)+\underset{i}{}P_{ii}(k).$$

A similar calculation shows that the cross spectrum between subclumps and mass is

$`P_\times (k)`$ $`=`$ $`{\displaystyle \underset{i}{}}[{\displaystyle \frac{M_s}{\overline{\rho }}}U(k)U_c(k)`$ (22)
$`+U_c(k)^2{\displaystyle \underset{ji}{}}{\displaystyle \frac{m_j}{\overline{\rho }}}u_j(k)+{\displaystyle \frac{m_i}{\overline{\rho }}}u_i(k)],`$
and the power spectrum of the subclumps is

$$P_{cc}(k)=\overline{n}\left(\frac{M}{\overline{\rho }}\right)^2U_c(k)^2\underset{i}{}\underset{ji}{}\frac{m_i}{M}\frac{m_j}{M}.$$
(23)
If we set $`U_c(k)=U(k)`$, and further assume that the subclumps are much smaller than the parent halo, then the power spectrum reduces to the sum of the power spectra of the individual components:

$`{\displaystyle \frac{P(k)}{\overline{n}}}`$ $``$ $`\left({\displaystyle \frac{M}{\overline{\rho }}}\right)^2\left[1{\displaystyle \underset{i}{}}\left({\displaystyle \frac{m_i}{M}}\right)^2\right]U(k)^2`$ (24)
$`+{\displaystyle \underset{i}{}}\left({\displaystyle \frac{m_i}{\overline{\rho }}}\right)^2u_i(k)^2.`$
This shows that whether or not the substructure component dominates the small scale power depends on the fraction of mass in the subclumps, and on how much more dense the subclumps are compared to the smooth component.

To see this a little more clearly, consider a specific example. Suppose that the density runs are Gaussians, and that all the $`N`$ subclumps have the same mass $`m_i`$ and size $`R_i`$. If $`R`$ denotes the characteristic scale of the density runs of the smooth component, then $`U=\mathrm{exp}(k^2R^2/2)`$ and $`u_i=\mathrm{exp}(k^2R_i^2/2)`$. The subclumps dominate the power if $`N(m_i/\overline{\rho })^2u_i^2>(M/\overline{\rho })^2U(k)^2`$. This happens on scales where $`k^2>[2\mathrm{ln}(M/M_{cl})+\mathrm{ln}N]/R^2/[1(R_i/R)^2]`$, where we have set $`M_{cl}Nm_i`$. Note that if the total mass in subclumps is small, $`M_{cl}/M1`$, then the subclumps dominate the power only at very large $`k`$; if this mass is divided up among many subclumps, the subclumps only dominate at even higher $`k`$. For a fixed mass ratio, the exact scale on which the subclumps dominate depends on the size ratio $`R_i/R`$. This suggests that a feature in $`P(k)`$ at large $`k`$ may provide information about the nature of the subclumps. Simulations suggest that $`M/M_{cl}10`$ and $`R_i\stackrel{<}{}R/10`$, which yields a critical value of $`k\stackrel{>}{}\sqrt{6+\mathrm{ln}N}/R`$. It is interesting that this is just beyond the reach of current simulations.

Fig. 1 illustrates this with slightly more realistic Navarro, Frenk & White (1997) density profiles. The various curves show $`\mathrm{\Delta }(k)4\pi k^3P(k)`$ for models in which the density field is made up of Poisson distributed NFW halos. The parent halos are truncated at their virial radii $`r_{\mathrm{vir}}`$ defined so that the average density within $`r_{\mathrm{vir}}`$ is 200 times that of the background. We set the NFW core radius to be $`a_{\mathrm{NFW}}=0.1r_{\mathrm{vir}}`$. The dashed line (same in both panels) shows $`\mathrm{\Delta }(k)`$ if there is no substructure. We then assumed that a fraction $`f`$ of the mass of each parent halo was in subclumps, each of mass $`m`$. For simplicity, we assumed that the distribution of subclumps around the halo center was given by the same NFW form, and that the distribution of particles around each subclump center was also NFW, with core radius $`a_{\mathrm{NFW}}=0.1r_{\mathrm{sc}}`$. The value of $`r_{\mathrm{sc}}`$ was set by requiring that the average density within the subclumps equal $`\delta _{\mathrm{sc}}`$ times the background density. The solid curves in the panel on the left show equation (21) when $`f=0.1`$, $`m=0.01M`$ and $`\delta _{\mathrm{sc}}=500`$ (lower) and $`\delta _{\mathrm{sc}}=5000`$ (upper). This shows that, all other things being equal, denser subclumps contribute more power. The solid curves in the panel on the right show results when $`f=0.9`$, $`\delta _{\mathrm{sc}}=500`$ and $`m/M=0.1`$ and $`m/M=0.01`$ (departure from dashed curve apparent at high and still higher $`k`$). Comparison with the panel on the left shows that increasing the fraction of mass in subclumps, $`f`$, increases the small scale power. The scale on which this increase becomes apparent depends on $`(m/M)^{1/3}(\delta _{\mathrm{vir}}/\delta _{\mathrm{sc}})^{1/3}`$, the typical radii of the subclumps.

## 3 Higher-order statistics

Higher order correlations at a given small scale are dominated more strongly by the one-halo term than is the two point function, so the effect of substructure on the one-halo term is of great interest for higher order correlations. In this section we will consider the effect of substructure on the 3-point correlation function. The extension to higher orders is obvious.

We begin with the expansion of the density $`\rho (r)`$ given in equation (2). The 3-point function in real space is then defined by

$$\xi _{123}\delta (𝒓_1)\delta (𝒓_2)\delta (𝒓_3)\frac{1}{\overline{\rho }^3}\rho (𝒓_1)\rho (𝒓_2)\rho (𝒓_3).$$
(25)
It is convenient to consider the expressions for higher order correlations in Fourier space, since we will find that in the limiting case of interest the only additional term due to substructure is a term involving products of the substructure profile.

The Fourier transform of the 3-point function is the bispectrum $`B`$ defined by

$$<\delta (𝒌_1)\delta (𝒌_2)\delta (𝒌_3)>=B(𝒌_1,𝒌_2,𝒌_3)(2\pi )^3\delta _D(𝒌_{123}),$$
(26)
where the Delta function indicates that $`𝒌_1+𝒌_2+𝒌_3=0`$.

It is the sum of contributions from triplets which are all in the smooth component, triplets all in the same halo, plus contributions from the various cross terms. Writing all the different terms explicitly gives

$`B(k_1,k_2,k_3)=\overline{n}[\left({\displaystyle \frac{M_s}{\overline{\rho }}}\right)^3U(k_1)U(k_2)U(k_3)`$
$`+{\displaystyle \underset{i}{}}{\displaystyle \frac{m_i}{\overline{\rho }}}u_i(k_1)\left({\displaystyle \frac{M_s}{\overline{\rho }}}\right)^2U(k_2)U(k_3)U_c(k_1)`$
$`+{\displaystyle \underset{i}{}}{\displaystyle \frac{M_s}{\overline{\rho }}}U(k_1)\left({\displaystyle \frac{m_i}{\overline{\rho }}}\right)^2u_i(k_2)u_i(k_3)U_c(k_1)`$
$`+{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}{\displaystyle \frac{m_im_j}{\overline{\rho }^2}}u_i(k_1)u_j(k_2){\displaystyle \frac{M_s}{\overline{\rho }}}U(k_3)U_c(k_1)U_c(k_2)`$
$`+{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}{\displaystyle \underset{kji}{}}{\displaystyle \frac{m_im_jm_k}{\overline{\rho }^3}}u_i(k_1)u_j(k_2)u_k(k_3)`$
$`\times U_c(k_1)U_c(k_2)U_c(k_3)`$
$`+{\displaystyle \underset{i}{}}{\displaystyle \underset{ji}{}}{\displaystyle \frac{m_i}{\overline{\rho }}}u_i(k_1)\left({\displaystyle \frac{m_j}{\overline{\rho }}}\right)^2u_j(k_2)u_j(k_3)U_c(k_1)^2`$
$`+{\displaystyle \underset{i}{}}\left({\displaystyle \frac{m_i}{\overline{\rho }}}\right)^3u_i(k_1)u_i(k_2)u_i(k_3)],`$ (27)
where $`𝒌_3=𝒌_1k_2`$. Further integrations over appropriate window functions yield the spatially smoothed skewness (see, e.g., Scoccimarro et al. 2001; Takada & Jain 2002; Cooray & Sheth 2002). Similar relations hold for the higher order correlations.

In the same small mass and size limits which we used when studying $`\xi (r)`$ and $`P(k)`$, the bispectrum is dominated by the first and last terms of equation (27).

## 4 Averaging over numbers and masses

So far we have assumed that the masses and numbers of subclumps were fixed. This section shows the result of allowing for a distribution of masses and numbers. Expressions for the power spectra are much simpler than for $`\xi (r)`$, so we will only consider $`P(k)`$ here. Analogously, the averaged bispectrum $`B`$ is much simpler to evaluate than $`\xi _{123}`$. For what is to follow, it is useful to write explicitly that the quantity in the previous section is computed at fixed values of $`M`$, $`𝒩`$ and $`𝒎\{m_1,\mathrm{},m_𝒩\}`$, where the masses and density profiles of the $`𝒩`$ subclumps may be different. The quantity we are after is

$$𝑑M\frac{dN(M)}{dM}p(𝒩|M)p(𝒎|𝒩,M)P(k|𝒎,𝒩,M),$$
(28)
where $`dN/dM`$ is the number density of parent $`M`$ halos (usually called the universal halo mass function), $`p(𝒩|M)`$ is the probability an $`M`$ halo has $`𝒩`$ subclumps, $`p(𝒎|𝒩,M)`$ is the probability that the $`𝒩`$ subclumps were $`\{m_1,\mathrm{},m_𝒩\}`$ given that there were $`𝒩`$ of them in the $`M`$ halo, $`P(k|𝒎,𝒩,M)`$ is the quantity we computed in the previous section, and the sum is over all values of $`𝒩`$ and $`𝒎`$. To proceed, we need models of these different distributions.

If all the subclumps are identical, and there are $`𝒩`$ of them, then each of the terms in the sums over $`i`$ and $`j`$ are the same, so the contribution to the power from the smooth–subclump cross-terms is $`𝒩`$ times the contribution from a single clump, and the contribution from the subclump–subclump terms is $`𝒩(𝒩1)U_c^2`$. In this case, the total power depends on the first two factorial moments of the $`p(𝒩|M)`$ distribution. If the subclumps are much smaller than their parents, then equation (24) implies that

$`P(k)`$ $``$ $`\overline{n}\left({\displaystyle \frac{M}{\overline{\rho }}}\right)^2\left[1𝒩(m/M)^2\right]U(k)^2`$ (29)
$`+𝒩\overline{n}\left({\displaystyle \frac{m}{\overline{\rho }}}\right)^2u(k)^2.`$
The factor $`𝒩\overline{n}`$ is simply the number density of the subclumps, so it may be useful to think of the two terms above separately. The first term requires knowledge of the mass function of the parent halos, whereas the second requires the ‘mass function’ of the subclumps. As discussed in Section 5, accurate formulae for the parent mass function are available; subclump mass functions are only just becoming available.

In practice, halos of the same $`M`$ have different substructure distributions. For example, we expect that there will be some scatter around the mean number of subclumps $`N_{cl}𝒩|M`$ in an $`M`$halo, as well as in the fraction of mass associated with the clumped component: $`M_{cl}/M`$. However, if the scatter around the typical subclump configuration is small, then we should get a reasonable estimate of the effects of substructure if we use a good model of the typical configuration, and neglect the fact that there is actually some scatter around it. This is our strategy.

Let $`dn(m|M)/dm`$ denote the typical subclump mass function. Then

$`P(k)`$ $`=`$ $`{\displaystyle 𝑑M\frac{dN(M)}{dM}\left(\frac{M_s}{\rho }\right)^2U(k|M)^2}`$ (30)
$`+`$ $`2{\displaystyle 𝑑M\frac{dN(M)}{dM}\frac{M_s}{\rho }U(k)}`$
$`\times {\displaystyle }dm{\displaystyle \frac{dn(m|M)}{dm}}{\displaystyle \frac{m}{\rho }}u(k|m)U_c(k|m)`$
$`+`$ $`{\displaystyle 𝑑M\frac{dN(M)}{dM}𝑑m_1\frac{dn(m_1|M)}{dm_1}\frac{m_1}{\rho }u(k|m_1,M)}`$
$`\times {\displaystyle }dm_2{\displaystyle \frac{dn(m_2|M,m_1)}{dm_2}}{\displaystyle \frac{m_2}{\rho }}`$
$`\times u(k|m_2,M,m_1)U_c(k|m_1,m_2)^2`$
$`+`$ $`{\displaystyle 𝑑M\frac{dN(M)}{dM}}`$
$`\times dm{\displaystyle \frac{dn(m|M)}{dm}}{\displaystyle \frac{m^2}{\rho ^2}}u(k|m,M)^2.`$
The third term is the tricky one: it requires the joint distribution of $`m_1`$\- and $`m_2`$-subclumps in an $`M`$halo. Following Sheth & Lemson (1999), we have chosen to write this as the typical number of $`m_1`$-subclumps in an $`M`$halo times the typical number of $`m_2`$-subclumps in $`M`$-halos which have an $`m_1`$-subclump. Similarly, the notation $`U_c(k|m_1,m_2)`$ is intended to allow for the possibility that the distribution of $`m_2`$-subclumps within the parent $`M`$-halo may depend on whether or not the parent contains an $`m_1`$-subclump.

In the limit in which the subclumps are much smaller than the parent, and contribute a small fraction of the total mass, the expression above becomes

$`P(k)`$ $``$ $`{\displaystyle 𝑑M\frac{dN(M)}{dM}\left(\frac{M_s}{\rho }\right)^2U(k|M)^2}`$
$`+`$ $`{\displaystyle 𝑑M\frac{dN(M)}{dM}𝑑m\frac{dn(m|M)}{dm}\frac{m^2}{\rho ^2}u(k|m,M)^2}.`$
The integral over $`m`$ is from zero to $`M`$. If $`u(k|m,M)`$ is only a function of $`m`$ (i.e., if the density profile of an $`m`$ subclump is independent of the mass of the parent in which it sits), then the order of the integrals above can be rearranged. If we define the mass function of the subclumps as

$$\frac{dN_c(m)}{dm}_m^{\mathrm{}}𝑑M\frac{dN(M)}{dM}\frac{dn(m|M)}{dm},$$
(31)
then

$`P(k)`$ $``$ $`{\displaystyle 𝑑M\frac{dN(M)}{dM}\left(\frac{M_s}{\rho }\right)^2U(k|M)^2}`$ (32)
$`+{\displaystyle 𝑑m\frac{dN_c(m)}{dm}\left(\frac{m}{\rho }\right)^2u(k|m)^2}.`$
In this limit, the total power is well approximated by adding to the power associated with the smooth parent halos the contribution which comes integrating over the subclump mass function.

## 5 Details

The results presented above are general. When used to model large-scale structure simulations, one will almost always use models of the halo density profiles and halo mass functions which have been found to provide good descriptions of numerical simulations. Specifically, the density profiles $`F`$ and $`f`$ are usually described using the functional form given by Navarro, Frenk & White (1997),

$$\frac{\rho (r)}{\overline{\rho }}=\frac{\mathrm{\Delta }_{\mathrm{vir}}}{3\mathrm{\Omega }}\frac{c^3f(c)}{(r/r_s)(1+r/r_s)^2}$$
(33)
where the profile is truncated at the virial radius $`r_{\mathrm{vir}}`$, $`cr_{\mathrm{vir}}/r_s`$ is called the concentration, and $`f(c)=1/[\mathrm{ln}(1+c)c/(1+c)]`$. The virial radius $`r_{\mathrm{vir}}`$ is defined by requiring that $`m=4\pi r_{\mathrm{vir}}^3\overline{\rho }\mathrm{\Delta }_{\mathrm{vir}}`$. For spatially flat universes with $`\mathrm{\Omega }_0=(1,0.3)`$ and $`\mathrm{\Lambda }=1\mathrm{\Omega }`$, $`\mathrm{\Delta }_{\mathrm{vir}}=(178,340)`$. The Fourier transform of the density run around such a halo of mass $`m`$ is

$`u(k|m)`$ $`=`$ $`f(c)[\mathrm{sin}\kappa (\mathrm{Si}[\kappa (1+c)]\mathrm{Si}(\kappa )){\displaystyle \frac{\mathrm{sin}(\kappa c)}{\kappa (1+c)}}`$ (34)
$`+\mathrm{cos}\kappa (\mathrm{Ci}[\kappa (1+c)]\mathrm{Ci}(\kappa ))]`$
(Scoccimarro et al. 2001), where $`\kappa kr_{\mathrm{vir}}/c`$, $`\mathrm{Si}(x)=_0^xdt\mathrm{sin}(t)/t`$ is the sine integral and $`\mathrm{Ci}(x)=_x^{\mathrm{}}dt\mathrm{cos}(t)/t`$ is the cosine integral function. The concentration parameter of the halos depends on halo mass; we use the parametrization of this dependence given by Bullock et al. (2001):

$$c(m)\frac{9}{1+z}\left(\frac{m}{m_{}(z)}\right)^{0.1}.$$
(35)
Tha parent halo mass function is well described by

$$\frac{M^2}{\rho }\frac{dN(M,z)}{dM}\frac{\mathrm{d}M}{M}=\nu f(\nu )\frac{\mathrm{d}\nu }{\nu },$$
(36)
where $`\nu \delta _{\mathrm{sc}}^2(z)/\sigma ^2(M)`$, $`\rho `$ is the background mass density, and

$$\nu f(\nu )=A(p)\left(1+(q\nu )^p\right)\left(\frac{q\nu }{2\pi }\right)^{1/2}\mathrm{exp}\left(\frac{q\nu }{2}\right),$$
(37)
with $`p0.3`$, $`A(p)=[1+2^p\mathrm{\Gamma }(1/2p)/\sqrt{\pi }]^10.3222`$, and $`q0.75`$ (Sheth & Tormen 1999). Here $`\delta _{\mathrm{sc}}(z)`$ is the critical density required for spherical collapse at $`z`$, extrapolated to the present time using linear theory, and

$$\sigma ^2(M)=\frac{4\pi }{(2\pi )^3}_0^{\mathrm{}}\frac{\mathrm{d}k}{k}k^3P_{\mathrm{Lin}}(k)W^2(kR_0),$$
(38)
where $`W(x)=(3/x^3)[\mathrm{sin}(x)x\mathrm{cos}(x)]`$ and $`R_0=(3M/4\pi \rho )^{1/3}`$. That is to say, $`\sigma (M)`$ is the rms value of the initial fluctuation field when it is smoothed with a tophat filter of comoving size $`R_0`$, extrapolated using linear theory to the present time. If $`p=1/2`$ and $`q=1`$, then $`dN/dM`$ is the same as the universal mass function first written down by Press & Schechter (1974). Note that the mass function is normalized so that

$$𝑑MM𝑑N(M)/𝑑M=\rho .$$
(39)
The subclump distribution is less well constrained. Therefore, we will consider two models. The first is motivated by the results of numerical simulations (Tormen et al. 1998; Ghigna et al. 1999) which suggest that, on average, a power law in mass is a reasonable model of the typical subclump mass function:

$$\frac{dn(m|M)}{dm}dm=N_0\left(\frac{M}{m}\right)^\mu \frac{dm}{m}\mathrm{where}\mu <1,$$
(40)
and $`N_0`$ is a normalization constant which is set by the fraction of mass in subclumps:

$$\frac{M_{cl}}{M}=_0^M𝑑m\frac{m}{M}\frac{dn(m|M)}{dm}=\frac{N_0}{1\mu }.$$

In this model, although the total number of subclumps may diverge, the mass contained in them does not. The mean square subclump mass, which is related to the total number of subclump pairs is

$$M^2_0^M𝑑m\left(\frac{m}{M}\right)^2\frac{dn(m|M)}{dm}=M^2\frac{N_0}{2\mu }=M_{cl}M\frac{1\mu }{2\mu }.$$

Simulations suggest $`\mu 0.9`$ and $`M_{cl}/M0.1`$. If $`P_{\mathrm{Lin}}(k)k^n`$ and we use the Press–Schechter form for $`dN/dM`$, then in this power law subclump model, the subclump mass function is

$`{\displaystyle \frac{m^2}{\rho }}{\displaystyle \frac{dN_c(m)}{dm}}`$ $`=`$ $`{\displaystyle \frac{2^{(\mu 1)/\alpha }}{\mathrm{\Gamma }(1/2)}}\left({\displaystyle \frac{M_{}}{m}}\right)^{\mu 1}`$ (41)
$`\times \mathrm{\Gamma }({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\mu 1}{\alpha }},{\displaystyle \frac{(m/M_{})^\alpha }{2}}),`$
where $`\alpha =(n+3)/3`$ (the value $`n3/2`$ is a reasonable approximation to most CDM models). At small masses, $`m/M_{}1`$, the subclump mass function is a power law: $`dN_c/dmm^{\mu 1}`$, whereas at $`m/M_{}1`$ it drops exponentially.

In practice, the simulations only have finite mass resolution, so a power law model with a minimum mass cutoff may be more appropriate. A low-mass cutoff may also be a useful model if one is only interested in subclumps which could host galaxies. In this case, the mass fraction in subclumps is

$$\frac{M_{cl}}{M}=_{M_{min}}^M𝑑m\frac{m}{M}\frac{dn(m|M)}{dm}=N_0\frac{1(M_{min}/M)^{1\mu }}{1\mu }.$$

The lower mass cutoff means that the number of subclumps no longer diverges:

$$N_{cl}=_{M_{min}}^M𝑑m\frac{dn(m|M)}{dm}=N_0\frac{(M/M_{min})^\mu 1}{\mu }.$$

Therefore, when $`M_{min}M`$, the first term is the dominant one, and the number of subclumps increases as a power law in $`M`$. In practice, $`\mu `$ is sufficiently close to unity that the mass fraction in subclumps is approximately $`N_0/(1\mu )`$ only for the most massive halos.

Our second model of $`dn(m|M)/dm`$ is to assume that it is given by the progenitor distribution at some higher redshift $`z_1`$, and all progenitors lost the same fraction of mass as they fell into the parent they now occupy, so they are all a fraction $`f`$ of the mass they were at $`z_1`$. In this case,

$$\frac{dN_c(m,z_0)}{dm}=_{m/f}^{\mathrm{}}𝑑M\frac{dN(M,z_0)}{dM}\frac{dn(m/f,z_1|M,z_0)}{dm}.$$
(42)
Mass conservation means that this integral reduces to the mass function at the higher redshift: $`dN_c(m,z_0)/dm=dN(m/f,z_1)/dm`$. The total amount of mass in these subclumps is

$$𝑑mm𝑑N_c(m)/𝑑m=f𝑑m(m/f)𝑑N(m/f)/𝑑m=f\rho .$$

This provides a convenient way of thinking about the subclump mass function. In practice, it is likely that the subclumps which have survived to the present come from a distribution of earlier epochs, rather than a single epoch $`z_1`$, and the distribution of $`z_1`$ might peak at different redshifts for different values of $`m/M`$. Ignoring this subtlety, which is what the model above does, should be a reasonable approximation if the distribution around the mean $`z_1`$ is not broad. In practice, for small values of $`m/M`$, the subclump mass function is reasonably well approximated by a power law, and so this model also gives a subclump mass function which is in reasonable agreement with the simulations.

The virtue of this model is that it allows a simple estimate of the shape of the correlation function in the limit in which the subclumps are each a small fraction of the mass of the parent halo, and the total mass in subclumps is a small fraction of the total mass. Recall that, in this limit, the total correlation function is well approximated by the sum of two terms. The first term is the contribution from pairs which are from the smoothly distributed component, and the second is from the subclumps. The subclump contribution, then, can be determined by rescaling the contribution from the smooth component at $`z_1`$. If the particles which were stripped from the halos present at $`z_1`$ were random particles, then the number of pairs is lower by a factor of $`f^2`$; since the halos at $`z_1`$ were a factor of $`(1+z_1)^3`$ denser, the contribution to the correlation function is different by a factor of $`f^2(1+z_1)^3`$. This yields a factor which is probably slightly smaller than unity. On the other hand, if the mass was stripped from the earlier halos in shells, much like layers off an onion, as simulations suggest, then a better model of the subclump contribution is to truncate the halos which were present at $`z_1`$ at a fraction $`f`$ of their virial radii when estimating how the number of pairs changes with scale. Although this changes the actual shape of the correlation function (e.g., the subclump pairs are shifted to scales which are a factor of $`f`$ smaller), the typical factor by which the correlations are affected is $`f^1(1+z_1)^3`$, which can be considerably larger than unity.

## 6 Discussion

We have shown how to incorporate the effects of substructure into the halo model description of the nonlinear density field. Accounting for this substructure is important on scales smaller than the virial radii of typical halos. The effects are more pronounced for statistics which treat the subclumps preferentially, such as the power spectrum measured in studies of weak galaxy–galaxy gravitational lensing. Substructure will also change the dynamics within halos. Although we have not done so here, it is straightforward to insert our model for substructure into the halo model of the cosmic virial theorem, and the mean pairwise velocity and velocity dispersion developed in Sheth et al. (2001).

The stable clustering limit is a physically appealing description of clustering on small scales (Peebles 1980). It has been argued that a model with smooth halos is inconsistent with this limit (Ma & Fry 2000; Scoccimarro et al. 2001). Substructure changes the shape of the small scale power spectrum (c.f., Fig. 1); at least in principle, it can bring the halo model predictions into agreement with the stable clustering solution. However, it is not obvious that stable clustering is, indeed, the correct physical limit. Smith et al. (2002) argue that the stable clustering assumption is inconsistent with the results of high resolution numerical simulations. They also find that the simulations do not follow the small scale scaling predicted by models in which halos are smooth. Once an accurate model of the subclump mass function is available, it will be interesting to compare the predictions of our description of substructure with their results.

Although we have focussed primarily on the implications of substructure for the halo model of nonlinear clustering, our results have a wide range of other applications. For example, excess power in the Fourier transforms of images of galaxies or distant clusters can be used to infer the existence of spiral arms or substructure. This is the subject of work in progress. Closely related is the question of what images of high redshift galaxies may look like. Observations through a filter which has a fixed wavelength range probe the emission from high redshift galaxies at shorter restframe wavelengths than for galaxies at low redshift. If obscuration by dust is not a problem, and the UV luminosity is dominated by patchy star forming regions, then the images of high redshift galaxies should show considerable substructure. Our results suggest that, in this case, the power spectrum obtained by Fourier transforming the image of a high redshift patch of sky should show an increase in small scale power.

In addition, although we have phrased the entire discussion of substrucutre in terms of spatial statistics, this is not really necessary. Large databases describing various observed characteristics of galaxies are now becoming available (e.g., the 2dFGRS and SDSS surveys). If some of $`n`$ observables are correlated with others, the data will not fill the full $`n`$dimensional space available: the data set itself can be thought of as being clumpy, and the various clumps in dataspace may themselves have substructure. The formalism developed here provides a way of discovering, quantifying and modelling such substructure.

We thank Masahiro Takada for helpful discussions. B.J. acknowledges financial support from a NASA-LTSA grant and a Keck foundation grant.