CosmoCoffee

[Guido Walter Pettinari]

Hi all!

I am trying to obtain the two-point correlation function ξ(r) from CAMB's output power spectrum using the following relation

.

I need to calculate ξ up to r = 10⁶.
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 ^{− 1}Mpc and some very noisy behaviour when r becomes bigger than ~10⁴h ^{− 1}Mpc. 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:

when k < 10 ^{− 6}
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 dependence at large k's. The exponent varies depending on k_max, i.e. I get −3.1 setting k_max=500 and −2.9 with k_max=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]

Andrew Hamilton's FFTLog code is generally good for this kind of thing.

[Antony Lewis]

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]

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 P(k_i), i=1, 2, ..., N Spline interpolation
I have got N−1 polynomials it this form: P_i(k) with k belonging to analytical integration of .

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]

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]

Hi all!

I tried the FFTLog algorithm. It is very fast, but the obtained correlation function shows an even more noisy behaviour for r>10⁴ h⁻¹ 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⁻² to r=10⁴ h⁻¹ Mpc, while the second one is the zoom between 10³ and 10⁴ . 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⁻¹ 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]

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]

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 instead of . 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 > r_n 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]

I need to think about it some more to explain the r⁻² 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 10⁵ Gpc?

[Guido Walter Pettinari]

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 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 r_max. I could solve my problem by just extrapolating ξ(r) with a r⁻⁴ 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, r_n, is inversely proportional to k_max: 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 λ₀ = 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.