text:
# Stellar Collisions and Pulsar Planets

## 1. Introduction

The subject of planet formation resulting from stellar encounters actually has a long history.111 For a proper history of theories of solar system formation, see Brush (1990). For many years a ‘tidal theory’ of planet formation in our own solar system held sway (Chamberlin 1901; Moulton 1905; Jeans222I’ve been unable to verify this particular reference. 1919; Jeffreys 1929) These postulated that the planets arose as condensations from filaments stripped off the sun or a companion during what we would now call an exchange interaction (Lyttleton 1936)333 Lyttleton’s paper would have been perfect for this meeting!. Eventually this theory fell in the face of angular momentum problems (Russell 1935), the inability of such hot filaments to condense rapidly enough to form planets (Spitzer 1939) and the presence of Deuterium in the planets which should have been destroyed while part of the sun (Cameron 1965).

Nevertheless, recent years have given us cause to reinvestigate some of these phenomena in the specific context of the pulsar planets. The discovery of two earth mass planets orbiting the pulsar B1257+12 (Wolszczan & Frail 1992; Wolszczan 1994) began a period of frenzied speculation as to the origin of these bodies in such an unusual location (see Podsiadlowski 1993 for a review of the scenarios). Most scenarios proposed harken back to the spirit of the ‘tidal theory’ of planet formation, primarily because the evolutionary history of the central neutron star precludes a monistic origin for these particular planets.

The purpose of this paper is to review what we think we know about the pulsar planet origins and to place this in the broader context of this meeting, ending with some speculation about variants on this theme to be expected from the broad range of stellar interactions under discussion here. In section 2. I will review some of the theories regarding the kinds of incidents that led to the formation of the pulsar planets. In section 3. I will describe an attempt to unify some of these scenarios in terms of a single theory of disk evolution differing only in global parameters (essentially an extension of Phinney & Hansen 1993). Finally in section 4. I will discuss a few other possible planetary systems resulting from mergers.

## 2. Possible Origins for Pulsar Planets

To properly assess the viability of the various proposed scenarios for the pulsar planets, we must review several pieces of pertinent information. The two outermost planets (the third and smallest planet is unimportant in terms of mass and angular momentum) are located at 0.36 and 0.47 Au, with masses of 3.4 and 2.8 $`M_{}`$ respectively (the true values being larger by a factor of 1/$`\mathrm{sin}i`$, where $`i`$ is the inclination angle of the planet to the line of sight). This then requires a minimum of $`1.3\times 10^{48}\mathrm{ergs}\mathrm{s}`$ of angular momentum. If these earth mass planets are rocky like our own, then an origin in solar composition material would require a mass budget $`1/Z_{}`$ larger, $`10^3M_{}`$.

The observed eccentricities are small, in striking contrast to some of the more ‘conventional’ planetary systems recently discovered. This is a problem for those who propose that these planets are primordial (Bailes, Lyne & Shemar 1991), in the sense that they formed like any other planetary system and survived the subsequent evolution, including the supernova that resulted in the pulsar. The fact that the pulsar has a 300 km/s space velocity makes such a scenario even more unlikely.

Another illuminating fact is that this is a millisecond pulsar (spin period 6.2 ms). While most pulsars are ‘young’, in the sense that they have resulted from a supernova that occurred within the last $`10^610^8`$ years, millisecond pulsars are thought to be on average somewhat444 Although I will argue later that this is not the case for this particular system. older (see Phinney & Kulkarni 1994 for a review). This is believed to result from an extended period of mass transfer from a companion, in which the accreted angular momentum revives the old, spun-down pulsar in a new incarnation. Thus, a natural setting for the pulsar planets is origin in a circumpulsar disk which is also accreting onto the central object and converting it into a millisecond pulsar (Memnonides scenarios, in the parlance of Phinney & Hansen). There is no shortage of such scenarios and it is these on which we will now focus.

While the circumpulsar disk is a natural ingredient for both millisecond pulsar formation and planet formation, the source of this disk is absent. Thus, we must turn to some kind of catastrophic incident to provide the end result of a neutron star surrounded by a disk, yet with no remaining donor.

1. One possibility is simply that binary mass transfer proceeds as usual but the donor is eventually disrupted or evaporated. Stevens, Rees & Podsiadlowski (1992) propose such a scenario in which a low mass companion is initially evaporated by pulsar radiation (such as appears to be happening in systems like PSR 1957+20 and PSR 1744-24). If the mass loss is rapid enough, the star expands adiabatically and may eventually suffer dynamical disruption (if sufficient mass is lost from the system), providing the necessary disk of material.
2. A variant on this (Tavani & Brookshaw 1992; Banit, Ruderman & Shaham 1992) is that material evaporated from the companion may settle into a circumbinary disk. The companion may then be relentlessly evaporated into oblivion, leaving the remnant disk to form planets. This scenario and the former have the attractive feature that they take a known phenomenon to a logical conclusion. On the other hand, the evaporation of the companion takes $`>10^8`$ years, based on the observed systems. The high velocity and small distance above the Galactic plane suggest that the B1257+12 system is ‘young’ $`0.6\mathrm{kpc}/300\mathrm{k}\mathrm{m}.\mathrm{s}^12\times 10^6`$ years. Thus, the most natural interpretation is that the formation of the planets and the spin-up of the pulsar to a millisecond period occurred soon after the supernova kick. We will use this estimate hereafter rather than the pulsar spin-down age ($`10^9`$ years) as spin-down ages are notoriously unreliable for millisecond pulsars (e.g. Camilo, Thorsett & Kulkarni 1994; Hansen & Phinney 1998)
3. Another possibility is if the kick received by the neutron star happens to be directed at a close binary companion, so that the two suffer a direct collision. The tidal disruption of the companion results in the required disk. While this scenario suffers from the fine tuning required to engineer a collision, it is the only one which naturally gives rise to rapid planet formation soon after the supernova event, as required above.
4. A completely different possibility is one in which the pulsar planet system is most closely linked with double neutron star systems. In a binary containing a neutron star and a massive Be star, the neutron star may capture some material from the Be star wind, forming the required disk (Fabian & Podsiadlowski 1991). The Be star eventually explodes in a supernova, unbinding the binary and leaving the first neutron star with the required disk. This scenario has the advantage of providing a lot of angular momentum but there is no natural method of providing a large kick for the final system, as the second supernova kick is applied to the companion, not the millisecond pulsar.
5. Finally, we have the case of those close, compact object binary systems which undergo gravitational wave induced orbital shrinkage and merger. A merger between two white dwarfs could form a pulsar by accretion-induced collapse and the resulting remnant disk could form planets (Podsiadlowski, Pringle & Rees 1991). This scenario has the twin advantages of having a known source population and providing metal-rich material for the formation of planets. The main problem with this scenario is the uncertain physics surrounding the issue of accretion-induced collapse (e.g. Canal et al 1990; Nomoto & Kondo 1991).

All of these scenarios reach a common end-point, namely the formation of a gaseous disk orbiting the neutron star. It is from this disk that we expect the planets to form. We may thus consider the viability of the various models above based on the mass and angular momentum they provide to make planets in the appropriate locations.

## 3. Gas Disk Evolution

In the scenarios discussed above, disks of substantial mass ($`0.1M_{}`$) may form, but generally on scales $`10^910^{11}`$ cm (thereby providing angular momenta in the range $`10^{49}10^{52}`$ ergs s, depending on the scenario), around the neutron star. This is well within the location of the planets, so the planets are formed from the part of the disk that expands to conserve angular momentum while most of the material is accreted (Lynden-Bell & Pringle 1974). Furthermore, we assume these bodies of earth mass are rocky in nature, as even the gas giants in our own solar system require a rocky core of $`10M_{}`$ before they can capture gaseous material (e.g. Mizuno 1980). This means that only the metallicity fraction of the material is of interest to us, thereby increasing the amount of original mass required (except, perhaps, for the case of disrupted white dwarfs).

The expansion of such disks under the influence of simple viscosity prescriptions has been studied by Ruden (1993) and Phinney & Hansen (1993). While the expanding disk contains enough mass and angular momentum to make the planets, by the time the disk becomes cool enough to reasonably form dust and rocks, the material is spread over several Au, so that the planet formation has to be both efficient and able to accumulate mass over several Au into a few planets within 1 Au.

### 3.1. Constant $`\alpha `$ Models

To demonstrate this, consider the evolution of a disk surface density $`\mathrm{\Sigma }`$ under the influence of viscosity $`\nu `$, given by mass and angular momentum conservation,

$$\frac{\mathrm{\Sigma }}{t}=\frac{3}{R}\frac{}{R}\left(R^{1/2}\frac{}{R}\left(\nu \mathrm{\Sigma }R^{1/2}\right)\right).$$
(1)
The viscosity is parameterised by an adjustable constant $`\alpha `$ in the standard way (Shakura & Sunyaev 1973) such that $`\nu =\alpha c_s^2/\mathrm{\Omega }`$, where $`c_s`$ is the gas sound speed in the disk and $`\mathrm{\Omega }`$ is the orbital angular frequency. Thermal balance is enforced by assuming that the energy dissipated by $`\nu `$ is radiated from the surface of the disk and vertical energy transport is by radiative diffusion, modulated by an opacity table due to Ruden & Pollack (1991). Of particular interest to us are those low temperature regions where the material is cool enough for dust to form. Once dust forms and agglomerates it may settle to the midplane and eventually form planetesimals to begin the process of planetary accumulation. Thus, in this simple model we will judge the success of the scenario by how much mass is in cool regions where material could condense in this manner.

The simple $`\alpha =0.1`$ disk model yields uncomfortable results (Figure 1) in that the material spreads rapidly over several Au. Even if material could condense into solid bodies in the hotter phases, it would have to do so rapidly, as the viscous times for the disk to expand beyond 1 Au are short,

$$t_\nu \frac{R^2}{\nu }\frac{6}{\alpha }\mathrm{years}\left(\frac{T}{3000K}\right)^1\left(\frac{R}{1Au}\right)^{1/2}.$$
(2)
Another way to demonstrate this is to examine the temporal evolution of the surface density at particular disk locations (Figure 2). Thus, any reasonable timescale for planet formation results in material spread over many Au, much like our own solar system. However, the planets are all located within 1 Au. This problem persists for the full range of potential global parameters (Phinney & Hansen 1993). Lowering $`\alpha `$ can prolong the period of high surface density inside 1 Au, but still leads to a lot of mass at large radii.

### 3.2. Physically Motivated Models

However, there is no reason why $`\alpha `$ should be constant throughout the disk. This convenient assumption is usually simply the result of having no better, physically motivated model. Despite considerable work over the last few decades (see Papaloizou & Lin 1995 for a review), there is only one angular momentum transport mechanism that is generally regarded as robust, namely the magnetorotational instability (Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991). This mechanism works only as long as the disk material is sufficiently ionized to couple magnetic field and gas motions. Thus, in the outer regions of protoplanetary disks, the instability breaks down (Gammie 1996; Sano & Miyama 1999; Wardle 1999). If there is no other angular momentum transport mechanism at work, then this material may be quiescent, except perhaps for a thin surface layer which is maintained at a sufficient level of ionization by cosmic rays (Gammie 1996) or illumination from the central star (Glassgold, Najita & Igea 1997). Thus, as an alternative model, we may consider the expansion of a disk using the same model as before, but where we set $`\alpha =0`$ when $`T<3000`$ K. The end result is that the disk tends to an asymptotic quiescent disk, whose profile is dictated by the transition temperature and whose extent is determined by the requirements of global angular momentum conservation. The profile corresponding to this transition temperature is

$$\mathrm{\Sigma }=1.1\times 10^4\mathrm{g}.\mathrm{cm}^2\left(\frac{\alpha }{0.1}\right)^{1/2}\left(\frac{R}{1Au}\right)^{3/4}.$$
(3)
Integrating over this profile provides a relationship between the outer radius of quiescent disk and the total angular momentum stored

$$R_{out}=1.61Au\left(\frac{J_{tot}}{10^{51}\mathrm{ergs}\mathrm{s}}\right)^{4/13}.$$
(4)
The evolution to this state is shown in Figure 3. Requiring the total angular momentum in the gas disk to be $`1/Z_{}`$ larger than the total planetary angular momentum today, we find an outer radius of 0.7 Au. The total rocky mass in this asymptotic disk is $`7M_{}`$. Thus, this hypothesis provides an excellent order-of-magnitude agreement with the observations.

In reality, the accretion onto the central object will result in copious X-rays, and the surface layers of the disk may be ionized sufficiently that the magnetorotational instability may operate and some accretion may continue. Thus, we now consider a model like the previous one, but in the spirit of Gammie (1996), in which a critical surface density $`\mathrm{\Sigma }=10g/cm^2`$ is assumed to remain sufficiently ionized to continue to evolve with $`\alpha =0.1`$. For those radii where the total $`\mathrm{\Sigma }>10g/cm^2`$, we assume the remainder is quiescent as before. This changes the evolution somewhat, since now both the inner (by virtue of their higher temperature) and the outer (by virtue of their lower surface density) layers remain active and the dead zone resides at intermediate radii. It is particularly encouraging that this model produces a dense, quiescent disk at precisely the radii where the pulsar planets reside (Figure 4).

The modifications to the original calculations now provide us with a good model, in accordance with the broad requirements of the observed planets. We can accumulate sufficient mass and angular momentum in cool material in the region at $`1Au`$ and keep it there without spreading too thin. The viscous timescales $`10^310^4`$ years are also short enough to be consistent with the requirements of a young age for the system. Note however this refers only to the gaseous phase of the evolution. The timescale requirements also demand that the heavy elements condense into solid form, accumulate into planetesimals and thence into planets within $`10^610^7`$ years. The formation of planetesimals is unlikely to be a problem in this regard as timescales $`10^310^4`$ years are expected (Weidenschilling & Cuzzi 1993).

### 3.3. Accumulation Into Planets

Assuming the formation of planetesimals occurs on timescales similar to the viscous spreading time, we must then consider the accumulation of planetesimals into larger bodies. A rapid initial accumulation leads to planetary embryos, isolated from each other in the sense that each has swept up the surrounding mass in an annulus of width appropriate to the embryo’s “Hill radius”, $`\mathrm{\Delta }R\mu ^{1/3}R`$, where $`\mu `$ is the mass ratio with respect to the central object. After this, longer range gravitational perturbations bring the isolated embryos into crossing orbits again and the accumulation proceeds. By taking the surface density of the quiescent disk in Figure 4, we can calculate the expected embryo mass as a function of radius

$$M_p28.4M_{}\left(\frac{R}{1Au}\right)^3\left(\frac{Z}{Z_{}}\right)^{3/2}\left(\frac{\mathrm{\Sigma }}{10^5g/cm^2}\right)^{3/2},$$
(5)
which is shown in Figure 5.

We see that the innermost planet is indicative of the kind of building blocks one expects for the final process of planetary accumulation. The outer two planets are larger, suggesting that further accumulation has most likely taken place. It is interesting to speculate about whether this final stage of accumulation has, in fact, ended in this particular system, given that estimates for our own solar system indicate this stage lasts for $`10^8`$ years (Wetherill 1990), i.e. somewhat longer than the likely lifetime of this system. The final clearing stage takes place on a timescale determined by the time it takes the largest bodies to scatter the smaller bodies onto parabolic orbits (Tremaine 1993)

$$t1.8\times 10^9years\left(\frac{M_p}{M_{}}\right)^2\left(\frac{a_p}{1Au}\right)^{3/2}$$
(6)
The mass versus semi-major axis relation required to accomplish this feat in $`10^7`$ years is shown by the curve labelled ‘Evolution’ in Figure 5. We see that, while the known bodies are not large enough to accomplish this feat in the expected time allotted, the largest planet is within a factor of two of the criterion, so that significant dynamical evolution of the embryo swarm is expected. However, some relic component of the embryo swarm may yet remain for future detection. It should also be noted that these bodies are not massive enough to eject most planetesimals, rather they will accrete the smaller bodies (Tremaine 1993) as long as

$$M_p<21.5M_{}\left(\frac{a_p}{1Au}\right)^{3/2}\left(\frac{\rho }{3g/cm^3}\right)^{1/2}\left(\frac{\mathrm{\Delta }i}{0.1rad}\right)^{3/4},$$
(7)
where $`\rho `$ is the average planet density and $`\mathrm{\Delta }i`$ is the average inclination of the planetesimal bodies. This curve is also shown in Figure 5, labelled ‘Collision’. The planets lie well within this bound, so that little rocky material has likely been ejected from the system (and no Oort cloud is expected).

### 3.4. Further Planets In This System

The preceding discussion is based solely on the three planets reported by 1994. Further timing of this system has revealed long term timing structure that could indicate the presence of additional planets in this system (Wolszczan et al 2000). Figure 6 shows the inferred planet parameters assuming a single additional body in a circular orbit (Joshi & Rasio 1997). There is considerable uncertainty in the planet parameters, stemming largely from the uncertainties in how much of the contribution to the pulsar first frequency derivative is due to the additional orbital motion. For bodies on the lower end of this curve, little adjustment need be made in the above arguments. However, bodies in the upper end of the allowed mass and distance range could well invalidate much of the above description. The problem is that the angular momentum required to generate such an object is larger (once augmented by the lost gas fraction) than many scenarios permit.

If the fourth body does turn out to be $`80M_{}`$ at $`30Au`$, it could perhaps represent the last remnant of a donor (such as in scenarios 1 and 2), in which case the angular momentum is primordial and is not included in the discussions above. On the other hand, if the anomalous timing is due to one or more smaller mass bodies in eccentric orbits, then it would be consistent with the low mass of the planets (and their inability to successfully clear the environment of smaller bodies) and the short lifetime of the system.

## 4. Speculations About Other Possible Systems

I have concentrated in the previous section on the scenario in which a solar metallicity disk expands and forms planets. While many scenarios can be analysed in a unified form by simply varying the global mass and angular momentum budgets (e.g. Phinney & Hansen 1993), compositional variations are not so easily attempted. It is not sufficient to simply vary Z, the metal content, in converting gas to solid bodies. The contribution of heavy elements to the opacity and thus the thermal evolution of the disk is critical, especially at lower temperatures where the molecular opacities dominate and which are the regions that especially concern us. Unfortunately, few opacity tables exist for the more exotic mixtures and we must content ourselves with a few more general speculations.

This is unfortunate, because many of the most interesting events fall into this category. Perhaps the most likely candidates are planets formed in the disks resulting from the merger of two white dwarfs. This is a source population we expect to find and of particular interest in the case of type Ia supernovae. Even if the central object does not collapse to form a pulsar, we might expect planets around massive white dwarfs (Livio, Pringle & Saffer 1992). How might the planetary systems we anticipate differ from those discussed above? The first effect worth noting is the larger expected opacity, which will maintain the disk at a higher temperature. On the other hand, Carbon has a higher ionization potential than Hydrogen and so a disk composed primarily of Carbon and Oxygen will recombine at higher temperatures, shutting off the viscosity earlier. Finally, nearly all of the gaseous disk is potentially available to be converted into planets. To get a general idea of the expected system, I repeated the calculations above with a toy model where each opacity was simply increased by a factor 50 ($`=1/Z_{}`$) and quiescence was assumed to set in below $`10^4`$ K. Figure 7 shows the equivalent of Figure 5 in this case, i.e. the expected size and semi-major axis of bodies expected from a disk of heavy elements. The fact that we expect significant mass ejection also suggests that the ejecting bodies will migrate inwards (Murray et al 1998). Thus, we expect any planets formed in such a disk to be $`30300M_{}`$ and located within $`a<0.2`$ Au.

Another possibility involves a disk of even heavier elements surrounding a neutron star. This could arise as in the recent ‘Fallback disk’ model of anomalous X-ray pulsars (Chatterjee, Hernquist & Narayan 2000) or, more appropriate to our conference, a remnant disk of neutron-rich material resulting from the merger of two neutron stars. Such a disk is susceptible to explosion due to the heat input from radioactive decay (Li & Paczynski 1998). However, given the deep potential well of the merger remnant, material inside a radius $`R_{bound}ϵ^1GM/c^2`$ (where $`ϵ`$ is the fraction of the nuclear rest mass energy released in the radioactive decay) will remain bound. For $`ϵ10^3`$, this is $`4\times 10^8`$ cm, yielding a significant disk with the concomitant possibility of planets. The thermal balance of this disk will be complicated not only by the unusual opacities, but also the additional heat source resulting from the radioactivity. Planets formed from such material may themselves be directly observable due to heating from radioactivity. While estimates are obviously wildly uncertain, order of magnitude numbers are

$$L0.5L_{}\frac{ϵ}{10^3}\frac{M_p}{10M_{}}\left(\frac{\tau }{10^6years}\right)^1$$
(8)
where $`ϵ`$ is the fraction of rest mass energy released in radioactivity on some characteristic timescale $`\tau `$. Effective temperatures in this case are $`40000`$ K, assuming a radius of $`10^9`$ cm. The observability will die away on the timescale $`\tau `$ as the heat capacity of such objects is small. This signature is likely to be quite unique; a body with an extremely strange spectrum consisting entirely of heavy elements, in close orbit about a dark compact object. Note that, if such planets do form, they do so only once the disk has expanded to radii $`>10^{11}`$ cm (the Roche limit for a rocky body around a 3 $`M_{}`$ compact object).

## 5. Conclusion

The models presented above are obviously quite simple, but it is encouraging that the physically most well-motivated models give numbers that are in good accord with the observations. However, this agreement is threatened by the appearance of further timing residuals in the data. If this turns out to be due to lower mass objects in eccentric orbits then it will represent a striking confirmation of the picture presented here. However, if it is due to a large object at Neptunian distances, then it could invalidate the preceding discussion.

Extrapolating our models to the systems of planets that could result from compact object mergers, we are limited by our lack of knowledge of the microphysics of gas of such unusual composition. Nevertheless, we find that quite massive planets in close orbits are possible. These could be found orbiting massive white dwarfs, neutron stars or black holes.

Support for this work was provided by NASA through Hubble Fellowship grant #HF-01120.01-99A, from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

## References

Bailes, M., Lyne, A. G. & Shemar, S. L., 1991, Nature, 352, 311

Balbus, S. A. & Hawley, J. F., 1991, ApJ, 376, 214

Banit, M., Ruderman, M. & Shaham, J., 1993, ApJ, 415, 779

Brush, S. G., 1990, Rev. Mod. Phys. 62, 43

Cameron, A. G. W., 1965, in “Introduction to Space Science”, ed. W.N.Hess, (Gordon & Breach, New York), p 553

Camilo, F., Thorsett, S. E. & Kulkarni, S. R., 1994, ApJ, 421, L15

Canal, R., Garcia, D., Isern, J. & Labay, J., 1990, ApJ, 356, L51

Chandrasekhar, S., 1960, Proc. Nat. Acad. Sci., 46, 253

Chamberlin, T. C., 1901, ApJ, 14, 17

