Abstract
We investigate the possibility of reconstructing the initial spectrum of density fluctuations from the cosmic microwave background (CMB) anisotropy. As a first step toward this program, we consider a spatially flat, CDM dominated universe. In this case, it is shown that, with a good accuracy, the initial spectrum satisfies a first order differential equation with the source determined by the CMB angular correlation function. The equation is found to contain singularities arising from zeros of the acoustic oscillations in the transfer functions. Nevertheless, we find these singularities are not fatal, and the equation can be solved nicely. We test our method by considering simple analytic forms for the transfer functions. We find the initial spectrum is reproduced within 5% accuracy even for a spectrum that has a sharp spike.
pacs:
PACS: 95.30.k; 98.80.EsOsaka University Theoretical Astrophysics
November 2001 OUTAP168
Cosmic Inversion
—Reconstructing primordial spectrum from CMB anisotropy—
Makoto Matsumiya,^{1}^{1}1Email: Misao Sasaki^{2}^{2}2Email: and Jun’ichi Yokoyama^{3}^{3}3Email: Department of Earth and Space Science, Graduate School of Science,
Osaka University, Toyonaka 5600043, Japan
I Introduction
Determination of the primordial spectrum of the density fluctuations is one of the most important issues in modern cosmology. The cosmic microwave background (CMB) anisotropy provides us with a great deal of information of the primordial fluctuations, and it is considered to be a powerful tool for studying the early universe [1]. This is because the physical processes that determine the CMB anisotropy are described by linear perturbation theory, and they have been well understood [2, 3].
Although the recent anisotropy observations are consistent with a flat universe with a scaleinvariant initial power spectrum, this does not exclude the possibilities of other models. In particular, we do not know how much the shape of the initial power spectrum is constrained. In most of previous investigations, when cosmological model parameters are estimated from the observational data by likelihood analysis, the initial spectrum is assumed to have a powerlaw shape [4]. It is true that a conventional slowroll inflation model [5], that has now become a ‘standard model’, gives a powerlaw spectrum which is almost scaleinvariant[6]. However, when analyzing the observed CMB anisotropy, it is much more desirable, and probably much healthier, to constrain the initial spectrum solely by observed data without any theoretical prejudices. For example, even within the context of inflationary cosmology, a variety of generation mechanisms for nonscaleinvariant perturbations have been proposed [7]. In this connection, recently several authors have discussed extraction of nonpower law features from the CMB observations [8], where the initial spectrum is allowed to have a piecewise powerlaw shape.
In this paper, we approach this issue in an entirely different way, namely, by formulating an inverse problem as faithful as possible. Such an approach will be eventually needed if we seriously want to constrain the initial power spectrum solely from observations of the CMB anisotropy. An approach to this inversion problem has been discussed recently [9].
As a first step, we consider a simple situation in which the transfer functions that relate the input power spectrum of the gravitational potential to the output CMB angular correlation function are given analytically. This is certainly a toy model. However, it has almost all the essential features a realistic model would have. In particular, unlike [9], our model takes account of not only the SachsWolfe (SW) effect[10] but also the Doppler effect . The latter, which gives rise to zero points in the transfer functions, is the main cause of the difficulty in this inversion problem.
The advantage of adopting this simple situation is that our method of inversion which we shall develop below may be easily tested at various stages of calculations. Since our primary concern here is to formulate the inversion problem, we fix the cosmological parameters and do not study the dependence of upon them.
The paper is organized as follows. In Sec. II, we give the basic equations that relate the primordial spectrum with the angular correlation function . In Sec. III, under some reasonable assumptions, we derive a differential equation for and develop a method to solve it. Then we test our method by applying it to several spectral shapes. We find our method is applicable even in the case of a spectrum with a sharp spike.
Ii Basic Equations
The angular correlation function of the CMB, is defined as
(1) 
where is the temperature fluctuation in the direction and the average is taken over all angles and all spatial positions. is expanded in the Legendre polynomials and related to the angular power spectrum, , as
(2) 
or are the fundamental observables of the CMB anisotropy. Our purpose is to reconstruct the initial power spectrum from the above observables. We calculate in the Newtonian gauge and follow the notation of [3].
Each Fourier mode of temperature perturbations, , obeys the Boltzmann equation [11],
(3) 
where the dot denotes a derivative with respect to the conformal time , is the comoving wave number, , is the differential Thomson optical depth and is the bulk velocity of baryons. and are the gaugeinvariant Newtonian potential and spatial curvature perturbation, respectively[11]. Here, we decompose into the multipole moments,
(4) 
where is the th multipole moment of . Integrating Eq. (3), we obtain
(5) 
where is the conformal time today and is the visibility function given by
(6) 
We have neglected the quadrupole term on the right hand side of Eq. (3) since its contribution is negligible in the tight coupling approximation. The visibility function has a sharp peak around the last scattering time so that we assume that recombination occurs instantaneously at . Then, the multipole moments of each mode is approximately given by
(7)  
where is a conformal distance from the present epoch to the last scattering surface (LSS). Equation (7) shows that there are four sources for the temperature anisotropy, namely, the intrinsic temperature variation, the SachsWolfe (SW) effect which is caused by the static gravitational potential at the last scattering surface and is dominant on large scales, the Doppler effect due to the fluid bulk velocity, and the integrated SachsWolfe (ISW) effect due to the time variation of the metric perturbation after decoupling. If the recombination occurs when the universe is matterdominated, and remain constant in time after decoupling until a stage at which either the curvature term or the cosmological constant term becomes significant. In particular, in the Einsteinde Sitter universe, they remain constant until today, hence the ISW effect is absent.
Conventionally, is expressed as
(8) 
From Eqs. (7) and (8), we find
(9) 
The angular correlation function is calculated from Eq. (9). Here, we define by
(10) 
This is the spatial distance between two points on the last scattering surface which are observed with the angular separation . Since the thickness of the LSS is neglected in Eq. (7), there is a onetoone correspondence between the observed temperature anisotropy and the perturbation variables on the LSS. Using the relation,
(11) 
the angular correlation function is given by
(12)  
In the linear perturbation theory, different modes do not couple with each other but evolve independently. The initial condition for each mode is directly reflected in the perturbation variables through transfer functions. The variables, , and , in the integrand of Eq. (12) are therefore linearly related to the initial condition and we can generally express the relation between and the initial spectrum as
(13) 
where we normalize the initial condition in terms of the curvature perturbation and define . In this case, introducing the transfer functions and defined by
(14) 
can be written as
(15) 
where are given by
(16) 
Our purpose is to solve the integral equation (13) for with the kernel given by Eq. (15). When there is only the monopole term, i.e., , Eq. (13) takes the familiar form,
(17) 
In this case, using the Fourier sine formula, we obtain
(18) 
It must be noted that has singularities on small scales since and are oscillatory functions reflecting the acoustic oscillations of the density perturbations. These singularities are inevitable as long as we take this approach and the onetoone correspondence between and breaks down at the singularities. However, we shall find in the next section that there is a way to resolve this difficulty.
Iii Inversion Method
In the above discussion, we have tacitly assumed that runs from zero to infinity. In reality, however, is bounded in the finite range . Furthermore, it is observationally impossible to determine on large scales due to the statistical ambiguity, i.e., the cosmic variance. However, the scales which we are interested in are and it is expected that modes with have little effect on these scales. We therefore neglect the terms proportional to and in Eq. (12). In this limit, we have
(19) 
where we have replaced and with and , respectively, for notational simplicity. Multiplying by and integrating by parts twice, we may reexpress Eq. (19) as
(20) 
where we have assumed that the boundary terms vanish and the integral converges. Employing the Fourier sine formula, we obtain
(21) 
This is a secondorder ordinary differential equation in which the source term carries the information of . It is of course necessary to fix the boundary conditions for in order to solve Eq. (21). However, we have no means to determine these conditions. The only knowledge we have is that the boundary conditions at and are restricted from the convergence of the integral in Eq. (19).
In order to avoid this difficulty, we consider the reduction of the order of the differential equation. If we take some appropriate combination of and its derivative, we can derive a differential equation of lower order. The simplest combination is given by
(22) 
Note that the power of in the integrand is reduced. No linear combination can reduce the order of the differential equation lower than . Integrating by parts, we obtain
(23) 
This is a first order differential equation. We discuss a method for solving this equation in the rest of this section.
We now give the expressions for and explicitly in order to examine the properties of Eq. (23), and to verify if thus obtained solution correctly reproduces the original spectrum. For a given model, and are determined by the coupled Einsteinfluid equations which are to be solved numerically. However, in order to understand the property of Eq. (23) and to establish a method which we can apply to general cases, we start with a toy model in which and are given analytically.
Before recombination, photons are tightly coupled with baryons through the Thomson scattering. With the tight coupling approximation, we can expand the Boltzmann equation for photons together with the continuity and Euler equations for baryons in the Thomson scattering time [2]. To the leading order, the dynamics of the photonbaryon fluid is described by the following simple equation for each Fourier mode [3],
(24) 
where is the sound speed of the photonbaryon fluid, , and is the driving term due to the gravitational potential,
(25) 
The dipole moment is obtained from the continuity equation,
(26) 
Note that Eqs. (24) and (26) are applicable not only to the Einsteinde Sitter universe but also to other cosmological models.
Here we solve Eqs. (24) and (26) under the following assumptions:

The gravitational potential fluctuation and the spatial curvature fluctuation are always independent of time, i.e., , and .

The baryon denisty is negligible, i.e., .

The perturbation is adiabatic.
The assumption (i) is valid only on large and intermediate scales, but fails on small scales (i.e., on scales much smaller than the Hubble horizon scale at decoupling) since the acoustic oscillations tend to destroy the gravitational potential on such scales. The assumption (ii) is, of course, quite unrealistic, while (iii) is a natural consequence of inflationary cosmology. However, these assumptions do not alter the essential features of the CMB anisotropy spectrum. The difficulty of solving Eq. (23) is caused by the acoustic oscillations of temperature fluctuations which give rise to zero points in the transfer functions and , and this is a common feature of the CMB anisotropy regardless of cosmological models. Furthermore, we do not take account of the diffusion damping [12] which becomes significant on small scales, since this effect do not significantly change the structure of Eq. (23) either.
We solve Eqs. (24) and (26) in this toy model. The solutions are given by
(27) 
where and we take the adiabatic initial condition as and . From Eq. (23), we obtain
(28) 
where and
(29) 
Now let us describe our method. When the data of are given, is computed from the Fourier sine transform of it. In actual situations, is to be obtained by observation. However, here we use obtained from Eq. (12) by giving an original spectrum by hand. Then we solve Eq. (28) and compare the solution with the original spectrum. Note that we do not use Eq. (19) for the evaluation of since it is an approximate formula obtained by taking the limit . This is because our purpose is to examine the accuracy of our method which uses several approximations. For the distance to the LSS, we choose .
The theoretical expression for includes contributions from of all ; . In reality, the monopole and dipole contributions of the anisotropy cannot be observed and the higher multipole moments are limited by the angular resolution of observation. Therefore, we subtract these contributions from and define the following function instead of Eq. (2),
(30) 
where is an upper limit of and we take in the following calculation. We denote the corresponding by . Since the , terms contribute to mainly at large , the source term is not expected to be modified significantly by neglecting these terms. As for the contributions of , they are negligible since decreases exponentially for large . In Fig. 1, we show and for a power law spectrum with a damping factor, . In this case, can be calculated analytically and we can estimate the effect of using the finite data set of . We also plot which is given by cutting off the high multipole contributions of from . The relative error between and the analytic is below 2 %, which justifies the neglect of the high multipole contributions. On the other hand, differs from substantially. However, since we are interested in scales much smaller than , this subtraction of modes is expected to be (and in fact found to be) harmless.
For a given spectrum, the exact can be evaluated by simply substituting the original into the lefthand side of Eq. (28). Since we are interested in the inversion problem, we need to know how accurately the source term can be calculated from . In particular, we must reproduce the exact one with a good accuracy on scales . We therefore examine the accuracy of the calculated by comparing it with the exact one.
Before we make this comparison, however, there is one more problem to be resolved. It is the problem of the range of definition of . The Fourier transform of cannot be performed in the exact sense since given by Eq. (10) is not defined for . To avoid this difficulty, we cut off by a smooth function which damps out steeply at some scale which is much smaller than but large enough to reproduce the spectrum over a sufficiently wide range. Although this lowers the amplitude of large modes, the function obtained by the Fourier transform has little deviation from the exact on scales of our interest. Both the exact and calculated are shown in Fig. 2. Since we take , the region where the calculated source term matches the exact one is below .
As mentioned before, Eq. (28) has singularities at which cause difficulties when we solve it numerically. We can, however, determine the values of at these singularities if we assume that the derivative of is finite. Then the first term on the lefthand side of Eq. (28) vanishes at , and the values of at these singularities are given by
(31) 
Once these values are given, we can solve Eq. (28) by expanding it around the singularities. We search the true solution which connect the adjacent singularities using the shooting method. We solve the equation until the 5th singularity, , for the following two cases of the original spectra:
(32)  
(33) 
The first one is a double powerlaw spectrum and the second is a single powerlaw spectrum with a spiky structure, either with a peak () or a dip (). The results for the choice of the parameters as , , , and are shown in Fig. 3. We find our method reproduces the original spectra with a good accuracy. In particular, even if the spectrum has a sharp peak or a dip, we can resolve such a local structure using this method. The numerical solution diverges as it approaches the singularities (indicated by the triangles), but the relative error except for the regions close to the singularities is below 5 %.
For comparison, in Fig. 4, we also show the corresponding CMB angular power spectra calculated by using the transfer functions given by Eq. (27). In the case of the spectrum with a peak, one can clearly see an enhancement of at , though the peak height is suppressed by a factor of as compared to the original peak in the primordial spectrum. As for the spectrum with a dip, the structure is even more suppressed and leaves only a small imprint in the angular spectrum. In a realistic case, a small structure may be present in the CMB spectrum as a consequence of an irregular feature in the primordial spectrum. The above result indicates that our method is capable of reproducing such a feature from the CMB spectrum.
Finally, it is worthwhile to comment that the presence of the singularities in the differential equation (23) may be regarded as an advantage, since the values of at the singularities can be estimated without solving the differential equation. In particular, if there is a good reason to believe that the spectrum should be a smoothly varying function, a qualitative feature of the spectrum can be obtained at once. For example, in the case of the double powerlaw spectrum (32), one can see that the original spectrum can be approximately recovered by simply interpolating between the adjacent triangles shown in Fig. 3.
Iv Conclusion
We have considered the problem of reconstructing the initial power spectrum of metric perturbation from data. As a first step, we have investigated a simple case, namely, the Einsteinde Sitter universe with negligible baryons and negligible thickness of the LSS. In this toy model, the observed temperature fluctuations are represented only by the perturbation variables on the LSS and the ISW effect is absent. Then, the relation between the initial spectrum and the angular correlation function is expressed in terms of an integral equation. We have shown that this equation can be transformed to a second order differential equation for . We have found that by forming an appropriate linear combination of the angular correlation function and its derivative, the order of the differential equation can be reduced to the first order. The resulting equation is found to have singularities that come from the acoustic oscillations of photons, hence their presence is inevitable in any cosmological models not restricted to our simple model. Fortunately, however, the presence of these singularities turns out to be not only harmless but rather advantageous. At the singularities the coefficient of vanishes, and the values of there can be obtained without solving the differential equation. By plotting the values of , one can obtain a rough estimate of the behavior of the spectrum. Then, to recover the precise features of the spectrum, the rest part of can be solved with its values at the singularities as the initial values. We have found our method can reproduce the original spectrum with a good accuracy even for a spectrum with a sharp, spiky structure.
The method presented here is applicable only to the Einsteinde Sitter universe in which the ISW effect is negligible. The ISW effect gives an important contribution even in a flat universe model if the cosmological constant is present. Such a model, called the CDM model is preferred by recent observations[13]. Our next step is to include the ISW effect in our method which is currently under study.
Acknowledgements.
This work was supported in part by the Monbukagakusho/JSPS GrantinAid for Scientific Research Nos. 12640269(MS) and 13640285(JY), by the Monbukagakusho/JSPS GrantinAid for Priority Area: Supersymmetry and Unified Theory of Elementary Particles (No. 707)(JY), and by the Yamada Science Foundation.References
 [1] M. J. White, D. Scott and J. Silk, Ann. Rev. Astron. Astrophys. 32, 319 (1994); W. Hu, N. Sugiyama and J. Silk, Nature 386, 37 (1997).
 [2] P. J. E. Peebles and J. T. Yu, Astrophys. J. 162, 815 (1970).
 [3] W. Hu and N. Sugiyama, Astrophys. J. 444, 489 (1995); W. Hu and N. Sugiyama, Phys. Rev. D 51, 2599 (1995).
 [4] J. R. Bond et al. Phys. Rev. Lett. 72, 13 (1994); L. Knox, Phys. Rev. D 52, 4307 (1995); G. Jungman, M. Kamionkowski, A. Kosowski and D. N. Spergel, Phys. Rev. D 54, 1332 (1995); J. R. Bond, G. Efstathiou and M. Tegmark, Mon. Not. R. Astron. Soc. 291, L33 (1997); D. J. Eisenstein, W. Hu and M. Tegmark, Astrophys. J. 518, 2 (1999).
 [5] For a review of inflation, see, e.g., A. D. Linde, Particle Physics and Inflationary Cosmology (Harwood, Chur, Switzerland, 1990); K. A. Olive, Phys. Rep. 190, 181 (1990); D. H. Lyth and A. Riotto, Phys. Rep. 314, 1 (1999).
 [6] S. W. Hawking, Phys. Lett. 115B, 295 (1982); A. A. Starobinsky, ibid 117B, 175 (1982); A. H. Guth and SY. Pi, Phys. Rev. Lett. 49, 1110 (1982).
 [7] A. D. Linde, Phys. Lett. 158B, 375 (1985); L. A. Kofman and A. D. Linde, Nucl. Phys. B282, 555 (1987); J. Silk and M. S. Turner, Phys. Rev. D 35, 419 (1987); H. M. Hodges and G. R. Blumenthal, Phys. Rev. D. 42, 3329 (1990); J. Yokoyama and Y. Suto, Astrophys. J. 379, 427 (1991); M. Sasaki and J. Yokoyama, Phys. Rev. 744, 970 (1991); A. A. Starobinsky, JETP Lett. 55, 489 (1992); [Pisma Zh. Eksp. Teor. Fiz. 55 (1992) 477]; J. Yokoyama, Astron. Astrophys. 318, 673 (1997); J. Yokoyama, Phys. Rev. D 58, 083510 (1998);59, 107303 (1999); S. M. Leach and A. R. Liddle, Phys. Rev. D 63, 043508 (2001) [arXiv:astroph/0010082]; S. M. Leach, M. Sasaki, D. Wands and A. R. Liddle, Phys. Rev. D 64, 023512 (2001) [arXiv:astroph/0101406].
 [8] T. Souradeep et. al., astroph/9802262; Y. Wang, D. N. Spergel and M. A. Strauss, astroph/9812291; Y. Wang, D. N. Spergel and M. A. Strauss, Astrophys. J. 510, 20 (1999); S. Hannestad, Phys. Rev. D 63, 043009 (2001).
 [9] A. Berera and P. A. Martin, Inverse Problems 15, 1393 (1999).
 [10] R. K. Sachs and A. M. Wolfe, Astrophys. J. 147, 73 (1967).
 [11] H. Kodama and M. Sasaki, Prog. Theor. Phys. Supp. 78, 1 (1984); H. Kodama and M. Sasaki, Int. J. Mod. Phys. A 1, 265 (1986); H. Kodama and M. Sasaki, Int. J. Mod. Phys. A 2, 491 (1987).
 [12] J. Silk, Astrophys. J., 151 459 (1968).
 [13] A. G. Riess et al. [Supernova Search Team Collaboration], Astron. J. 116, 1009 (1998) [arXiv:astroph/9805201]; S. Perlmutter et al. [Supernova Cosmology Project Collaboration], Astrophys. J. 517, 565 (1999) [arXiv:astroph/9812133].