[astroph/0503277] Speeding Up Cosmological Boltzmann Codes
Authors:  Michael Doran 
Abstract:  We introduce a novel strategy for cosmological Boltzmann codes leading to an increase in speed by a factor of \sim 30 for small scale Fourier modes. We (re)investigate the tight coupling approximation and obtain analytic formulae reaching up to the octupoles of photon intensity and polarization. This leads to accurate results reaching optimal precision, while still being simple. Damping rapid oscillations of small scale modes at later times, we simplify the integration of cosmological perturbations. We obtain analytic expressions for the photon density contrast and velocity as well as an estimate of the quadrupole from after last scattering until today. These analytic formulae hold well during reionization and are in fact negligible for realistic cosmological scenarios. However, they do extend the validity of our approach to models with very large optical depth to the last scattering surface. 
[PDF] [PS] [BibTex] [Bookmark] 

 Posts: 1618
 Joined: September 23 2004
 Affiliation: University of Sussex
 Contact:
[astroph/0503277] Speeding Up Cosmological Boltzmann Codes
In this paper the author considers some useful ways to speed up the calculation of linear CMB power spectra and the matter power spectrum.
It first looks at tight coupling, suggesting the possibly useful idea of switching off tightcoupling at different times for different terms. It also considers order [tex](k\tau_c)^2[/tex] terms for the photon modes and claim a significant improvement. Comment: CMBFAST is actually setting the anistropic stress (and E quadrupole) to zero during tight coupling  I suspect most of the improvement comes from including these first order terms rather than the higher terms (CAMB only includes first order terms  the (fairly standard) first order equations being used are publicly documented in the notes linked from the web page).
The idea for 'curing' rapid oscillations is a good one, and something that I've been wanting to get round to trying for a while. Schemes like this should give a significant time saving for the matter power spectrum calculation. For the CMBonly calculation this is less of a problem as the evolution is usually stopped just after recombination (because the late time sources are negligible).
I find Fig. 6 rather odd: obviously there's no actual gauge ambiguity in matter density observables. The 'synchronous gauge' P(k) (i.e. the density power spectrum in the rest frame of the linear CDM perturbations) can be calculated unambiguously in any other gauge. Perhaps he is actually (rather perversely) plotting different P(k)'s in the two cases?
PS. the main ref for CAMB is astroph/9911177, not the CosmoMC ref astroph/0205436.
It first looks at tight coupling, suggesting the possibly useful idea of switching off tightcoupling at different times for different terms. It also considers order [tex](k\tau_c)^2[/tex] terms for the photon modes and claim a significant improvement. Comment: CMBFAST is actually setting the anistropic stress (and E quadrupole) to zero during tight coupling  I suspect most of the improvement comes from including these first order terms rather than the higher terms (CAMB only includes first order terms  the (fairly standard) first order equations being used are publicly documented in the notes linked from the web page).
The idea for 'curing' rapid oscillations is a good one, and something that I've been wanting to get round to trying for a while. Schemes like this should give a significant time saving for the matter power spectrum calculation. For the CMBonly calculation this is less of a problem as the evolution is usually stopped just after recombination (because the late time sources are negligible).
I find Fig. 6 rather odd: obviously there's no actual gauge ambiguity in matter density observables. The 'synchronous gauge' P(k) (i.e. the density power spectrum in the rest frame of the linear CDM perturbations) can be calculated unambiguously in any other gauge. Perhaps he is actually (rather perversely) plotting different P(k)'s in the two cases?
PS. the main ref for CAMB is astroph/9911177, not the CosmoMC ref astroph/0205436.

 Posts: 41
 Joined: November 22 2004
 Affiliation: ITP Heidelberg
 Contact:
[astroph/0503277] Speeding Up Cosmological Boltzmann Codes
Hi Antony,
as far as tight coupling is concerned, it is important to have a
good estimate of the shear. Not so much for the slip, but much more
so for the evolution equations for the baryon and photon velocities.
Now anything beyond the CMBFAST approximation shear=0 is of
course a great improvment. My approach was to gain as
much as possible while keeping the algebra and numerics simple.
Including the Octupole and the time derivative of the shear to get an
improved estimate of the shear is algebraically and numerically inexpensive,
which is one reason why I did it. The other reason is that it truely
improves on the precision of the analytic result. Maybe I should add a plot
to the paper comparing the two. The point is that the result of the
shear including the higher order corrections is very precise as far as an
expansion in [tex](k\tau_c)[/tex] is concerned and additionally
improves on the assumption that the time derivative of the shear is zero
(in the analytic expression for the shear).
As far as curing the rapid oscillations is concerned, I was rather surprised
that it worked. For several days, I was playing with all kinds of ideas to
get rid of the oscillations (like using different variables, projecting out
some things, averaging, etc.) and in the end, I thought: why not get
rid of the evolution alltogether ? As you might have anticipated, it
works. For high precision and considerable optical depth, though, it
was necessary to get a grip on [tex]\delta_\gamma[/tex], and an
upper bound for the shear which is why I discuss the evolution during
reionization.
You are of course perfectly right that this doesn't really speed up
the lowl cmb calculations. Personally, however, I most frequently
use the CMB together with LSS and for this it helps. The true benefit,
however lies in the speed improvements for modes [tex]k > 0.1 Mpc^{1}[/tex],
from then on, it really becomes incredibly fast compared to the
standard strategy
As far as Figure 6 is concerned, I have to reject the idea that I did something
perverse :) I did plot the power spectrum of the synchronous gauge density
contrast together with the power spectrum in the gauge invariant
[tex]D_g = \delta_{longit} + 3(1+w)\Phi[/tex]. That's what I mean by "gauge
ambiguities". Had I plotted the synchronous gauge density contrast infered
from the gauge invariant one, the two curves would of course fall on top of
each other. I just like contrasting the two, because I find it nice to see when
the two "gauges" start to coincide. That's personal taste, of course. But
hopefully not perverse :)
Finally, I'm sorry for citing the wrong paper. I just had our first conversation at
Moriond (was it 2002 ?) in mind, where you said that you didn't write a paper
on CAMB really (my mistake, I must have misunderstood it). As I wanted to
give a reference, I took the cosmomc paper. But I'll add the reference now
that I know that's the official one :)
Cheers, Michael
as far as tight coupling is concerned, it is important to have a
good estimate of the shear. Not so much for the slip, but much more
so for the evolution equations for the baryon and photon velocities.
Now anything beyond the CMBFAST approximation shear=0 is of
course a great improvment. My approach was to gain as
much as possible while keeping the algebra and numerics simple.
Including the Octupole and the time derivative of the shear to get an
improved estimate of the shear is algebraically and numerically inexpensive,
which is one reason why I did it. The other reason is that it truely
improves on the precision of the analytic result. Maybe I should add a plot
to the paper comparing the two. The point is that the result of the
shear including the higher order corrections is very precise as far as an
expansion in [tex](k\tau_c)[/tex] is concerned and additionally
improves on the assumption that the time derivative of the shear is zero
(in the analytic expression for the shear).
As far as curing the rapid oscillations is concerned, I was rather surprised
that it worked. For several days, I was playing with all kinds of ideas to
get rid of the oscillations (like using different variables, projecting out
some things, averaging, etc.) and in the end, I thought: why not get
rid of the evolution alltogether ? As you might have anticipated, it
works. For high precision and considerable optical depth, though, it
was necessary to get a grip on [tex]\delta_\gamma[/tex], and an
upper bound for the shear which is why I discuss the evolution during
reionization.
You are of course perfectly right that this doesn't really speed up
the lowl cmb calculations. Personally, however, I most frequently
use the CMB together with LSS and for this it helps. The true benefit,
however lies in the speed improvements for modes [tex]k > 0.1 Mpc^{1}[/tex],
from then on, it really becomes incredibly fast compared to the
standard strategy
As far as Figure 6 is concerned, I have to reject the idea that I did something
perverse :) I did plot the power spectrum of the synchronous gauge density
contrast together with the power spectrum in the gauge invariant
[tex]D_g = \delta_{longit} + 3(1+w)\Phi[/tex]. That's what I mean by "gauge
ambiguities". Had I plotted the synchronous gauge density contrast infered
from the gauge invariant one, the two curves would of course fall on top of
each other. I just like contrasting the two, because I find it nice to see when
the two "gauges" start to coincide. That's personal taste, of course. But
hopefully not perverse :)
Finally, I'm sorry for citing the wrong paper. I just had our first conversation at
Moriond (was it 2002 ?) in mind, where you said that you didn't write a paper
on CAMB really (my mistake, I must have misunderstood it). As I wanted to
give a reference, I took the cosmomc paper. But I'll add the reference now
that I know that's the official one :)
Cheers, Michael

 Posts: 1618
 Joined: September 23 2004
 Affiliation: University of Sussex
 Contact:
Re: [astroph/0503277] Speeding Up Cosmological Boltzmann Co
OK, thanks. I'd be interested to know what precision you get on the matter power spectrum (in consistent gauges!) using the parameters given in the paper for turning off the oscillations?
In case anyone is interested, in the synchronous gauge, during matter domination but before reionization, the nonoscillatory parts of the radiation densities and heat fluxes of adiabatic modes evolve as
[tex] \Delta_\gamma \approx \Delta_\nu \approx \frac{2\kappa S^2}{3k^2} (\rho_c\Delta_c+\rho_b\Delta_b)[/tex]
[tex] q_\gamma \approx q_\nu \approx
\frac{2k}{3}\frac{\Delta_\gamma}{\sqrt{\kappa S^2[\rho_b+\rho_c]/3}}[/tex]
(S=scale factor). In particular the synchronous gauge radiation velocities (3/4 q) are significant and needs to be included, unlike in the Newtonian gauge where they become small. Using these analytic approximations and setting all other terms to zero seems to work reasonably well after matter domination for the synchronous gauge. Note that [tex]\Delta[/tex] for radiation differs by a factor of three from the Newtonian result.
One more question: does the Newtonian gauge code work for isocurvature modes? I thought the usual reason people used synchronous gauge was because there are nasty instabilities in the Newtonian equations for nonadiabatic modes.
In case anyone is interested, in the synchronous gauge, during matter domination but before reionization, the nonoscillatory parts of the radiation densities and heat fluxes of adiabatic modes evolve as
[tex] \Delta_\gamma \approx \Delta_\nu \approx \frac{2\kappa S^2}{3k^2} (\rho_c\Delta_c+\rho_b\Delta_b)[/tex]
[tex] q_\gamma \approx q_\nu \approx
\frac{2k}{3}\frac{\Delta_\gamma}{\sqrt{\kappa S^2[\rho_b+\rho_c]/3}}[/tex]
(S=scale factor). In particular the synchronous gauge radiation velocities (3/4 q) are significant and needs to be included, unlike in the Newtonian gauge where they become small. Using these analytic approximations and setting all other terms to zero seems to work reasonably well after matter domination for the synchronous gauge. Note that [tex]\Delta[/tex] for radiation differs by a factor of three from the Newtonian result.
One more question: does the Newtonian gauge code work for isocurvature modes? I thought the usual reason people used synchronous gauge was because there are nasty instabilities in the Newtonian equations for nonadiabatic modes.

 Posts: 41
 Joined: November 22 2004
 Affiliation: ITP Heidelberg
 Contact:
[astroph/0503277] Speeding Up Cosmological Boltzmann Codes
Hi Antony,
I modified figure 6 to compare the gauge invariant cdm power spectrum using the old gauge invariant implementation to the new strategy in gauge invariant variables. The relative deviation between the two is roughly 0.02 per cent, in other words nothing :)
Concerning instabilities for isocdm, I must say that I didn't see any yet. Whenever I compare the gauges (I do that quite often to cross check), they coincide. The exception to the rule is that for very large k, i.e. very small modes (70 Mpc^1), the gauge invariant implementation seems to become unstable. I didn't have time to look into the reason, though and didn't notice before the speedup paper, because I never went to such high kvalues. Just to clarify: the code does not run in longitudinal gauge but uses gauge invariant variables.
As said, I would have to spend more time investigating, why the CDM spectrum becomes unstable for k > 70 Mpc^1. If anyone has a good explanation (my guess are the potentials: density contrasts>potentials>poisson equation>gravitational infall for densities), please go ahead and post. Or better still: maybe someone sees a way to circumvent the numerical instability. Most propably it's cancellation at work.
To be honest, I use synchronous gauge for practically all my numerical work. But I like the concept of gauge invariance and in particular newtonian gauge for analytic calculations.
I modified figure 6 to compare the gauge invariant cdm power spectrum using the old gauge invariant implementation to the new strategy in gauge invariant variables. The relative deviation between the two is roughly 0.02 per cent, in other words nothing :)
Concerning instabilities for isocdm, I must say that I didn't see any yet. Whenever I compare the gauges (I do that quite often to cross check), they coincide. The exception to the rule is that for very large k, i.e. very small modes (70 Mpc^1), the gauge invariant implementation seems to become unstable. I didn't have time to look into the reason, though and didn't notice before the speedup paper, because I never went to such high kvalues. Just to clarify: the code does not run in longitudinal gauge but uses gauge invariant variables.
As said, I would have to spend more time investigating, why the CDM spectrum becomes unstable for k > 70 Mpc^1. If anyone has a good explanation (my guess are the potentials: density contrasts>potentials>poisson equation>gravitational infall for densities), please go ahead and post. Or better still: maybe someone sees a way to circumvent the numerical instability. Most propably it's cancellation at work.
To be honest, I use synchronous gauge for practically all my numerical work. But I like the concept of gauge invariance and in particular newtonian gauge for analytic calculations.