Chatterjee, P., Hernquist, L. & Narayan, R., 2000, ApJ, 534, 373

Fabian, A. C., & Podsiadlowski, P., 1991, Nature, 353, 801

Gammie, C. F., 1996, ApJ, 457, 355

Glassgold, A. E., Najita, J. & Igea, J., 1997, ApJ, 480, 344

Hansen, B. M. S. & Phinney, E. S., 1998, MNRAS, 294, 569

Jeans, J. H., 1919, “Problems of Cosmogony & Stellar Dynamics”, Adams Prize Essay for 1917, Cambridge University Press, Cambridge

Joshi, K. J. & Rasio, F. A., 1997, ApJ, 479, 948

Jeffreys, H., 1929, “The Earth”, Cambridge Univ. Press, Cambridge, p27

Li, L.-X., & Paczynski, B., 1998, ApJ, 507, L59

Livio, M., Pringle, J. E. & Saffer, R. A., 1992, MNRAS, 257, P15

Lynden-Bell, D. & Pringle, J. E., 1974, MNRAS, 168, 603

Lyttleton, R. A., 1936, MNRAS, 96, 559

Mizuno, H., 1980, Prog. Theor. Phys. 64, 544

Moulton, F. R., 1905, ApJ, 22, 165

Murray, N., Hansen, B., Holman, M. & Tremaine, S., 1998, Science, 279, 69

Nomoto, K. & Kondo, Y., 1991, ApJ, 367, L19

Papaloizou, J. C. B. & Lin, D. N. C., 1995, ARA&A, 33, 505

Phinney, E. S. & Hansen, B. M. S., 1993, in “Planets around Pulsars”, ASP Conf. Ser. 36, ed. J.A.Phillips, S. E. Thorsett & S. R. Kulkarni, p371

Phinney, E. S. & Kulkarni, S. R., 1994, ARA&A, 32, 591

Podsiadlowski, P., 1993, in “Planets around Pulsars”, ASP Conf. Ser. 36, ed. J.A.Phillips, S. E. Thorsett & S. R. Kulkarni, p149

Podsiadlowski, P., Pringle, J. E. & Rees, M. J., 1991, Nature, 352, 783

Ruden, S. P., 1993, in “Planets around Pulsars”, ASP Conf. Ser. 36, ed. J.A.Phillips, S. E. Thorsett & S. R. Kulkarni, p197

Ruden, S. P. & Pollack, J. B., 1991, ApJ, 375, 740

Russell, H. N., 1935, “The Solar System & it’s Origins”, New York: McMillan, p112

Sano, T. & Miyama, S. M., 1999, ApJ, 515, 776

Shakura, N. I. & Sunyaev, R. A., 1973, A&A, 24, 337

Spitzer, L., 1939, ApJ, 90, 675

Stevens, I. R., Rees, M. J. & Podsiadlowski, P., 1992, MNRAS, 254, 19

Tavani, M. & Brookshaw, L., 1992, Nature, 356, 320

Tremaine, S., 1993, in “Planets around Pulsars”, ASP Conf. Ser. 36, ed. J.A.Phillips, S. E. Thorsett & S. R. Kulkarni, p335

Velikhov, E. P., 1959, J. Exp. Theor. Phys. (U.S.S.R.) 36, 1398

Wardle, M., 1999, MNRAS, 307, 849

Weidenschilling, S. J. & Cuzzi, J. N., 1993, in “Protostars and Planets III”, ed. E.H.Levy, & J.I.Lunine, (Tucson: University of Arizona Press),p 1031

Wetherill, G. W., 1990, Ann. Rev. Earth Planet. Sci., 18, 205

Wolszczan, A., Doroshenko, O., Konacki, M., Kramer, M., Jessner, A., Wielebinski, R., Camilo, F., Nice, D. J. & Taylor, J. H., 2000, ApJ, 528, 907

Wolszczan, A. & Frail, D. A., 1992, Nature, 355, 145

Wolszczan, A., 1994, Science, 264, 538


format: text

text:
# Quasi-sparse eigenvector diagonalization and stochastic error correction

## 1 Quasi-sparse eigenvector method

Quasi-sparse eigenvector (QSE) diagonalization is a new computational method which finds approximate low-lying eigenvalues and eigenvectors for a general quantum field Hamiltonian $`H`$ . It handles the exponential increase in the dimension of Fock space by exploiting the sparsity of the Hamiltonian. The method is most effective when the splitting between low-lying eigenvalues is not too small compared to the size of the Hamiltonian off-diagonal entries. In such cases the low-lying eigenvectors are quasi-sparse, meaning that the vector is dominated by a small fraction of its largest components. The QSE algorithm is applied using the following steps:

1. Select a subset of orthonormal basis vectors $`\{e_{i_1},\mathrm{},e_{i_n}\}`$ and call the corresponding subspace $`S`$.
2. Diagonalize $`H`$ restricted to $`S`$ and find one eigenvector $`v.`$
3. Sort the basis components $`e_{i_j}|v`$ according to their magnitude and remove the least important basis vectors.
4. Replace the discarded basis vectors by new basis vectors. These are selected at random from a pool of candidate basis vectors which are connected to the old basis vectors through non-vanishing matrix elements of $`H`$.
5. Redefine $`S`$ as the subspace spanned by the updated set of basis vectors and repeat steps 2 through 5.

If the subset of basis vectors is sufficiently large, the exact low energy eigenvectors will be stable fixed points of the QSE update process. We can show this as follows. Let $`|i`$ be the eigenvectors of the submatrix of $`H`$ restricted to the subspace $`S`$, where $`S`$ is the span of the subset of basis vectors after step 3 of the QSE algorithm. Let $`|A_j`$ be the remaining orthonormal basis vectors in the full space not contained in $`S`$. We can represent $`H`$ as

$$\left[\begin{array}{cccccc}\lambda _1& 0& \mathrm{}& 1\left|H\right|A_1& 1\left|H\right|A_2& \mathrm{}\\ 0& \lambda _2& \mathrm{}& 2\left|H\right|A_1& 2\left|H\right|A_2& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ A_1\left|H\right|1& A_1\left|H\right|2& \mathrm{}& E\lambda _{A_1}& A_1\left|H\right|A_2& \mathrm{}\\ A_2\left|H\right|1& A_2\left|H\right|2& \mathrm{}& A_2\left|H\right|A_1& E\lambda _{A_2}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right].$$
(1)
We have written the diagonal terms for the basis vectors $`|A_j`$ with an explicit factor $`E`$ for reasons to be explained shortly. We let $`|1`$ be the approximate eigenvector of interest and shift the diagonal entries so that $`\lambda _1=0.`$ Our starting hypothesis is that $`|1`$ is close to some exact eigenvector of $`H`$ which we denote as $`|1_{\text{full}}`$. More precisely we assume that the components of $`|1_{\text{full}}`$ outside $`S`$ are small enough so that we can expand in inverse powers of the introduced parameter $`E.`$

We now expand the eigenvector as

$$|1_{\text{full}}=\left[\begin{array}{c}1+\mathrm{}\\ c_2^{}E^1+\mathrm{}\\ \mathrm{}\\ c_{A_1}^{}E^1+\mathrm{}\\ c_{A_2}^{}E^1+\mathrm{}\\ \mathrm{}\end{array}\right]$$
(2)
and the corresponding eigenvalue as

$$\lambda _{\text{full}}=\lambda _1^{}E^1+\mathrm{}.$$
(3)
In (2) we have constrained $`|1_{\text{full}}`$ to have unit norm. From the eigenvalue equation $`H|1_{\text{full}}=\lambda _{\text{full}}|1_{\text{full}}`$ we find at lowest order

$$c_{A_j}^{}=\frac{A_j\left|H\right|1}{\lambda _{A_j}}.$$
(4)
We see that at lowest order the component of $`|1_{\text{full}}`$ in the $`|A_j`$ direction is independent of the other vectors $`|A_j^{}`$. If $`|1`$ is sufficiently close to $`|1_{\text{full}}`$ then the limitation that only a fixed number of new basis vectors is added in step 4 of the QSE algorithm is not important. At lowest order in $`E^1`$ the comparison of basis components in step 3 (in the next iteration) is the same as if we had included all remaining vectors $`|A_j`$ at once. Therefore at each update only the truly largest components are kept and the algorithm will converge to some optimal approximation of $`|1_{\text{full}}`$. This is consistent with the actual performance of the algorithm for a wide range of examples.

## 2 Stochastic error correction

The quasi-sparse eigenvector method approximates the low energy eigenvalues and eigenvectors of a Hamiltonian $`H`$ by finding and diagonalizing optimized subspaces $`S`$. This procedure however has an intrinsic error due to the omission of basis states not in $`S`$. The goal of stochastic error correction is to sample the contribution of these remaining basis vectors in a systematic fashion. Due to space limitations we consider here only the series expansion method. Several other techniques are discussed in a forthcoming paper .

The series method for stochastic error correction is based on the $`E^1`$ expansion introduced in the previous section,

$$|1_{\text{full}}\left[\begin{array}{c}1\\ c_2^{}E^1+c_2^{\prime \prime }E^2+\mathrm{}\\ \mathrm{}\\ c_{A_1}^{}E^1+c_{A_1}^{\prime \prime }E^2+\mathrm{}\\ c_{A_2}^{}E^1+c_{A_2}^{\prime \prime }E^2+\mathrm{}\\ \mathrm{}\end{array}\right],$$
(5)
$$\lambda _{\text{full}}=\lambda _1^{}E^1+\lambda _1^{\prime \prime }E^2\mathrm{}.$$
(6)
This time it is more convenient to choose the normalization of the eigenvector such that the $`|1`$ component remains 1. The convergence of the expansion is controlled by the proximity of $`|1`$ to $`|1_{\text{full}}`$. If $`|1`$ is not at all close $`|1_{\text{full}}`$ then it will be necessary to use a more robust method, and such methods are discussed in .

At first order in $`E^1`$ we find

$$c_{A_j}^{}=\frac{A_j\left|H\right|1}{\lambda _{A_j}}$$
(7)
$$\lambda _1^{}=\underset{j}{}\frac{1\left|H\right|A_jA_j\left|H\right|1}{\lambda _{A_j}}$$
(8)
$$c_j^{}=\frac{1}{\lambda _j}\underset{k}{}\frac{j\left|H\right|A_kA_k\left|H\right|1}{\lambda _{A_k}}.$$
(9)
The second order contributions are

$$c_{A_j}^{\prime \prime }=\underset{kj}{}\frac{A_j\left|H\right|A_kA_k\left|H\right|1}{\lambda _{A_j}\lambda _{A_k}}\underset{l1}{}\underset{k}{}\frac{A_j\left|H\right|ll\left|H\right|A_kA_k\left|H\right|1}{\lambda _{A_j}\lambda _l\lambda _{A_k}}$$
(10)
$$\lambda _1^{\prime \prime }=\underset{j}{}\underset{kj}{}\frac{1\left|H\right|A_jA_j\left|H\right|A_kA_k\left|H\right|1}{\lambda _{A_j}\lambda _{A_k}}\underset{j}{}\underset{l1}{}\underset{k}{}\frac{1\left|H\right|A_jA_j\left|H\right|ll\left|H\right|A_kA_k\left|H\right|1}{\lambda _{A_j}\lambda _l\lambda _{A_k}}$$
(11)
$`c_m^{\prime \prime }`$ $`={\displaystyle \underset{j}{}}{\displaystyle \underset{kj}{}}\frac{m\left|H\right|A_jA_j\left|H\right|A_kA_k\left|H\right|1}{\lambda _m\lambda _{A_j}\lambda _{A_k}}`$ (12)
$`+{\displaystyle \underset{j}{}}{\displaystyle \underset{l1}{}}{\displaystyle \underset{k}{}}\frac{m\left|H\right|A_jA_j\left|H\right|ll\left|H\right|A_kA_k\left|H\right|1}{\lambda _m\lambda _{A_j}\lambda _l\lambda _{A_k}}`$
$`\frac{1}{\lambda _m}\left[{\displaystyle \underset{j}{}}\frac{1\left|H\right|A_jA_j\left|H\right|1}{\lambda _{A_j}}\right]\frac{1}{\lambda _m}\left[{\displaystyle \underset{k}{}}\frac{m\left|H\right|A_kA_k\left|H\right|1}{\lambda _{A_k}}\right].`$
These contributions are calculated by a Monte Carlo sampling process. Let $`P(A_{\text{trial}(i)})`$ denote the probability of selecting $`A_{\text{trial}(i)}`$ on the $`i^{\text{th}}`$ trial. For example the first order correction to the eigenvalue can be calculated using

$`\lambda _1^{}`$ $`={\displaystyle \underset{j}{}}\frac{1\left|H\right|A_jA_j\left|H\right|1}{\lambda _{A_j}}`$ (13)
$`=\underset{N\mathrm{}}{lim}\frac{1}{N}{\displaystyle \underset{i=1,\mathrm{},N}{}}\frac{1\left|H\right|A_{\text{trial}(i)}A_{\text{trial}(i)}\left|H\right|1}{\lambda _{A_{\text{trial}(i)}}P(A_{\text{trial}(i)})}.`$ (14)
Application of the series method for stochastic error correction to the Hubbard model was discussed by Nathan Salwen at the Heidelberg meeting and is presented in another article in these proceedings.

### Acknowledgments

The author thanks the organizers of the Light-Cone Meeting at Heidelberg and collaborators on the original papers cited here. Financial support provided by the National Science Foundation.


text:
# Secondary CMB Anisotropies from Cosmological Reionization

## 1 Introduction

The reionization of the intergalactic medium (hereafter, IGM) stands out as one of the most important physical processes that have taken place in the early universe. Recent observational progress measuring the angular spectrum of anisotropies in the Cosmic Microwave Background radiation (hereafter CMB) on sub-degree angular scales (de Bernardis et al. 2000; Hanany et al. 2000), and in obtaining the absorption spectra of the most distant quasars (Stern et al. 2000; Zheng et al. 2000; Fan et al. 2000a) allows us to limit the redshift of reionization to somewhere between 6 and 30 (Griffiths, Barbosa, & Liddle 1998; Bond & Jaffe 1999; Tegmark & Zaldarriaga 2000). Combined with recent theoretical and numerical breakthroughs in our understanding how reionization proceeds in the inhomogeneous universe (Giroux & Shapiro 1996; Haiman & Loeb 1997; Gnedin & Ostriker 1997; Madau, Haardt, & Rees 1999; Miralda-Escudeé, Haehnelt, & Rees 1999; Ciardi et al. 1999; Valageas & Silk 1999; Chiu & Ostriker 2000; Gnedin 2000), these data finally place a study of reionization on much more solid footing.

Among the many effects of reionization, the secondary CMB anisotropies have been the subject of active study for a long time, beginning with the pioneering works of Ostriker & Vishniac (1986) and Vishniac (1987) and up to the most recent careful investigations (Hu & White 1996; Aghanim et al. 1996; Knox, Scoccimarro, & Dodelson 1998; Gruzinov & Hu 1998; Jaffe & Kamionkowski 1998; Peebles & Juszkiewicz 1998; Haiman & Knox 1999; Hu 2000; Bruscoli et al. 2000; da Silva et al. 2000; Refregier et al. 2000; Seljak, Burwell, & Pen 2000; Springel, White, & Hernquist 2000). However, most of the recent studies have not benefited from the recent progress on reionization, and were based on over-simplified ad hoc models. In addition, many of these studies only considered the formerly favored “Standard” Cold Dark Matter (CDM) model. However, there appears finally to be convergence toward a low density CDM model with the cosmological constant (the so-called CDM+$`\mathrm{\Lambda }`$ model; c.f. White, Scott, & Pierpaoli 2000; Tegmark & Zaldarriaga 2000; Bridle et al. 2000; Hu et al. 2000; Jaffe et al. 2000), and thus it would make sense to calculate the secondary anisotropies for this currently most favored model.

In this paper we reconsider the effect of reionization on the CMB, using the latest numerical simulations of reionization (Gnedin 2000). By combining our numerical results with analytical calculations and with recent numerical calculation by Springel et al. (2000) on larger angular scales, we are able to extend the range of angular scales over which predictions are reliable to over more than six orders of magnitude. We also take special care to make sure that our results are not contaminated by numerical artifacts, the dominant of which is the periodicity of the simulation boundary conditions.

## 2 Theory

### 2.1 Secondary CMB Anisotropies

Thomson scattering of CMB photons by ionized gas is usually called the Sunyaev-Zel’dovich (SZ) effect, which is often divided into “thermal” and “kinetic” SZ effects. The “thermal” SZ effect is the upscattering of CMB photons after compton scattering by hot gas (most pronounced in clusters of galaxies). The “kinetic” SZ effect is the temperature fluctuation induced by bulk motions of the gas. Both effects play a role in generating secondary CMB anisotropies during and after reionization.

The fractional temperature perturbation in the direction $`\widehat{\theta }`$ induced by bulk motions (“kinetic effect”) is

$$\frac{\mathrm{\Delta }T}{T}(\stackrel{}{\theta })=_0^{\eta _0}n_e\sigma _Te^\tau [\widehat{\theta }\stackrel{}{v}(w\widehat{\theta };w)]a(w)𝑑w,$$
(1)
where $`n_e`$ is the electron density along the line of sight, $`\stackrel{}{v}(\stackrel{}{w};w)`$ is the bulk peculiar velocity at position $`\stackrel{}{w}`$ at a conformal time $`\eta _0w`$, $`\sigma _T`$ is the Thomson cross section, $`\tau `$ is the optical depth from us to $`w`$, and the subscript 0 denotes a quantity at the current moment. The factor of $`a(w)`$ arises because the physical time is $`dt=a(w)dw`$. Note that $`\widehat{\theta }`$ represents a three-dimensional unit vector along the line of sight, whereas $`\stackrel{}{\theta }`$ will refer to a dimensionless two-dimensional vector in the plane perpendicular to the line of sight—i.e., for directions near $`\widehat{z}`$, $`\widehat{\theta }=(\theta _1,\theta _2,\sqrt{1\theta _1^2\theta _2^2})(\theta _1,\theta _2,1)`$, whereas $`\stackrel{}{\theta }=(\theta _1,\theta _2,0)`$.

It is convenient to rewrite equation (1) in the following dimensionless form:

$$\frac{\mathrm{\Delta }T}{T}(\stackrel{}{\theta })=\tau _0_0^{\eta _0}\chi e^\tau (\widehat{\theta }\frac{\stackrel{}{v}}{c})\frac{dw}{a^2\eta _0},$$
(2)
where $`\tau _0c\eta _0\sigma _T(\overline{n}_{\mathrm{H},0}+2\overline{n}_{\mathrm{He},0})`$ and

$$\chi \frac{n_e}{\overline{n}_\mathrm{H}+2\overline{n}_{\mathrm{He}}}=x_e(1+\delta ).$$

Here $`x_e`$ is the ionization fraction, $`\delta `$ is cosmic overdensity, and $`\overline{n}_\mathrm{H}`$ and $`\overline{n}_{\mathrm{He}}`$ are the mean hydrogen and helium number densities.

Analogously, the “thermal” SZ effect is described by the following integral (for observations in the Rayleigh-Jeans regime):

$$\frac{\mathrm{\Delta }T}{T}(\stackrel{}{\theta })=2\tau _0_0^{\eta _0}\chi e^\tau \frac{k_B(T_gT_0)}{m_ec^2}\frac{dw}{a^2\eta _0},$$
(3)
where $`T_g`$ and $`T_0`$ are the gas and the CMB temperatures respectively.

For a particular realization, we can simply calculate integral (1) or (3) along each line of sight and make a map of temperature anisotropy. In this paper we use numerical simulation of cosmological reionization (described in detail in the following section) to directly calculate the anisotropy map on the patch on the sky subtended by the computational box.

Often however it is just the two-point correlation function of the anisotropies, $`C(\theta )`$, or its spherical-harmonic transform, the power spectrum, $`C_{\mathrm{}}`$, which are of interest. In this case we can directly calculate the correlation function from equation (1),

$$C(\theta )\frac{\mathrm{\Delta }T}{T}(\stackrel{}{\theta }_1)\frac{\mathrm{\Delta }T}{T}(\stackrel{}{\theta }_2)_{\mathrm{cos}\theta =\stackrel{}{\theta }_1\stackrel{}{\theta }_2},$$
(4)
which can be reduced to the following form:

$$C(\theta )=\tau _0^2_0^{\eta _0}\frac{dw_1}{a_1^2\eta _0}e^{\tau _1}_0^{\eta _0}\frac{dw_2}{a_2^2\eta _0}e^{\tau _2}$$

$$\times C_{ee}(\eta _1,\eta _2,\theta )C_{vv}(\eta _1,\eta _2,\theta ),$$
(5)
where $`C_{ee}`$ and $`C_{vv}`$ are the electron density and velocity correlation functions,

$`C_{ee}(\eta _1,\eta _2,\theta )`$ $``$ $`\chi (\stackrel{}{r}_1,\eta _1)\chi (\stackrel{}{r}_2,\eta _2),`$
$`C_{vv}(\eta _1,\eta _2,\theta )`$ $``$ $`\widehat{\theta }_1{\displaystyle \frac{\stackrel{}{v}}{c}}(\stackrel{}{r}_1,\eta _1)\widehat{\theta }_2{\displaystyle \frac{\stackrel{}{v}}{c}}(\stackrel{}{r}_2,\eta _2),`$ (6)
with $`\mathrm{cos}\theta =(\stackrel{}{r}_1\stackrel{}{r}_2)/(r_1r_2)`$, provided we make a usual assumption (c.f. Bruscoli et al. 1999) that the electron density and velocity are uncorrelated, i.e.,

$$C_{ev}(\eta _1,\eta _2,\theta )\chi (\stackrel{}{r}_1,\eta _1)\widehat{\theta }_2\frac{\stackrel{}{v}}{c}(\stackrel{}{r}_2,\eta _2)=0.$$
(7)
We will return to the accuracy of this assumption below.

### 2.2 Terminology for Secondary Anisotropies

In reality we have to deal with only one universe, and thus we can only observe the total CMB anisotropies given by equation (1), generically labeled as the kinetic Sunyaev-Zel’dovich effect. Historically, this was first investigated by restricting to a much simpler, unphysical situation: homogeneous reionization and the linear evolution of velocity and density fields. This is the “Ostriker-Vishniac effect” (OV; Ostriker & Vishniac 1986; Jaffe & Kamionkowski 1998), also sometimes called the “Vishniac effect.” In this case, we observe the autocorrelation of the total gas density and velocity fields, related to the power spectrum of density fluctuations by linear theory.

In this work, we will further distinguish between the “Linear Ostriker-Vishniac Effect” (LOV) and the “Nonlinear Ostriker-Vishniac effect” (NLOV). The former is the usual calculation of OV, using the linear relationship between velocity and density. This comes into the needed calculation of the four-point function $`v\delta v\delta `$, in which it is assumed that (1) the linear relationship between $`v`$ and $`\delta `$ holds; and (2) that the four-point function is appropriate for an underlying Gaussian density field, which only holds in the fully linear regime (and only with Gaussian initial conditions). The NLOV assumes a homogeneous ionization fraction, but the full nonlinear gas density and velocity fields.

