CosmoCoffee

 FAQ   Search  SmartFeed   Memberlist    Register Profile   Log in Arxiv New Filter | Bookmarks & clubs | Arxiv ref/author:

 CAMB & Two-point correlation function: weird behaviour
Author Message
Guido Walter Pettinari

Joined: 20 May 2008
Posts: 7
Affiliation: University of Bologna

 Posted: February 17 2009 Hi all! I am trying to obtain the two-point correlation function ξ(r) from CAMB's output power spectrum using the following relation $\xi (r)=\frac{1}{2\pi ^2}\int_{0}^{\infty}k^2 P(k) \frac{\sin(k\cdot r)}{kr} dk$. I need to calculate ξ up to r = 106. I used several methods to solve the integral, including the dftint routine from Numerical Recipes, Romberg quadrature and some simpler algorithms (e.g. trapezoidal method). The result is always the same: I get an anticorrelation ξ(r) < 0 for r > 120h − 1Mpc and some very noisy behaviour when r becomes bigger than ~104h − 1Mpc. By noisy I mean that the correlation function goes quickly to zero (it has r - 4 strength) but by doing so it oscillates wildly. I first tought that the oscillation could be due to the finite extent of CAMB's P(k), so I extrapolated the power spectrum in this way: $P(k) \propto k$ when k < 10 − 6 $P(k) \propto k^{-3}$ when k > 500 In this way I can do the integration analytically on the tails. The resulting ξ(r) is a little less noisy, but the problem persists. I will be grateful to anybody who could answer the following questions: 1) Am I doing some terrible mistake I am not aware of? 2) Does the anticorrelation at big r's physically make sense? 2) Most importantly, are you aware of an algorithm specifically designed to obtain the two point correlation function from a model power spectrum? (and hopefully eliminate the noisiness?) Thank you very much. Cheers, Guido P. S. I noticed that CAMB's power spectrum shows a $k^{-2.8\div-3.2}$ dependence at large k's. The exponent varies depending on kmax, i.e. I get −3.1 setting kmax=500 and −2.9 with kmax=10000. I realize that going to such scales is a physical nonsense; nonetheless, shouldn't P(k) show a k − 3 dependency at large k's?
Patrick McDonald

Joined: 06 Nov 2004
Posts: 17
Affiliation: CITA

 Posted: February 17 2009 Andrew Hamilton's FFTLog code is generally good for this kind of thing.
Antony Lewis

Joined: 23 Sep 2004
Posts: 675
Affiliation: University of Sussex

 Posted: February 17 2009 CAMB includes baryons, so you see the effect of them on very small scales (though not at high precision). For direct evaluation (not tried FFTLog, but sounds like a good idea..): In the tail when the sine is oscillating fast compared to variations in P you can also do a series exapansion of P, then integrate analytically over the oscillations to give a result in terms of the derivatives of P that is much more numerically stable. Other tricks include writing the integral in terms of a sine offset by half a phase, then integrating the sum of that and the original, so the integrand is now something small.
Guido Walter Pettinari

Joined: 20 May 2008
Posts: 7
Affiliation: University of Bologna

 Posted: February 17 2009 Hi Patrick, Antony, thank you very much for the quick answers. I will check FFTLog out and post my results. I have a question concerning the first trick you suggested me, Antony. Is it to some extent equivalent to do a spline interpolation of P(k), extracts the coefficients of the piecewise polynomial, and then sum the contributes of the now analitycally integrable function over all the intervals? Schematically: CAMB $\longrightarrow$ P(ki), i=1, 2, ..., N $\longrightarrow$ Spline interpolation $\longrightarrow$ I have got N−1 polynomials it this form: Pi(k) with k belonging to $[k_{i}, \, k_{i+1}]$ $\longrightarrow$ analytical integration of $\sum_{i=0}^{N-1}\int_{k_i}^{k_{i+1}}k^2P(k)\frac{\sin{kr}}{kr}dk$. The results are the same as when using Romberg quadrature or dftint, but somewhat quicker. However, the "ringing" effect is still there. Cheers, Guido
Anze Slosar

Joined: 24 Sep 2004
Posts: 198
Affiliation: Brookhaven National Laboratory

 Posted: February 18 2009 Correlation funtion is subject to integral constraint, so if you want to have positive correlation at small scales, then you MUST have negative correlation at large scales... The problems you're seeing are connected with oscillating integrand as other people have suggested. One simple solution is to use NAG's D01AUF which can integrate oscillating functions by using magic.
Guido Walter Pettinari