Recently, the LOV calculations have been complimented by consideration of the distribution of the ionized hydrogen in the universe as it undergoes reionization, so-called “patchy” or “inhomogeneous” reionization. Workers in the field have considered models for the distribution of the ionization fraction in the linear regime (e.g., Knox, Scoccimarro, & Dodelson 1998; Gruzinov & Hu 1998), as well as full simulations or semianalytic models of the ionized gas, as in this work (see also Bruscoli et al. (2000), Benson et al. (2000), and Springel et al (2000)).

This separation is obviously somewhat artificial. Nevertheless, for academic purposes, we can denote those two effects in the following fashion:

$$\frac{\mathrm{\Delta }T}{T}|_{OV}=\tau _0_0^{\eta _0}\overline{x}_e(1+\delta )e^\tau (\widehat{\theta }\frac{\stackrel{}{v}}{c})\frac{dw}{a^2\eta _0}$$
(8)
for the NL Ostriker-Vishniac effect, and

$$\frac{\mathrm{\Delta }T}{T}|_{PR}=\tau _0_0^{\eta _0}(x_e\overline{x}_e)(1+\delta )e^\tau (\widehat{\theta }\frac{\stackrel{}{v}}{c})\frac{dw}{a^2\eta _0}$$
(9)
for patchy reionization.

Thus, the NL Ostriker-Vishniac effect describes anisotropies generated in the universe in which the ionization fraction is homogeneous in space (“homogeneous reionization”), and patchy reionization includes the rest of the anisotropies. We again emphasize that this separation is artificial and unphysical, and we adopt it purely for historical reasons, as the two effects were often considered to be “distinct” effects.

Using the above definition, we can also define the electron density correlation functions due to NL Ostriker-Vishniac and patchy reionization effects,

$`C_{ee}^{OV}(\eta _1,\eta _2,\theta )`$ $``$ $`\overline{x}_e^2\delta (\stackrel{}{r}_1,\eta _1)\delta (\stackrel{}{r}_2,\eta _2),`$
$`C_{ee}^{PR}(\eta _1,\eta _2,\theta )`$ $``$ $`\mathrm{\Delta }\chi (\stackrel{}{r}_1,\eta _1)\mathrm{\Delta }\chi (\stackrel{}{r}_2,\eta _2),`$ (10)
where $`\mathrm{\Delta }\chi (x_e\overline{x}_e)(1+\delta )`$. The total correlation function $`C_{ee}`$ is then a sum of these two functions plus the cross-correlation,

$$C_{ee}=C_{ee}^{OV}+C_{ee}^{PR}+2C_{ee}^{OVPR},$$

where

$$C_{ee}^{OVPR}(\eta _1,\eta _2,\theta )\overline{x}_e\delta (\stackrel{}{r}_1,\eta _1)\mathrm{\Delta }\chi (\stackrel{}{r}_2,\eta _2).$$
(11)
## 3 Method

### 3.1 Simulation

We use a cosmological simulation of reionization reported in Gnedin (2000). The simulation includes 3D radiative transfer (in an approximate implementation) and other physical ingredients required for modeling the process of cosmological reionization.

The simulation of a representative CDM+$`\mathrm{\Lambda }`$ cosmological model111With the following cosmological parameters: $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.7`$, $`\mathrm{\Omega }_b=0.04`$, $`n=1`$, $`\sigma _8=0.91`$, where the amplitude and the slope of the primordial spectrum are fixed by the COBE and cluster-scale normalizations. was performed in a comoving box with the size of $`4h^1\text{ Mpc}`$ with the mass resolution of $`5\times 10^5\text{ M}_{}`$ in baryons and the comoving spatial resolution of $`1h^1\text{ kpc}`$.

The simulation was stopped at $`z=4`$ because at this time the rms density fluctuation in the computational box is about 0.25, and at later times the box ceases to be a representative region of the universe.

As was noted in Gnedin (2000), this simulation is still insufficient to give the fully converged numerically result. As we show below, the computational box size of this simulation is not large enough to allow for accurate computation of the CMB anisotropies, as perturbations at $`z<4`$ contribute about 20% of the signal. We therefore must emphasize here that our calculations are valid only on a semi-qualitative basis, within a factor of 1.5-2, and a still larger simulation is required to accurately model all the relevant scales present in the problem. But since the observations of the secondary anisotropies are some years away, theorists have time to improve upon their models and perform larger, highly accurate simulations.

Our simulation assumes that all cosmological reionization occurs via radiation from stellar sources, with star formation parameterized by the phenomenological Schmidt law as discussed in Gnedin (2000). We ignore alternative and complimentary scenarios in which the bulk of reionization is produced by Active Galactic Nuclei. In light of recent high-redshift quasar counts (e.g., Fan et al. 2000b), a scenario in which the universe is reionized by bright optically selected QSOs seems unlikely in any event. There remains however the possibility that low brightness AGNs contributed significantly (if not dominantly) to the reionization of the universe. In the latter case however they will be clustered on scales similar to the stellar sources, and thus will not result in a qualitatively different reionization scenario, although the epoch of reionization in this case is less well-determined.

### 3.2 Numerical Issues

Given a cosmological simulation, we are now faced with the task of calculating the secondary CMB anisotropies accurately from the simulation. Since perturbations are generated over a considerable redshift range, ideally we would like to have a simulation with the box size that extends from $`z=0`$ until $`z30`$ in one direction. However, with currently feasible simulations this is impossible, and we utilize the cubic box with the $`4h^1\text{ Mpc}`$ side used in Gnedin (2000).

If we take the output of a simulation as is, the artificial periodicity of the universe will amplify the anisotropies by a large factor, simply because in real universe the signal will be averaged over many randomly placed Hii regions, while in a simulated periodic universe a photon, on its way along the line of sight, will encounter the same structure over and over again. This case is depicted in the first two rows of Figure 1. In order to avoid this artificial amplification, we randomly flip and transpose the computational box around any of its six edges, and in addition shift it by a random distance in a random direction. This “randomization” procedure ensures that the next periodic image of the computational box is not correlated with the given image, and thus we lose the correlations over the scales larger than half the size of the computational box along the line-of-sight direction. (In addition to the flipping, transposition and shifting of the simulation box, we have also rotated it by a random amount, but this has negligible effect). This is sketched in the third row of Fig. 1. Thus, our procedure actually underestimates true anisotropies by ignoring large-scale signal. However, we can add the missing large-scale signal using the linear theory calculation, as we discuss below.

In addition, since we cannot unambiguously eliminate periodicity in the plane of the sky, we restrict our image to the image of one periodic simulation box, which corresponds to the angular size of

$$\theta _{\mathrm{max}}=\frac{L}{R(t_{\mathrm{im}})}$$

where $`L=4h^1\text{ Mpc}`$ is the comoving size of our computational box, $`R(t)`$ is the comoving angular-diameter distance from time $`t`$ to the current moment,

$$R=c_t^{t_0}\frac{dt}{a(t)},$$

and we choose to scale all the computational boxes along the line-of-sight to the box taken at $`a_{\mathrm{im}}a(t_{\mathrm{im}})=1/(1+z_{\mathrm{im}})=0.115`$, the epoch just before reionization which contributes most to the CMB anisotropies. For our calculations, $`\theta _{\mathrm{max}}60`$ arcsec corresponding to a spherical harmonic multipole of $`\mathrm{}10^4`$.

Because in the simulation the bulk velocity of the computational box is set to zero, whereas a volume with the size of our box in real universe will have some non-vanishing bulk velocity, we add to each box a random bulk velocity drawn from the Gaussian distribution with the dispersion equal to the linear rms velocity on the scale of the computational box at each redshift. Because this effect depends on the correlation of density with velocity, our results on all angular scales depend on such an appropriate realization of the velocities.

In summary, our procedure to calculate the CMB anisotropies from the cosmological simulation with a periodic computational box can be described as follows:

1. We select the outputs from the simulation spaced with $`\mathrm{\Delta }a`$ in scale factor.
2. For each pair of outputs at $`a_i`$ and $`a_{i+1}`$, where index $`i`$ numbers the outputs, we calculate the number of computational boxes along the line of sight between $`a_i`$ and $`a_{i+1}`$.
3. Then for all boxes between $`a_i`$ and $`(a_i+a_{i+1})/2`$ we fill them in with the simulation output data at $`a_i`$, and the boxes between $`(a_i+a_{i+1})/2`$ and $`a_{i+1}`$ we fill in with the data at $`a_{i+1}`$.
4. Then each of thus-filled-in boxes (including the original simulation boxes at $`a_i`$ and $`a_{i+1}`$) is randomized by a random shift in the plane of the sky and by a random flip and transposition, removing correlations on scales larger than half the box size. In addition, each box is assigned a bulk velocity along the line-of-sight according to the linear theory.
5. We repeat steps 2-4 for all pairs of outputs in the descending order from the last output at $`z=4`$ ($`a=0.2`$) until the start of the simulation at $`100`$.
6. After arranging all the boxes, we calculate the CMB map on a square patch on the sky (corresponding to the image of our cubic computational box at $`z=7.7`$, the epoch when about half the box is ionized) with a given resolution (we achieve the highest resolution of $`512^2`$ pixels).

This procedure allows us to construct an image on the sky. However, it is customary to present the CMB anisotropies on the sky in terms of the power spectrum $`C_{\mathrm{}}`$, which is the spherical-harmonic transform of the correlation function $`C(\theta )`$. In the small-angle limit, the spherical-harmonic transform reduces to the Hankel transform of zeroth order,

$$C_{\mathrm{}}=2\pi _0^{\mathrm{}}𝑑\theta \theta J_0(\theta \mathrm{})C(\theta ),$$
(12)
where $`J_0`$ is the zero order Bessel function. However, because of the numerical errors and artificial ringing due to a finite box size, the power spectrum calculated from equation (12) is noisy. An alternative method of calculating the power spectrum is a direct Fourier transform of the CMB anisotropies on the sky,

$$C_{\mathrm{}}=\theta _{\mathrm{max}}^2|\mathrm{\Delta }_{\stackrel{}{\mathrm{}}}|^2,$$
(13)
where $`\theta _{\mathrm{max}}`$ is the angular size of a square patch on the sky (obviously, $`\theta _{\mathrm{max}}1`$), and

$$\mathrm{\Delta }_{\stackrel{}{\mathrm{}}}\frac{1}{\theta _{\mathrm{max}}^2}\frac{\mathrm{\Delta }T}{T}(\stackrel{}{\theta })e^{i\stackrel{}{\theta }\stackrel{}{\mathrm{}}}d^2\theta $$

is the Fourier transform of the CMB anisotropies.

Figure 2 shows a comparison between the two methods for computing the power spectrum $`C_{\mathrm{}}`$. One can see that the correlation function method produces an extremely noisy estimate of the power spectrum, whereas the one computed directly is well defined up to the maximal possible multipole

$$\mathrm{}_{\mathrm{max}}=\sqrt{2}\mathrm{}_{\mathrm{Ny}}=\sqrt{2}N\frac{\pi }{\theta _{\mathrm{max}}},$$

where $`\mathrm{}_{\mathrm{Ny}}`$ is the Nyquist frequency of a square image on the sky and $`N`$ is the number of pixels along one dimension in the image.

The process of pixelization - binning of the computed image on the square patch on the sky of a given size - produces small-scale smearing similar to the one produced by assigning density on a regular mesh from a particle distribution in cosmological $`N`$-body simulations. However, in our case the precise form of this smoothing cannot be computed from the first principles, because the computational box subtends somewhat different angular size at different redshifts. Instead, we correct for the small-scale smearing empirically. Let $`C_{\mathrm{}}`$ be the true power spectrum. Since the smearing is caused by pixelization, the computed power spectrum $`\stackrel{~}{C}_{\mathrm{}}`$ will be a product of the true power spectrum and the window function $`W`$ which depends only on the ratio of the multipole $`\mathrm{}`$ and the maximal possible multipole $`\mathrm{}_{\mathrm{max}}`$ for a given pixel size, or, in other words, on the product of $`\mathrm{}`$ and the number of pixels along one dimension $`N`$. Thus,

$$\stackrel{~}{C}_{\mathrm{}}=W(\mathrm{}/\mathrm{}_{\mathrm{max}})C_{\mathrm{}},$$

and if we take the ratio of the power spectra from two different pixelizations, this ratio will depend only on the window function,

$$\frac{\stackrel{~}{C}_{1,\mathrm{}}}{\stackrel{~}{C}_{2,\mathrm{}}}=\frac{W(\mathrm{}/\mathrm{}_{1,\mathrm{max}})}{W(\mathrm{}/\mathrm{}_{2,\mathrm{max}})}.$$
(14)
Using equation (14), we find the following approximate form for the window function

$$W(\mathrm{}/\mathrm{}_{\mathrm{max}})=\frac{\left[1+10/3(\mathrm{}/\mathrm{}_{\mathrm{max}})^2\right]^2}{1+(\mathrm{}/\mathrm{}_{\mathrm{max}})^4}.$$
(15)
Figure 3 shows the uncorrected and corrected power spectra for four resolutions: $`N=64`$, $`128`$, $`256`$, and $`512`$. One can see that the corrected power spectra agree with each other up to their respective $`\mathrm{}_{\mathrm{max}}`$. We also notice that the power spectrum seems to continue as a power law with index -3 to the smallest scales. Thus, we are unable to achieve convergence on small angular scales, which implies that even with resolution of $`512^2`$ pixels per 2.2 arcmin patch on the sky, there exists structure on the unresolved scales. As we mentioned above, this structure is due to small-scale density inhomogeneity in the ionized high density regions. Our simulation has a dynamical range of 4000, and thus in principle we can measure the power spectrum of the secondary anisotropies from our simulation up to $`l3\times 10^7`$ (compared to $`l=3.5\times 10^6`$ for our highest $`512^2`$ resolution). However, making CMB maps with resolution higher than $`512^2`$ is beyond the limit of computational resources available to us (doing a $`4000^2`$ image would take more computer time than the simulation itself consumed). We thus conclude that the cut-off in the power spectrum at $`l3\times 10^5`$ reported in Bruscoli et al. (2000) is not real and due to the finite angular resolution ($`256^2`$) of their final CMB map.

The outputs from our simulation are only stored at discreet moments in time (which we choose to parametrize with the cosmological scale factor $`a1/[1+z]`$). Integration over the path of a photon is done numerically using the outputs from the simulation, and thus the numerical value of the integral depends on the interval $`\mathrm{\Delta }a`$ between the consequent outputs, converging to the exact value in the limit $`\mathrm{\Delta }a0`$. Of course, we cannot reach this limit as it would require an infinite number of outputs, and thus we must ensure that our choice of $`\mathrm{\Delta }a`$ is sufficiently small to give an accurate value for the integral over the photon path. Figure 4 shows the power spectra of the CMB anisotropies for two choices of $`\mathrm{\Delta }a`$. Both values give convergent results for the CMB anisotropies, and so we adopt $`\mathrm{\Delta }a=0.005`$ as our choice for this parameter (which is equivalent to having about 40 outputs until $`z=4`$). It also becomes prohibitively expensive to calculate the integral with a smaller $`\mathrm{\Delta }a`$. We also show in Fig. 4 with the long-dashed line the CMB anisotropies calculated without randomization (i.e., in an unrealistic periodic universe), which are off by up to a factor of 10.

We would like to note here that a recent similar analysis by Bruscoli et al. (2000) did not use randomization in calculating the CMB anisotropies, and by adopting a periodic universe, overestimated the correct result by a factor of 3-10.

Finally, since our simulation is stopped at $`z=4`$, we need to estimate the contribution from the lower redshifts. Since the universe is fully ionized at $`z=4`$, this contribution is entirely due to a homogeneous ionization fraction, $`x_e=1`$ (the NL Ostriker-Vishniac effect). We can calculate it by extrapolating from higher redshifts. Figure 5a shows the rms CMB temperature anisotropy as a function of the final redshift of calculation $`z`$ from our simulation (up to $`z=4`$) with the solid line. One can see that even at $`z=4`$ the rms temperature anisotropy continues to rise (despite the fact that the universe is fully ionized by then), and thus we miss a portion of the total signal. In order to estimate the missing signal, we extrapolate the signal to $`z<4`$ assuming that the density and velocity structure remains fixed in the comoving coordinates (in other words, taking the output of the simulation at $`z=4`$ and adopting it for all lower redshifts). Then we find that

$$\frac{\mathrm{\Delta }T}{T}|_{\mathrm{RMS}}(z=0)=1.25\frac{\mathrm{\Delta }T}{T}|_{\mathrm{RMS}}(z=4),$$

i.e., we miss about 20% of the signal. This number is likely to be an overestimate, as the gas pressure effects will erase the structure on spatial scales up to $`500\text{ kpc}`$ by $`z=0`$, which corresponds to angular scales which are much larger than the ones dominating the signal. Nevertheless, since we cannot continue the simulation beyond $`z=4`$, we use the extrapolated power spectrum as our final result. The two spectra, the extrapolated one and the one directly computed from the simulation, are shown in Fig. 5b. The difference between the two is larger than, albeit comparable to, the statistical error due to the finite box size of the simulation, and is the largest uncertainty of our calculation on small angular scales (must less than the size of the box). On the angular scales comparable to the box size, an additional uncertainty is introduced by the missing large-scale power.

### 3.3 Perturbations Outside of the Simulation Box

Since the size of the simulation box is fixed to $`4h^1\text{ Mpc}`$ in comoving units, we are unable to calculate directly anisotropies due to perturbations on larger scales. However, we can include them analytically using the linear theory. Since we are not able to follow the ionization state of the gas on these scales, we assume that in the linear regime the ionization fraction is spatially homogeneous. This is not a bad approximation after all, since the size of a typical Hii region is smaller than the simulation box size, and thus on larger scales we average over a sufficiently large volume of the universe. The latter implies that in the linear regime we can consider only the Ostriker-Vishniac effect – the NLOV effect is the same as the LOV effect. As we argue below, it is the dominant contribution to the total anisotropies on all angular scales.

To calculate the pure LOV effect (linear fluctuations and homogeneous reionization), we use the approach described in Jaffe & Kamionkowski (1998). We use a small-angle Fourier-space version of Limber’s equation which decouples the line-of-sight projection from the three-dimensional power spectra. As discussed above, we also assume the linear relationship between density and velocity, uniform reionization, and Gaussian statistics. This reduces the calculation to a multiple integral over a product of the linear power spectrum of density perturbations with itself.

To check our methods against this analytic calculation, we perform an additional “simulation” of the LOV effect in which we assume that the density and velocity fields evolve only according to linear theory, and we assign the ionization fraction at a given redshift uniformly over the computational box with the value given by the volume average ionization fraction in the full numerical simulation. We then make an image and calculate the CMB power spectrum in precisely the same fashion as in the full simulation. Thus, the spectrum of the anisotropies from a such “linear” simulation includes all the finite-box effects present in the full numerical simulation.

The full nonlinear power spectrum is shown with the bold solid line in Fig. 6. In addition, we show the output from our linear “simulation” along with the analytic OV calculation, which match over all scales computed. We emphasize that this agreement requires the imposition of appropriate velocities on scales larger than the simulation box, as discussed above. We also note that our full nonlinear prediction (bold solid line) matches the linear theory on the angular scale corresponding to the box size. This is of course is an artifact of our procedure (and the limited size of our simulation box), which does not include nonlinear perturbations outside the simulation box. In reality, there exists additional power on scales which we can not resolve, simply because the nonlinear scale at $`z=0`$ is about twice larger than our whole box. To account for this power, we show the nonlinear power spectrum computed by Springel et al. (2000). They assumed precisely the same cosmological model and the redshift of reionization. In addition, their mass resolution almost exactly corresponds to 1/8 of the total mass of our simulation box (corresponding to 1/2 the box size - the largest scale which can be resolved in our simulation). Thus, Springel et al. (2000) begin precisely at the scale where our calculation ends, and we will use their results to compliment ours below. (Conversely, this means that we unfortunately cannot actually check the calculations against one another.)

## 4 Results

### 4.1 Maps and non-Gaussianity

It is customary to characterize the CMB temperature anisotropies on the sky by the power spectrum $`C_{\mathrm{}}`$. The power spectrum is sufficient to fully describe properties of a physical quantity only if this quantity is Gaussian distributed. For the primary temperature anisotropies this is indeed the case, since they are generated in the linear regime. Reionization however is a non-linear process, and secondary anisotropies we consider in this paper are not necessarily Gaussian.

This is illustrated in Figure 7, which shows the actual maps of the CMB anisotropies on the square patch on the sky in the simulation. The left panel shows the actual anisotropies on 2.2 arcmin patch on the sky (the angular size of our computational box at $`z=8`$, just before the full overlap of Hii regions). The right panel, however, shows anisotropies with the same power spectrum but distributed as a Gaussian. We notice that the actual map is highly non-Gaussian: it is much more structured than the Gaussian version and contains higher fluctuations than would be expected in a Gaussian field.

To illustrate this further, we show in Figure 8 the one-point distribution of the actual temperature fluctuations, and a Gaussian distribution for comparison. We notice that the temperature anisotropies are non-Gaussian, and thus more information than simply $`C_{\mathrm{}}`$ can be extracted from the actual distribution on the sky. In addition to temperature anisotropies on the sky, we also show in Fig. 8 the one-point distribution function of spherical multipole amplitudes $`a_\mathrm{}m`$, which are defined in the usual fashion by expansion over spherical harmonics,

$$\frac{\mathrm{\Delta }T}{T}(\theta ,\varphi )=\underset{\mathrm{},m}{}a_\mathrm{}mY_{\mathrm{}}^m(\theta ,\varphi ).$$

The one-point distribution of the amplitudes is actually quite similar to a Gaussian, which indicates that it is indeed the the phase correlations between different multipoles what drive the non-Gaussian patterns. In order to characterize this non-Gaussianity, one would have to measure the higher-order moments (bispectrum, trispectrum, etc.) at non-zero lags, or other quantities such as the Minkowski functionals, but that is outside the scope of this paper.

### 4.2 Correlation Functions

We now investigate the behavior of the correlation functions. In particular, the electron density correlation function $`C_{ee}(\eta _1,\eta _2,\theta )`$ is of highest interest. We show in Figure 9 the electron density correlation function at three redshifts: before the overlap of Hii regions at $`z=10`$, during the overlap at $`z=7`$, and after the overlap at $`z=5`$. One can notice that there is no drastic change in the shape of the correlation function despite the qualitative change in the structure of the ionized regions; only the correlation length increases with time. This indicates that the signal is dominated by the ionized high density regions at all times, and is not sensitive to the ionization state of the low density intergalactic medium. We will elaborate on this below.

Figure 10 shows both the electron density correlation function $`C_{ee}`$ and the velocity correlation function $`C_{vv}`$ as a function of angle (measured as comoving distance at this redshift) at $`z=7`$. We also show with the dotted line the quantity $`C_{ev}^2/C_{vv}`$, which is much smaller than $`C_{ee}`$ if the cross-correlation between the electron density and velocity is not significant, and is equal to $`C_{ee}`$ is they are completely correlated. One can see that the cross-correlation only becomes significant on very small scales, about 0.1 arcsec ($`\mathrm{}6\times 10^6`$), and thus the usual assumption that velocity and electron density are uncorrelated is highly accurate on all scales of interest.

### 4.3 Kinetic vs Thermal SZ Effect

It is interesting to compare the kinetic SZ effect given by equation (1) with the thermal SZ effect, given by equation (3). One can expect a priori that the kinetic SZ effect, generated in the gas moving with some $`300\text{ km}/\text{ s}`$ is about 1000 times more important (in terms of the power spectrum) than the thermal SZ effect, generated due to thermal motions in the $`10^4\text{ K}`$ hot gas. Indeed, this is demonstrated in Figure 11. Thus, for all practical purposes, the kinetic effect is the only one that needs to be taken into account on these scales. Of course, for clusters of galaxies, the situation is reversed and the temperature fluctuation induced by even the several-hundred km/s bulk motion of the cluster is dwarfed by the thermal effect, simply because the cluster gas temperature is some four orders of magnitude higher than the IGM temperature, whereas the bulk gas velocity in a cluster is only about a few times that of the IGM. We return to the comparison between the two effects below.

### 4.4 The Ostriker-Vishniac Effect versus Patchy Reionization

We now focus on a comparison between the NL Ostriker-Vishniac effect and patchy reionization. Figure 12 shows the comparison between the two correlation functions $`C_{ee}^{PR}`$ and $`C_{ee}^{OV}`$, which are defined by equations (2.2). The dotted line, marking the cross-correlation term in comparison with $`C_{ee}^{PR}`$, only falls significantly below $`C_{ee}^{PR}`$ at large radii ($`R>0.1h^1\text{ Mpc}`$, where the inhomogeneous distribution of the ionized fraction due to expanding Hii regions make the NL Ostriker-Vishniac effect and patchy reionization effect uncorrelated. However, as the thin dashed line shows, in the high density regions (i.e., on small scales) the effect of patchy reionization is almost equivalent to having $`x_e=1`$, i.e., to having almost all high density regions ionized. This implies that the bulk of the signal comes from high density regions and correlations between them.

In order to investigate the relationship between the two effects further, we have constructed a simulation which contains only the “Ostriker-Vishniac” effect by assigning a uniform ionization fraction (equal to the volume average at a given redshift) to the outputs of our full numerical simulation, and performing line-of-sight integration as described above. Thus, this “NLOV simulation” has the density and velocity structure of our full numerical simulation, and has the same evolution of the volume averaged ionization fraction, but has uniform ionization at all densities and spatial locations. In particular, the Patchy Reionization correlation function $`C_{ee}^{PR}`$ is equal to zero in that simulation, and its NL Ostriker-Vishniac correlation function ($`C_{ee}^{OV}`$) is the same as in the full simulation. Thus, the NLOV simulation allows us to construct an image of the CMB anisotropies which correspond to just the $`C_{ee}^{OV}`$ correlation function, whereas the full simulations gives us anisotropies corresponding to the sum of the NL Ostriker-Vishniac and Patchy Reionization effects.

Figure 13 now shows the rms temperature anisotropy on the sky (for the $`256^2`$ pixelization) for the total effect, and separately for the NL Ostriker-Vishniac effect (as computed from the “NLOV simulation”) and for the patchy reionization effect (computed as the difference between the two). Note, that by $`z=4`$ almost all the signal comes from the NL Ostriker-Vishniac effect, and the contribution of the patchy reionization actually slightly declines at $`z<7`$, after the overlap of Hii regions. This is due to the fact that after the overlap the topology of the ionized phase flips over: reionization starts with isolated ionized regions expanding into the neutral medium, and ends with the isolated neutral regions being ionized from outside inward, which produces an anti-correlation between the ionized and neutral phases and leads to suppression of anisotropies due to patchy reionization.

Figure 14 presents now the main result of this paper: the power spectrum of the secondary CMB anisotropies from the kinetic SZ effect, from the thermal SZ effect, and from the patchy reionization effect. For the former two, we combined our results with the calculations by Springel et al. (2000). In the range $`3\times 10^4\mathrm{}2\times 10^5`$ our calculation for the kinetic SZ effect disagrees with the extrapolation from the Springel at al. (2000), and for the thermal SZ effect the range of disagreement is even larger. Since our results can suffer from the missing large-scale power at low redshift, and small scale extrapolation of Springel et al. (2000) is likely to be severely affected by limited numerical resolution, it is impossible to give a reliable prediction for the CMB power spectra in this range of angular scales, and we emphasize this fact by using the bold lines to show the calculations where they are reliable, and by using the thin lines in the range of scales where consider calculations to be unreliable.

We also show the prediction from an analytical calculation of the patchy reionization effect by Gruzinov & Hu (1998) for our cosmological model. Note that this prediction is a factor of $`10^2`$ lower than the one given in Gruzinov & Hu (1998) for the standard CDM model, because they adopted the value for the rms velocity at $`z=0`$ of $`1200\text{ km}/\text{ s}`$, whereas in our model this number is about 3 times smaller, in agreement with observations (Baker, Davis, & Lin 2000). There is an additional reduction in the amplitude of the effect by a factor of 10 due to the redshift of reionization being 7 in our simulation compared to an adopted value of 30 in Gruzinov & Hu (1998). Including all the relevant factors, we notice that Gruzinov & Hu (1998) estimate gives roughly the right prediction for the patchy reionization effect on large scales. Of course, since the analytical estimate of Gruzinov & Hu (1998) does not include the density structure on small scales, they are unable to reproduce the anisotropies on small angular scales.

### 4.5 Dependence on the Redshift of Reionization

Up to know we assumed a fixed value for the redshift of reionization, since we only used one simulation. In order to investigate the dependence of the power spectrum on the redshift of reionization properly, we would need to run several simulation with different reionization histories, which is not realistic at the current moment due to required computational resources. However, in order to illustrate the dependence of anisotropies on the redshift of reionization, and using the fact that at $`z_{\mathrm{REI}}=7`$ the NL Ostriker-Vishniac dominates the contribution to the anisotropies due to Patchy Reionization, we have computed three additional “NLOV” simulations where we changed the redshift of reionization by simply rescaling the volume averaged ionization fraction as a function of redshift. In other words, if $`n_{e,0}(a)`$ is the electron density and $`\overline{x}_e(a)`$ is the volume averaged ionization fraction as a function of the scale factor for the original “NLOV” simulation, for rescaled “NLOV” simulations we assumed

$$n_e(a)=\frac{\overline{x}_e(8a_{\mathrm{REI}}a)}{\overline{x}_e(a)}n_{e,0}(a),$$

where the factor 8 is simply $`1/0.125=1/a_{\mathrm{REI}}=1+z_{\mathrm{REI}}`$ for our simulation.

Figure 15 shows the CMB power spectra for four different values for the redshift of reionization in our rescaled “NLOV” simulations. We note that in linear theory, $`\delta a`$ and $`va^{1/2}`$ at high redshift, and therefore the total anisotropy

$$\frac{\mathrm{\Delta }T}{T}^ta^{3/2}\overline{n}_b𝑑t\mathrm{log}(t)$$

goes only logarithmically as a function of the moment of reionization (which is the upper limit in the integral  for a model with sudden reionization). The dependence that we find in the nonlinear calculation is somewhat steeper than that from the linear theory, but not by much.

We also note here that the analytical model of patchy reionization by Gruzinov & Hu (1998) predicts that the spectrum of the CMB anisotropies increases as $`(1+z_{\mathrm{REI}})^{3/2}R_{\text{H}\text{ii}}`$, where $`R_{\text{H}\text{ii}}`$ is the characteristic comoving size of Hii regions, which is likely to go down as the redshift of reionization increases. Thus, the total dependence on the redshift of reionization of the patchy reionization effect is somewhat slower than $`(1+z_{\mathrm{REI}})^{3/2}`$ but perhaps not as slow as the nearly logarithmic dependence of the NL Ostriker-Vishniac effect. Thus, if reionization occurred early, the role of patchy reionization (as compared to the NL Ostriker-Vishniac effect) would be greater, albeit the current observational limits of $`z_{\mathrm{REI}}30`$ will not allow a regime when patchy reionization becomes dominant.

## 5 Conclusions

We use a cosmological numerical simulation that models time-dependent and spatially-inhomogeneous cosmological reionization to compute secondary CMB temperature anisotropies in a representative CDM+$`\mathrm{\Lambda }`$ cosmology. We compare our numerical results with analytical calculations and compliment them with the numerical simulations of Springel et al. (2000) (whose mass resolution almost exactly compliments ours, and accidentally both simulations have the same values of cosmological parameters). We thus are able to compute the power spectrum of secondary anisotropies for the kinetic SZ effect and for the patchy reionization effect over a wide range of angular scales. The thermal SZ effect, dominated at medium $`\mathrm{}`$ by massive clusters of galaxies, is more difficult to compute on all scales due to the extreme non-linearity of the problem.

We find that the role of patchy reionization per se (i.e., evolution of Hii regions around the sources of ionization) in generating the secondary anisotropies is subdominant. This conclusion is however only valid when the sources of reionization are clustered on scales comparable to the matter correlation length, in which case one can expect to find a source (or multiple sources) in almost every dark matter halo. If, in the opposite case, the sources of reionization are rare (like bright quasars), the relative role of patchy reionization will be greater, as many high density regions will remain neutral for a long time, and may not be negligible. The latter scenario is, however, less likely given the current observational limits on the abundances of bright quasars (Fan et al. 2000b), and thus CMB observations on small angular scales may not be as useful as hoped to study the details of reionization.

A real window onto reionization is provided by the detailed morphology of the fluctuations, that is, the non-Gaussianity discussed above. The angular scale of the fluctuations traces the three-dimensional separation scale of the ionized regions. Of course, a full understanding of this will require not just a detection of broad-band fluctuation power at very small angular scales, but a detailed, high signal-to-noise map of these very small-scale fluctuations with an angular resolution of arc seconds, which is definitely some years away.

While the details of reionization may not be easily discernible from the broad-band power spectrum, the amplitude of the fluctuations, which is directly related to the redshift of reionization, surely is. However, as our comparison with Springel et al. (2000) shows, on medium angular scales ($`l10^5`$) the signal is dominated by lower redshifts ($`z24`$), and only on small angular scales (several arc seconds) the amplitude of the anisotropies is related to the redshift of reionization. In addition, this dependence is quite weak, unless reionization occurred close to the currently allowed observational limit of $`z_{\mathrm{REI}}30`$. Nevertheless, observations on those angular scales can put significant constraints on the redshift of reionization. Such arcsecond-scale CMB observations are best accomplished with interferometers, and indeed the recent development of compact submillimeter-wave arrays has paved the way for such observations. For example, the ATCA telescope finds a limit of $`\mathrm{}(\mathrm{}+1)C_{\mathrm{}}/(2\pi )<2\times 10^{10}`$ at $`\mathrm{}4000`$ (Subrahmanyan et al. 2000), well above the values presented here. In addition, these observations are inevitably plagued by point source confusion, especially dangerous when the signal also has a compact galaxy-like component.

We would like to acknowledge helpful conversations with Marc Kamionkowski, Joe Silk, Andrew Liddle, Naoshi Sugiyama and Ue-Li Pen. AHJ acknowledges support from NSF KDI grant 9872979 and NASA LTSA grant NAG5-6552.


text:
# Error and attack tolerance of complex networks

Many complex systems display a surprising degree of tolerance against errors. For example, relatively simple organisms grow, persist and reproduce despite drastic pharmaceutical or environmental interventions, an error tolerance attributed to the robustness of the underlying metabolic network. Complex communication networks display a surprising degree of robustness: while key components regularly malfunction, local failures rarely lead to the loss of the global information-carrying ability of the network. The stability of these and other complex systems is often attributed to the redundant wiring of the functional web defined by the systems’ components. In this paper we demonstrate that error tolerance is not shared by all redundant systems, but it is displayed only by a class of inhomogeneously wired networks, called scale-free networks. We find that scale-free networks, describing a number of systems, such as the World Wide Web (www), Internet, social networks or a cell, display an unexpected degree of robustness, the ability of their nodes to communicate being unaffected by even unrealistically high failure rates. However, error tolerance comes at a high price: these networks are extremely vulnerable to attacks, i.e. to the selection and removal of a few nodes that play the most important role in assuring the network’s connectivity. Such error tolerance and attack vulnerability are generic properties of communication networks, such as the Internet or the www, with complex implications on assuring information readiness.

The increasing availability of topological data on large networks, aided by the computerization of data acquisition, has lead to major advances in our understanding of the generic aspects of network structure and development. The existing empirical and theoretical results indicate that complex networks can be divided into two major classes based on their connectivity distribution $`P(k)`$, giving the probability that a node in the network is connected to $`k`$ other nodes. The first class of networks is characterized by a $`P(k)`$ that is peaked at an average $`k`$ and decays exponentially for large $`k`$. The most investigated examples of such exponential networks are the random graph model of Erdős and Rényi and the small-world model of Watts and Strogatz, both leading to a fairly homogeneous network, in which each node has approximately the same number of links, $`kk`$. In contrast, results on the world-wide web (www), Internet and other large networks indicate that many systems belong to a class of inhomogeneous networks, referred to as scale-free networks, for which $`P(k)`$ decays as a power-law, i.e. $`P(k)k^\gamma `$, free of a characteristic scale. While the probability that a node has a very large number of connections ($`k>>k`$) is practically prohibited in exponential networks, highly connected nodes are statistically significant in scale-free networks (see Fig. 1).

We start by investigating the robustness of the two basic network models, the Erdős-Rényi (ER) model that produces a network with an exponential tail, and the scale-free model with a power-law tail. In the ER model we first define the $`N`$ nodes, and then connect each pair of nodes with probability $`p`$. This algorithm generates a homogeneous network (Fig. 1), whose connectivity follows a Poisson distribution peaked at $`k`$ and decaying exponentially for $`k>>k`$.

The inhomogeneous connectivity distribution of many real networks is reproduced by the scale-free model that incorporates two ingredients common to real networks: growth and preferential attachment. The model starts with $`m_0`$ nodes. At every timestep $`t`$ a new node is introduced, which is connected to $`m`$ of the the already existing nodes. The probability $`\mathrm{\Pi }_i`$ that the new node is connected to node $`i`$ depends on the connectivity $`k_i`$ of that node, such that $`\mathrm{\Pi }_i=k_i/_jk_j`$. For large $`t`$ the connectivity distribution is a power-law following $`P(k)=2m^2/k^3`$.

The interconnectedness of a network is described by its diameter $`d`$, defined as the average length of the shortest paths between any two nodes in the network. The diameter characterizes the ability of two nodes to communicate with each other: the smaller $`d`$ is, the shorter is the expected path between them. Networks with a very large number of nodes can have a rather small diameter; for example the diameter of the www, with over $`800`$ million nodes, is around $`19`$, while social networks with over six billion individuals are believed to have a diameter of around six. To properly compare the two network models we generated networks that have the same number of nodes and links such that $`P(k)`$ follows a Poisson distribution for the exponential, and a power-law for the scale-free network.

Error tolerance— To address the networks’ error tolerance, we study the changes in the diameter when a small fraction $`f`$ of the nodes is removed. The malfunctioning (absence) of a node in general increases the distance between the remaining nodes, since it can eliminate some paths that contribute to the system’s interconnectedness. Indeed, for the exponential network the diameter increases monotonically with $`f`$ (Fig. 2a), thus, despite its redundant wiring (Fig. 1), it is increasingly difficult for the remaining nodes to communicate with each other. This behavior is rooted in the homogeneity of the network: since all nodes have approximately the same number of links, they all contribute equally to the network’s diameter, thus the removal of each node causes the same amount of damage. In contrast, we observe a drastically different and surprising behavior for the scale-free network (Fig. 2a): the diameter remains unchanged under an increasing level of errors. Thus even when as high as $`5\%`$ of the nodes fail, the communication between the remaining nodes in the network is unaffected. This robustness of scale-free networks is rooted in their extremely inhomogeneous connectivity distribution: since the power-law distribution implies that the majority of nodes have only a few links, nodes with small connectivity will be selected with much higher probability, and the removal of these ”small” nodes does not alter the path structure of the remaining nodes, thus has no impact on the overall network topology.

Attack survivability— An informed agent that attempts to deliberately damage a network, such as designing a drug to kill a bacterium, will not eliminate the nodes randomly, but will rather target the most connected nodes. To simulate an attack we first remove the most connected node, and continue selecting and removing nodes in the decreasing order of their connectivity $`k`$. Measuring the diameter of an exponential network under attack, we find that, due to the homogeneity of the network, there is no substantial difference whether the nodes are selected randomly or in decreasing order of connectivity (Fig. 2a). On the other hand, a drastically different behavior is observed for scale-free networks: when the most connected nodes are eliminated, the diameter of the scale-free network increases rapidly, doubling its original value if $`5\%`$ of the nodes are removed. This vulnerability to attacks is rooted in the inhomogeneity of the connectivity distribution: the connectivity is ensured by a few highly connected nodes (Fig. 1b), whose removal drastically alters the network’s topology, and decreases the ability of the remaining nodes to communicate with each other.

Network fragmentation— When nodes are removed from a network, clusters of nodes, whose links to the system disappear, can get cut off from the main cluster. To better understand the impact of failures and attacks on the network structure, we next investigate this fragmentation process. We measure the size of the largest cluster, $`S`$, shown as a fraction of the total system size, when a fraction $`f`$ of the nodes are removed either randomly or in an attack mode. We find that for the exponential network, as we increase $`f`$, $`S`$ displays a threshold-like behavior such that for $`f>f_c0.28`$ we have $`S0`$. A similar behavior is observed when we monitor the average size $`s`$ of the isolated clusters (i.e. all the clusters except the largest one), finding that $`s`$ increases rapidly until $`s2`$ at $`f_c`$, after which it decreases to $`s=1`$. These results indicate the following breakdown scenario (Fig. 4): For small $`f`$, only single nodes break apart, $`s1`$, but as $`f`$ increases, the size of the fragments that fall off the main cluster increases, displaying a singular behavior at $`f_c`$. At $`f_c`$ the system practically falls apart, the main cluster breaking into small pieces, leading to $`S0`$, and the size of the fragments, $`s`$, peaks. As we continue to remove nodes ($`f>f_c`$), we fragment these isolated clusters, leading to a decreasing $`s`$. Since the ER model is equivalent to the infinite dimensional percolation, the observed threshold behavior is qualitatively similar to the percolation critical point.

However, the response of a scale-free network to attacks and failures is rather different (Fig. 3b). For random failures no threshold for fragmentation is observed, rather the size of the largest cluster slowly decreases. The fact that $`s1`$ for most $`f`$ indicates that the network is deflated by nodes breaking off one by one, the increasing error level leading to the isolation of single nodes only, not clusters of nodes. Thus, in contrast with the catastrophic fragmentation of the exponential network at $`f_c`$, the scale-free network stays together as a large cluster for very high values of $`f`$, providing additional evidence of the topological stability of these networks under random failures. This behavior is consistent with the existence of an extremely delayed critical point (Fig. 3), the network falling apart only after the main cluster has been completely deflated. On the other hand, the response to attack of the scale-free network is similar (but swifter) to the response to attack and failure of the exponential network (Fig. 3b): at a critical threshold $`f_c^{sf}0.18`$, smaller than the value $`f_c^e0.28`$ observed for the exponential network, the system breaks apart, forming many isolated clusters (Fig. 4).

While great efforts are being made to design error tolerant and low yield components for communication systems, little is known about the effect of the errors and attacks on the large-scale connectivity of the network. To demonstrate the impact of our model based studies to these systems, next we investigate the error and attack tolerance of two networks of increasing economic and strategic importance: the Internet and the www.

Recently Faloutsos et al. investigated the topological properties of the Internet at the router and inter-domain level, finding that the connectivity distribution follows a power-law, $`P(k)k^{2.48}`$. Consequently, we expect that it should display the error tolerance and attack vulnerability predicted by our study. To test this, we used the latest survey of the Internet topology, giving the network at the inter-domain (autonomous system) level. Indeed, we find that the diameter of the Internet is unaffected by the random removal of as high as $`2.5\%`$ of the nodes (an order of magnitude larger than the failure rate ($`0.33\%`$) of the Internet routers), while if the same percentage of the most connected nodes are eliminated (attack), $`d`$ more than triples (Fig. 2b). Similarly, the large connected cluster persists for high rates of random node removal, but if nodes are removed in the attack mode, the size of the fragments that break off increases rapidly, the critical point appearing at $`f_c^I0.03`$ (Fig. 3b).

The www forms a huge directed graph whose nodes are documents and edges are the URL hyperlinks that point from one document to another, its topology determining the search engines’ ability to locate information on it. The www is also a scale-free network: the probabilities $`P_{out}(k)`$ and $`P_{in}(k)`$ that a document has $`k`$ outgoing and incoming links follow a power-law over several orders of magnitude, i.e. $`P(k)k^\gamma `$, with $`\gamma _{in}=2.1`$ and $`\gamma _{out}=2.45`$. Since no complete topological map of the www is available, we limited our study to a subset of the web containing $`325,729`$ nodes and $`1,469,680`$ links ($`k=4.59`$). Despite the directedness of the links, the response of the system is similar to the undirected networks we investigated earlier: after a slight initial increase, $`d`$ remains constant in the case of random failures, while it increases for attacks (see Fig. 2c). The network survives as a large cluster under high rates of failure, but the behavior of $`s`$ indicates that under attack the system abruptly falls apart at $`f_c^w=0.067`$ (Fig. 3c).

In summary, we find that scale-free networks display a surprisingly high degree of tolerance against random failures, a property not shared by their exponential counterparts. This robustness is probably the basis of the error tolerance of many complex systems, ranging from cells to distributed communication systems. It also explains why, despite frequent router problems, we rarely experience global network outages or, despite the temporary unavailability of many webpages, our ability to surf and locate information on the web is unaffected. However, the error tolerance comes at the expense of attack survivability: the diameter of these networks increases rapidly and they break into many isolated fragments when the most connected nodes are targeted. Such decreased attack survivability is useful for drug design, but it is less encouraging for communication systems, such as the Internet or the www. While the general wisdom is that attacks on networks with distributed resource management are less successful, our results indicate that the topological weaknesses of the current communication networks, rooted in their inhomogeneous connectivity distribution, have serious effects on their attack survivability, that could be exploited by those seeking to damage these systems.