Joined: 20 May 2008
Posts: 7
Affiliation: University of Bologna

 Posted: February 19 2009 Hi all! I tried the FFTLog algorithm. It is very fast, but the obtained correlation function shows an even more noisy behaviour for r>104 h−1 Mpc. In the attached plots you can see it very clearly. The three curves represent three different methods I adopted: Hamilton's fftlog, Numerical Recipes' dftint and the splines semi-analytical approach I explained in my last post in this board. The first graph spans a range from r=10−2 to r=104 h−1 Mpc, while the second one is the zoom between 103 and 104 . I re-scaled the curves to make them fit in one graph; being a loglog graph, I had to plot xi * (−1) for r>120 h−1 Mpc (the curves change color at that point). Is this oscillating behaviour at small scales entirely due to the baryons? @Anze: Thank you very much for your suggestion. Alas I have no access to the NAG library; instead I used the dftint routine for Fourier sine integrals from Numerical Recipes. By the way, why is the correlation function subject to integral constraint? Is it some sort of normalization? Thank you all for your attention. Cheers, Guido
Chris Orban

Joined: 11 Dec 2008
Posts: 2
Affiliation: The Ohio State University

 Posted: June 18 2009 Hi Guido, The noisiness at large r is likely from evaluating the integrand at very high k values where P(k) may not be evaluated perfectly (e.g. k = 50 h /Mpc, well beyond the usual range of interest for cosmologists). The solution is to add an exponential damping factor in your integrand, exp(-k*lambda0), where lambda0 is some constant, say 0.1 Mpc/h. See section 42 of Peebles Large Scale Structure to see what I mean. The only consequence of introducing this cutoff is that xi (r) will be wrong on scales smaller than lambda0. And, as I said, this should make the oscillations go away. Also see section 42 for a discussion of how zero crossings are an expected feature in xi (r). -Chris
Guido Walter Pettinari

Joined: 20 May 2008
Posts: 7
Affiliation: University of Bologna

 Posted: June 21 2009 Thank you very much for your suggestion and for the reference, Chris. I had a look at section 42 of Peebles' book and I read the zero-crossing part. I also found out that the r − 4 behaviour is expected. I implemented the damping exponential function in my code and I obtained a great improvement. The noisiness is still there, but now it starts at $r_n \simeq 10^8 \text{Mpc}/h$ instead of $r_n \simeq 10^4 \text{Mpc}/h$. As damping scales I used 10, 100, 1000 Mpc/h (my input P(k) file arrives up to 15k) but, apart from a slightly modified low r behaviour, the output does not change much. However, I noted that for r > rn the correlation function shows a r − 2 dependence instead of the predicted r − 4. This is easily seen in the attached figure, where I plotted the fit to the noisy part as well. I am convinced that this behaviour is a numerical artifact, possibly ringing http://en.wikipedia.org/wiki/Ringing_artifact), but I do not know how to remove it. Do you have any idea? Cheers, Guido
Chris Orban

Joined: 11 Dec 2008
Posts: 2
Affiliation: The Ohio State University

 Posted: June 22 2009 I need to think about it some more to explain the r−2 behavior in detail. If you really need it to be smooth out there, my only advice is to try different lambda0 values (usually the smaller the better) and to try double precision instead of single precision (i.e. "float" in C programming) if those plots were constructed with single precision. However I'm surprised that you need to calculate xi (r) out that far. That's well beyond the range probed with Large Scale structure surveys and (I think) beyond the range one would need to compare with CMB observations. Do you really need it to be accurate all the way out to 105 Gpc?
Guido Walter Pettinari

Joined: 20 May 2008
Posts: 7
Affiliation: University of Bologna

 Posted: June 22 2009 Hi Chris, you are absolutely right, it does not make much physical sense to go at such high separations. However, I need my correlation function to satisfy the integral constraint $\int_0^\infty r^2 \xi (r) = 0$ as much as possible, since any residual makes the successive step in my project less precise. This is achieved only by having a very high upper integration limit, i.e. a very large rmax. I could solve my problem by just extrapolating ξ(r) with a r−4 curve, but I would rather prefer avoiding this. By the way, I am already using double precision. Actually, I implemented the MPFR library in my code, so that I can attain 128 bit precision, but alas with no difference at all. A last remark: the r − 2 dependence of the correlation function at large r is not peculiar to the exponential damping technique you suggested me since I can reproduce it with a un-dumped input power spectrum as well. The only difference is that with a plain P(k) the r − 2 behaviour starts at lower r's. Moreover, it seems that the starting point of the r − 2 dependence, rn, is inversely proportional to kmax: the more extended the power spectrum, the less noisy the correlation function. Thank you again for your help! Cheers, Guido EDIT: I ran my code with λ0 = 1 and 0.1, thus obtaining the same result: the r − 2 behaviour is still there. In this way the damping is so large that at large k's the damped P(k) is 0 within machine precision.
 Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT + 5 Hours Page 1 of 1

 Jump to: Select a forum Arxiv paper discussion----------------arXiv papers Topic discussion----------------Early UniverseCosmological ModelCosmological Observations  Projects and Resources----------------Computers and softwareTeaching, Papers and PresentationsResearch projectsiCosmo Coming up----------------Job vacanciesConferences and meetings Management----------------CosmoCoffee
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum