# Cosmological perturbations without the Boltzmann hierarchy

Marc Kamionkowski  Department of Physics and Astronomy, Johns Hopkins
University, 3400 N. Charles St., Baltimore, MD 21218

###### Abstract

Calculations of the evolution of cosmological perturbations generally involve
solution of a large number of coupled differential equations to describe the
evolution of the multipole moments of the distribution of photon intensities
and polarization. However, this “Boltzmann hierarchy” communicates with the
rest of the system of equations for the other perturbation variables only
through the photon-intensity quadrupole moment. Here I develop an alternative
formulation wherein this photon-intensity quadrupole is obtained via solution
of two coupled integral equations—one for the intensity quadrupole and another
for the linear-polarization quadrupole—rather than the full Boltzmann
hierarchy. This alternative method of calculation provides some physical
insight and a cross-check for the traditional approach. I describe a simple
and efficient iterative numerical solution that converges fairly quickly. I
surmise that this may allow current state-of-the-art cosmological-perturbation
codes to be accelerated.

##  I Introduction

Linear-theory calculations of the evolution of primordial density
perturbations provide the foundation for the interpretation of cosmic
microwave background and large-scale-structure measurements. They are thus an
essential tool in the construction of our current cosmological model and in
the continuing quest for new cosmological physics.

The calculations, which trace back over 50 years earlywork , involve time
evolution of a set of coupled differential equations 1980s  for the metric
perturbations and for the dark-matter, baryon, neutrino, and photon density
and velocity perturbations. There is also a (nominally infinite) “Boltzmann
hierarchy” of differential equations for the higher moments (quadrupole,
octupole, etc.) of the photon-intensity and photon-polarization and neutrino-
momentum distributions. The photon hierarchies can be truncated at some
maximum multipole moment $l_{\rm max}\simeq 30$ to provide sufficient
precision for the monopole, dipole, and octupole from which the higher-order
moments (which provide the CMB temperature/polarization power spectra) can be
obtained through a line-of-sight integral Seljak:1996is . Higher-order
extensions to the tight-coupling approximation (TCA) CyrRacine:2010bk ;
Blas:2011rf , improved numerical integrators, and novel approximations to
free-streaming relativistic particles Blas:2011rf ) have provided incredible
code acceleration to what is still a fairly complicated numerical calculation.
At present, virtually all work in cosmology now relies on two publicly
available codes, CAMB Lewis:1999bs  and CLASS Blas:2011rf , which combine
speed and precision with model flexibility.

These codes are now extremely efficient and reliable. However, modern
cosmological analyses, which employ Markov chain Monte Carlos to map the
likelihood in a multi-dimensional parameter space, require these codes to be
run repeatedly, thus employing signficant computational resources. It is thus
worthwhile to explore new numerical approaches. New approaches can also often
provide new insights into the physics and may perhaps provide tools that can
be applied to other problems.

It was realized that for primordial tensor perturbations (i.e., gravitational
waves), the Boltzmann hierarchy can be replaced by a small set of integral
equations (IEs) Weinberg:2003ur ; Baskaran:2006qs , an approach used in Refs.
Flauger:2007es ; Pritchard:2004qp  A similar approach was discussed for scalar
perturbations (primordial density perturbations) in Ref. Weinberg:2006hh , but
not implemented numerically.

Here, I re-visit this integral-equation approach for primordial density
perturbations. I discuss simplifications to the equations in Ref.
Weinberg:2006hh  and describe a specific implementation where the Boltzmann
hierarchy for all photon intensity/polarization multipole moments from the
quadrupole ($l=2$) and higher are replaced by two IEs, one for the photon
quadrupole, and another for the polarization quadrupole. I discuss the
numerical solution of these integral equations and how the initial conditions
for the IEs are set from an early-time solution obtained with the TCA. I
describe an iterative algorithm to solve these integral equations
simultaneously with the differential equations for the other perturbation
variables. I show results from two simple numerical codes that are identical
except for the replacement of the Boltzmann hierarchy in the first with the
two integral equations in the second. Numerical experiments with these codes
suggest that this iterative IE algorithm may, with further work, allow current
state-of-the-art codes to be accelerated.

This paper is organized as follows. Section II, presents and discusses the
integral equations. Section III provides the differential equations for the
other perturbation variables (i.e., for neutrinos, dark matter, baryons, and
the metric) and describe how the two integral equations are combined with
these other equations. Section IV describes a simple algorithm to solve the
integral equations numerically and how the initial conditions for the IE
solver are obtained from the tight-coupling approximation at early times. This
Section also describes an iterative algorithm to solve them in tandem with the
differential equations. Section VI describes the two rudimentary codes to
evolve the Boltzmann hierarchy and the IE equations. I then present and
discuss results of the calculation. Section VII concludes with a discussion of
possible concerns and ideas for further steps. Appendix A provides the photon
Boltzmann equations in the notation used here, and Appendix B provides details
of the algorithm to solve the integral equation. The codes are provided at
https://github.com/marckamion/IE for readers interested to follow up on
calculational details that cannot be inferred from the presentation here.

##  II Formalism

If we have a spectrum of initial curvature fluctuations with power spectrum
$P_{\cal R}(k)=\left\langle|{\cal R}_{\vec{k}}|^{2}\right\rangle$, then the
CMB temperature/polarization power spectra are

$C_{l}^{\rm XX^{\prime}}=(2\pi^{2})^{-1}\int\,k^{2}\,dk\,P_{\cal
R}(k)\Delta_{kl}^{\rm X}(\tau_{0})\Delta_{kl}^{\rm X^{\prime}}(\tau_{0}),$ (1)

for X,X’$=$T,E with “T” the temperature and “E” the E-mode of the
polarization. The transfer functions $\Delta_{kl}^{\rm X}(\tau)$ are obtained
through solution of differential equations for the time evolution of the
relativistic gravitational potentials, the baryon, dark-matter, photon, and
neutrino densities and bulk velocities, and the higher moments of the photon
and neutrino momentum distributions. The moments of the intensity distribution
of photon momenta are the transfer functions $\Delta^{\rm T}_{kl}(\tau)$ and
the moments of the distribution of photon polarizations are $\Delta^{\rm
E}_{kl}(\tau)$.

The temperature transfer functions can be written as111The notation here
resembles largely that in Ref. Blas:2011rf . The differences are that (i) the
photon $\Delta^{\rm T}_{kl}$ here is one quarter of theirs; (ii) the $R$ here
is the inverse of their $R$; (iii) the $\dot{\kappa}$ here is their
$\tau_{C}^{-1}$; (iv) the $\alpha$ here is their $h+6\eta$. The $\Pi$ here is
the same as that in Ref. Seljak:1996is  and is $\Pi=(F_{\gamma 2}+G_{\gamma
0}+G_{\gamma 2})/4$ in terms of the variables in Ref. Blas:2011rf .

$\displaystyle\Delta^{\rm T}_{kl}(\tau)$ $\displaystyle=$
$\displaystyle\int_{\tau_{i}}^{\tau}\,d\tau^{\prime}\,g(\tau,\tau^{\prime})\left\\{\left[-\frac{1}{6}\frac{\dot{h}_{k}(\tau^{\prime})}{\dot{\kappa}(\tau^{\prime})}+\Delta^{\rm
T}_{k0}(\tau^{\prime})\right]j_{l}(x)\right.$ $\displaystyle-$
$\displaystyle\left.\left[\frac{1}{3}\frac{\dot{\alpha}_{k}(\tau^{\prime})}{\dot{\kappa}(\tau^{\prime})}+\frac{1}{2}\Pi_{k}(\tau^{\prime})\right]R^{\rm
LL}_{l}(x)+\theta_{bk}(\tau^{\prime})\frac{j_{l}^{\prime}(x)}{k}\right\\},$

where $x=k(\tau-\tau^{\prime})$; a dot denotes a partial derivative with
respect to $\tau$; and
$g(\tau,\tau^{\prime})=(d/d\tau^{\prime})e^{-\kappa(\tau,\tau^{\prime})}=\dot{\kappa}(\tau^{\prime})\,e^{-\kappa(\tau,\tau^{\prime})}$
is the visibility function. The initial conformal time $\tau_{i}$ must be
taken to be deep in the tight-coupling regime and will be discussed more
below. Here, $\dot{\kappa}(\tau)=d\kappa/d\tau$ is the opacity, the derivative
of the Thomson-scattering optical depth with respect to conformal time, and

$\kappa(\tau,\tau^{\prime})=\int_{\tau^{\prime}}^{\tau}\,d\tau_{1}\,\dot{\kappa}(\tau_{1}).$
(3)

Also, $R^{\rm
LL}_{l}(x)=-\frac{1}{2}\left[j_{l}(x)+3j_{l}^{\prime\prime}(x)\right]$
Hu:1997hp ; Dai:2012bc  in terms of spherical Bessel functions $j_{l}(x)$, and
$\theta_{bk}(\tau)$ is the baryon velocity. It is related to the photon
velocity (suppressing hereafter the subscript $k$ for notational economy)
$\theta_{\gamma}(\tau)=3k\Delta^{\rm T}_{k1}(\tau)$ through

$\dot{\theta}_{b}=-{\cal
H}\theta_{b}+c_{s}^{2}k^{2}\delta_{b}+\frac{\dot{\kappa}}{R}(\theta_{\gamma}-\theta_{b}),$
(4)

where ${\cal H}(\tau)\equiv\dot{a}/a$ and
$R(\tau)\equiv(3/4)\rho_{b}(\tau)/\rho_{\gamma}(\tau)$, the scale factor in
units of $3/4$ of that at matter-baryon equality ($\rho_{b}(\tau)$ and
$\rho_{\gamma}(\tau)$ are mean baryon and photon energy densities,
respectively). The baryon sound speed $c_{s}$ is increasingly important on
small scales but has little effect on the larger distance/angular scales
relevant for CMB fluctuations. Here, $h(\tau)$ is the standard synchronous-
gauge perturbation variable, and $\alpha(\tau)=h(\tau)+6\eta(\tau)$ in terms
of the commonly used synchronous-gauge variable $\eta(\tau)$.

The function $\Pi(\tau)$ is a linear combination of the photon-intensity and
polarization quadrupoles; for simplicity, I refer to it here as the
polarization quadrupole. It can also be written as an IE,

$\Pi(\tau)=\Delta^{\rm T}_{2}(\tau)+9E_{2}(\tau),$ (5)

with

$E_{l}(\tau)=\int_{\tau_{i}}^{\tau}\,d\tau^{\prime}\,g(\tau,\tau^{\prime})\frac{j_{l}(k(\tau-\tau^{\prime}))}{\left(k(\tau-\tau^{\prime})\right)^{2}}\Pi(\tau^{\prime}).$
(6)

The CMB E-mode transfer function is then $\Delta^{\rm
E}_{l}(\tau)=(3/4)\sqrt{(l+2)!/(l-2)!}E_{l}(\tau)$.

Figure 1: Flow charts for the perturbation calculation with (a) the Boltzmann
hierarchy and (b) the integral-equation approach. An arrow points from an
element that appears in the differential equation for the element it points
to. Ingredients that appear in the integral equation for a given quantity are
indicated in (b) with an integral sign. As both figures indicate, the higher
moments ($l\geq 3$ for $\Delta_{l}^{\rm T}$ and $\Delta_{l}^{\nu}$ and $\geq
2$ for $E_{l}$) communicate to the rest of the system of equations only
through the quadrupole ($l=2$). The diagrams also indicate that in both cases,
the photon-intensity quadrupole $\Delta_{2}^{\rm T}$ feeds into the rest of
the system of equations only through the photon velocity $\theta_{\gamma}$,
and similarly for the neutrino quadrupole.

A derivation of Eqs. (LABEL:eqn:Tfluctuation) and (5) will be provided in Ref.
inprep  using the total-angular-momentum formalism Dai:2012bc , but it is
easily verified that they agree with Eq. (18) in Ref. Zaldarriaga:1996xe ,
Eqs. (74) and (77) in Ref. Hu:1997hp , and with the IEs in Ref.
Weinberg:2003ur . It can also be verified, using the relation,
$(2l+1)j_{l}^{\prime}(x)=lj_{l-1}(x)-(l+1)j_{l+1}(x)$ (which $R^{\rm
LL}_{l}(x)$ and $j_{l}^{\prime}(x)$ also satisfy), that differentiation of
these two IEs recovers the usual Boltzmann hierarchy as given, for example, in
Eqs. (2.4) of Ref. Blas:2011rf  or Eq. (63) of Ref. Ma:1995ey . Thus, these
two IEs are formally equivalent to the Boltzmann hierarchy. For completeness,
the Boltzmann hierarchy is provided in the notation/conventions used here in
Appendix A.

##  III Implementation

The left flowchart in Fig. 1 shows the interdependency between the different
perturbation variables in the differential equations for their evolution. In
the middle are the metric-perturbation variables $h$ and $\alpha$. These are
sourced by the baryon, dark-matter, neutrino, and photon densities and bulk
velocities. Apart from the baryon-photon coupling that connects
$\theta_{\gamma}$ and $\theta_{b}$, the only communication between the
different matter components is through the metric perturbations. The neutrino
velocity is connected to the neutrino quadrupole $\Delta^{\nu}_{2}$ which is
then connected to an infinite tower of Boltzmann equations for the higher-
order neutrino moments $\Delta^{\nu}_{l}$ for $l\geq 3$. The same can be said
for the photon velocity, except that there are two infinite Boltzmann
hierarchies for the higher photon-intensity and photon-polarization moments.
When considered in tandem, the photon monopole and dipole equations combine
into a second-order differential equation that resembles that for a driven
simple harmonic oscillater (discussed below); this describes oscillations of
the amplitude of the photon-baryon fluid driven by changes in the metric
perturbations and in the photon quadrupole.

In the line-of-sight approach Seljak:1996is , the Boltzmann hierarchy is
solved up to a maximum multipole $l_{\rm max}\sim 30$ to obtain the photon
monopole, dipole, and quadrupole, and $\Pi$ to reasonable accuracy. The
$C_{l}$ are then obtained by evaluating the integrals in Eqs.
(LABEL:eqn:Tfluctuation) and (5).

As Fig. 1 illustrates, the two (nominally) infinite towers of photon
differential equations—one for the temperature moments ($\Delta_{l}^{\rm T}$
for $l\geq 3$) and polarization moments ($E_{l}$ for $l\geq 2$)—communicate
with the rest of the system of equations only through the photon-intensity
quadrupole $\Delta_{2}^{\rm T}$. Thus, one can replace the two photon
Boltzmann hierarchies with a pair of integral equations, one for
$\Delta_{2}^{\rm T}$ and another for $\Pi$. The rest of the system of
equations is then exactly the same as in the Boltzmann approach.

In this approach we retain the two lowest-order equations, for the photon
monopole ($l=0$) and dipole ($l=1$). These equations are,

$\dot{\Delta}^{\rm
T}_{0}=-\frac{1}{3}\theta_{\gamma}-\frac{1}{6}\dot{h},\quad\dot{\theta}_{\gamma}=k^{2}(\Delta^{\rm
T}_{0}-2\Delta^{\rm T}_{2})-\dot{\kappa}\Theta_{\gamma b},$ (7)

with $\Theta_{\gamma b}(\tau)\equiv\theta_{\gamma}(\tau)-\theta_{b}(\tau)$.
These equations are supplemented by those,

$\dot{\delta}_{b}=-\theta_{b}-\frac{1}{2}\dot{h},\quad\dot{\theta}_{b}=-{\cal
H}\theta_{b}+c_{s}^{2}k^{2}\delta_{b}+\frac{\dot{\kappa}}{R}\Theta_{\gamma
b},$ (8)

for the baryon density and velocity, respectively. There is also an equation,
$\dot{\delta}_{c}=-\frac{1}{2}\dot{h}$, for the CDM-density perturbation (the
CDM peculiar velocity vanishes in synchronous gauge).

The photon quadrupole $\Delta^{\rm T}_{2}(\tau)$ in Eq. (7) is obtained at
early times by the TCA (up to second order in $\dot{\kappa}^{-1}$, as
described in Refs. CyrRacine:2010bk ; Blas:2011rf  for improved
speed/precision). The two equations for the early-time evolution of
$\theta_{\gamma}$ and $\theta_{b}$ can also be replaced by their TCA, with the
slip $\dot{\Theta}_{\gamma b}$ evaluated (again, up to second order
$\dot{\kappa}^{-1}$) CyrRacine:2010bk ; Blas:2011rf .

At later times, the quadrupole is obtained from Eq. (LABEL:eqn:Tfluctuation)
with $l=2$, along with Eq. (5) for the time evolution of $\Pi(\tau)$. With
this approach, the equations in Eq. (7) combine to describe a driven
oscillator damped by the photon quadrupole semianalytic . The photon
quadrupole is provided at early times by the TCA and at later times from the
integral equation.

For completeness, the Einstein equations are

$\ddot{h}+\frac{\dot{a}}{a}\dot{h}=-8\pi Ga^{2}\left[\delta\rho_{\rm
tot}+3\delta p_{\rm tot}\right],$ (9) $\frac{1}{3}(\dot{h}-\dot{\alpha})=8\pi
Ga^{2}\left[\frac{4}{3}\bar{\rho}_{\gamma}\theta_{\gamma}+\frac{4}{3}\bar{\rho}_{\nu}\theta_{\nu}+\bar{\rho}_{\rm
b}\theta_{\rm b}\right],$ (10)

Note that the the Einstein equations are written here in terms of the energy
and momentum densities, but not the anisotropic stress. In this way, the
photon-intensity quadrupole $\Delta_{2}^{\rm T}(\tau)$ communicates with the
rest of the set of perturbation equations only through Eq. (7). The IEs for
massless neutrinos are obtained from those for photons, but setting
$\Pi=\dot{\kappa}=0$. These IEs have come into play in the development of an
effective ultra-relativistic-fluid approximation Blas:2011rf .

##  IV Numerical solution of the integral equations

The IEs here are Volterra equations of the second kind, which are typically
solved as follows Volterrabook ; Press:1992zz . A pair of such equations has
the form,

$f^{\alpha}(t)=\int_{a}^{t}K^{\alpha\beta}(t,s)f^{\beta}(s)ds+g^{\alpha}(t).$
(11)

with $\alpha,\beta=1,2$ (and implied sum over repeated $\alpha,\beta$ not not
$ij$). They are solved on a mesh of $N$ uniformly spaced time steps
$t_{i}=a+ih$ with $i=1,2,\ldots,N$, with $h=(t-a)/N$. The integrals are then
evaluated with the trapezoidal rule. The solution to the IEs are then
$f_{\alpha,0}=g_{\alpha,0}$ and

$\left(\delta^{\alpha\beta}-\frac{1}{2}hK^{\alpha\beta}_{ii}\right)f^{\beta}_{i}=h\left(\frac{1}{2}K^{\alpha\beta}_{i0}f^{\beta}_{0}+\sum_{j=1}^{i-1}K^{\alpha\beta}_{ij}f^{\beta}_{j}\right)+g^{\beta}_{i}.$
(12)

For the pair of Volterra equations we deal with here, the $2\times 2$ matrix
on the left-hand side must be inverted at each time step Press:1992zz . The
ordinary differential equations, which must be solved simultaneously, are
simply stepped forward in time (i.e., Euler integration).

This algorithm works well if the kernels $K^{\alpha\beta}(t,s)$ are smooth and
slowly varying. The visibility function in our integrands are smoothly varying
after decoupling begins to occur, at redshifts $z\lesssim 1400$ ($\tau\gtrsim
230$ Mpc). The perturbation variables that multiply it, as well as the radial
eigenfunctions, are also relatively smooth. The trapezoidal-rule integration
therefore works reasonably well. However, for early conformal times
($\tau\lesssim 230$ Mpc) during tight coupling, when $\dot{\kappa}\gg{\cal
H}$, the visibility function is very sharply peaked at $\tau^{\prime}\to\tau$.
The trapezoidal rule will therefore be inaccurate (unless we take a huge
number of time steps).

To remedy this, and to improve the transition from tight coupling, we replace
the trapezoidal rule in $\Delta\tau^{\prime}$ with one in
$de^{-\kappa(\tau,\tau^{\prime})}$. More precisely, we write the integrand in
terms of the visibility function,
$(d/d\tau^{\prime})e^{-\kappa(\tau,\tau^{\prime})}$, times the more slowly-
varying perturbation variables. The integrals can then be written,

$\displaystyle I(\tau)$ $\displaystyle=$
$\displaystyle\int^{\tau}\,d\tau^{\prime}f(\tau^{\prime})\frac{d}{d\tau^{\prime}}\left[e^{-\kappa(\tau,\tau^{\prime})}\right]f(\tau)$
(13) $\displaystyle\simeq$
$\displaystyle\sum_{n=1}\int_{\kappa_{n}}^{\kappa_{n-1}}d\left(e^{-\kappa(\tau,\tau^{\prime})}\right)\left[f_{n-1}\right.$
$\displaystyle\left.\ \ \
+\left(\frac{df}{d\kappa^{\prime}}\right)_{n-1}(\kappa-\kappa^{\prime})\right],$

where $\kappa_{n}=\kappa(\tau-nh)$, and $h$ is the small conformal-time step.
The remaining $\kappa^{\prime}$ integrals can then be done analytically and
the derivative $df/d\kappa^{\prime}$ approximated by differencing. Details are
provided in Appendix B.

By expanding the integrand $f(\tau)$ to linear order, as in Eq. (13), we
obtain a result that is exact for variations of $f(\tau)$ that are up to
linear in $\kappa$. At early times, this then reproduces the first-order TCA
(to order $\dot{\kappa}^{-1}$), even for one step that is not necessarily
small compared with $\dot{\kappa}^{-1}$. The second-order TCA is then
recovered by evaluating the IE with two time steps. This allows a smooth
transition from the TCA approximation to the IE algorithm in Appendix B as
long as the TCA values for the perturbation variables are stored for at least
two time steps. At late times, the visibility function in Eq. (13) can be
Taylor expanded to linear order in $\Delta\kappa$. Doing so then recovers the
trapezoidal scheme in Eq. (12).

The formula in Eq. (12) requires for each time step $i$ a sum over all earlier
timesteps $j<i$. However, given the $e^{-\kappa(\tau,\tau^{\prime})}$ factor
in the visibility function in the integrand, the sum can for all practical
purposes be started, for any given $\tau_{i}$ at some $j$ such that
$\kappa(\tau_{i},\tau_{j})\leq\Delta\tau_{\rm max}\simeq 10-20$. If the other
factors in the integrand are slowly varying, this yields a precision
degradation of $\lesssim e^{-\Delta\tau_{\rm max}}$.

When the IE solver first begins, the photon-baryon fluid is still tightly
coupled, and so the visibility function has support only over values of
$\tau^{\prime}$ fairly close to $\tau$; i.e., $(\tau-\tau^{\prime})\lesssim
N\dot{\kappa}^{-1}$. The argument $x=k(\tau-\tau^{\prime})$ of the radial
eigenfunctions in Eq. (LABEL:eqn:Tfluctuation) is thus small, and so the
radial eigenfunctions can be approximated as $j_{2}(x)\simeq x^{2}/15$,
$R^{\rm LL}_{2}(x)\simeq-1/5$, $j_{2}^{\prime}(x)\simeq(2/15)x$. The integrand
cannot, however, be approximated simply by the $R^{\rm LL}(x)$ term, because
$\Pi$ is ${\cal O}(\dot{\kappa}^{-1})$ times $\theta_{b}$. The third (i.e.,
the $\theta_{b}$) term contributes, at lowest order in the TCA.

##  V Iterative solution of integral and differential equations

The next step is to consider how to solve simultaneously the differential
equations for the rest of the system. This includes those for the metric-
perturbation variables, and the baryon and dark-matter densities and
velocities. It also in principle includes the neutrino perturbation variables;
here, however, I will assume that these can be obtained with a generalized-
dark-matter Hu:1998kj  or ultrarelativistic-fluid approximation (UFA)
Blas:2011rf , both of which have been made fairly effective. In principle, the
integral-equation techniques described for photons here can be applied to the
neutrino sector as well. For clarity, I focus here, though, on the photon
sector.

In trying to do so, however, the coupling between the IEs and the DEs pose a
chicken-and-egg problem: The differential equations for the rest of the system
require knowledge of $\Delta_{2}^{\rm T}(\tau)$, but the IEs for
$\Delta_{2}^{\rm T}(\tau)$ cannot be obtained without the solution to the DEs.
One possibility is to solve the IEs and DEs simultaneously by simply stepping
the differential equations forward—i.e., Euler integration. This, however,
requires very fine time steps, especially toward the end of the TCA, and thus
eliminates the advantages of the early-time IE algorithm described above.
Another possibility is to step the IEs forward on a coarse time grid, and then
integrate the DEs forward (using an extrapolation of the IE solutions from
earlier time steps) using an off-the-shelf adaptive-time-step DE solver.

However, the IEs and DEs can be solved very efficiently with a simple
iterative algorithm. Here, we start with some initial anzatz for
$\Delta_{2}^{\rm T}(\tau)$ and $\Pi(\tau)$ and then solve the DEs for all the
other perturbation variables with this ansatz. We then integrate the IEs using
the solutions to those DEs to obtain new values of $\Delta_{2}^{\rm T}(\tau)$
and $\Pi(\tau)$. We then iterate. Of course, there is no guarantee a priori
that this iterative procedure will converge to the correct answer, but some
simple numerical experiments show that this procedure converges, and does so
fairly quickly, even for a lousy (e.g., $\Delta_{2}^{\rm T}(\tau)=\Pi=0$)
initial ansatz for the IE solutions.

Figure 2: The CMB visibility function $\dot{\kappa}(\tau_{0},\tau)$ as a
function of conformal time. It is shown to indicate the range of conformal
times, peaked at $\tau\simeq 280$ Mpc, that contribute to the observed CMB
power spectra from recombination. Figure 3: The transfer function
$\Delta_{2}^{\rm T}(\tau)$ for the CMB photon-intensity quadrupole as a
function of conformal time $\tau$ for a Fourier mode of wavenumber $k=0.2$ Mpc
(which corresponds roughly to a CMB multipole moment $l\sim 3000$). The black
curve shows the results of the full Boltzmann hierarchy as a function of
conformal time. The other curves show results of the iterative integral-
equation solution, taking $\Delta_{2}^{\rm T}(\tau)=0=\Pi(\tau)$ as an initial
ansatz. The yellowish curve shows the result for $\Delta_{2}^{\rm T}(\tau)$
after the first iteration—i.e., after integrating the differential equations
for all perturbation variables except $\Delta_{2}^{\rm T}(\tau)$ and
$\Pi(\tau)$ and then integrating the integral equations for $\Delta_{2}^{\rm
T}(\tau)$ and $\Pi(\tau)$ using the results of the differential equations. The
red curve shows results after three iterations, and the blue curve after five
iterations. The thickness of the curves is such that if two are
indistinguishable, the agreement between the two is $O(0.1\%$).

##  VI Numerical results

I have written a rudimentary C code to calculate the transfer functions for
the perturbation variables with the iterative numerical implementation
described here. To simplify, I approximate neutrinos (taken to be massless) as
a generalized-dark-matter component with $w=c_{s}^{2}=c_{\rm vis}^{2}=1/3$
Hu:1998kj . I stop the code at redshift $z\simeq 560$, after recombination but
before reionization, and use an analytic approximation (which takes into
account only radiation and nonrelativistic matter at these times) for the
expansion history. I use an ionization history from HyRec-2 hyrec . To compare
this IE approach with the standard Boltzmann hierarchy, I also wrote a second
code that is identical in every way except that it swaps out the integral
equations for $\Delta_{2}^{\rm T}(\tau)$ and $\Pi(\tau)$ for the complete
photon Boltzmann hierarchy. The code uses an off-the-shelf differential-
equation solver Shampine  with adaptive step size, although not necessarily
optimized for stiff equations.

In the IE code, the handoff from the TCA to the IE solver takes place at
$\tau=160$ Mpc. The Boltzmann code uses the same TCA at early times and then
starts the full Boltzmann hierarchy at $\tau=160$ Mpc. The Boltzmann code
follows the Boltzmann hierarchy up to $l_{\rm max}=50$ (which I found was
required to keep the perturbation variables stable over the $\tau$ range
considered here). The results are similar, and the code a bit quicker, for
smaller $l_{\rm max}$. The differential-equation solver in the Boltzmann code
runs with a relative error requirement of $10^{-5}$ and absolute error of
$10^{-4}$. The integral equations are evolved on a time grid that has spacing
$\Delta\tau=1.0$ from $160\,{\rm Mpc}\leq\tau\leq 240\,{\rm Mpc}$ and
$350\,{\rm Mpc}\leq\tau\leq 450\,{\rm Mpc}$, and $\Delta\tau=0.5$ Mpc for
$240\,{\rm Mpc}\leq\tau\leq 350\,{\rm Mpc}$, for a total of 401 grid points.
The time required for the IE part of the calculation scales as the square of
the number of grid points.

Fig. 2 shows the visibility function, which indicates the conformal-time
regime, $250\,{\rm Mpc}\lesssim\tau\lesssim 400\,{\rm Mpc}$, over which the
source functions for the CMB power spectra are evaluated.

Fig. 3 illustrates the results of the numerical experiment. Shown there are
results for the photon-intensity quadrupole $\Delta_{2}^{\rm T}(\tau)$ of the
Boltzmann code and the iterative integral-equation results, starting from a
naive initial ansatz $\Delta_{2}^{\rm T}(\tau)=\Pi(\tau)=0$. Results are shown
for $k=0.2$ Mpc, which corresponds roughly to CMB multipole moments $l\sim
3000$, near the upper limit of current measurements. The frequency of
oscillations in the transfer function are reduced at smaller $k$, and so the
numerical algorithm should, if anything, work even better at lower $k$.

The results are shown for one iteration (yellow), three iterations (red) and
(five iterations) blue. The iterative solutions converge first at early times
and then require more iterations to converge at later times. The overlap
between the black and blue (5 iterations) curves indicates that the agreement
is at the $O(0.1\%)$ level over the conformal-time range that contributes to
the observed CMB power spectra. This IE code takes $\sim 0.15$ times as long
to run as the Boltzmann code, implying that each iteration can be completed in
$\sim 1/30$ the time required for the Boltzmann code. Both codes are fairly
rudimentary, and so these time comparisons should be taken with a grain of
salt. Still, these results suggest that this may provide a route to speeding
up the standard Boltzmann codes.

There may be room for even further improvement. The results shown in Fig. 3
are obtained using the most naive possible initial ansatz for $\Delta_{2}^{\rm
T}(\tau)$ and $\Pi(\tau)$. The number of iterations required for convergence
to the required precision can be reduced if one starts with a better initial
guess for these quantities. It should be possible to derive a simple semi-
analytic ansatz that interpolates between the well-understood early-time TCA
behavior and the late-time behavior, which comes from the Sachs-Wolfe effect.

One should, however, be able to do even better. These calculations are not
performed in isolation. In cosmological MCMC analyses, the Boltzmann codes are
run repeatedly to map the likelihood functions in a multidimensional
cosmological-parameter space. Thus, each time the calculation is done, it has
presumably already been done for a nearby point in that cosmological parameter
space. Thus, it should be possible to start the iterative algorithm by using
the results for $\Delta_{2}^{\rm T}(\tau)$ and $\Pi(\tau)$ from a previous
run. To test this, I ran the code using as the initial ansatz the results for
$\Delta_{2}^{\rm T}(\tau)$ and $\Pi(\tau)$ from a prior run with $\Omega_{b}$
reduced by 2%. This code converges to $O(0.1\%)$ after just one iteration.

##  VII Conclusions and ideas for future work

I have presented an alternative formulation of the equations for the evolution
of cosmological perturbations in which the infinite Boltzmann hierarchy for
the photon distribution function is replaced by a pair of integral equations.
There is no new physics here—it is simply a recasting of the equations in a
way that may lead to physical insight and alternative schemes for numerical
solution. As was known from the line-of-sight approach Seljak:1996is , CMB
fluctuations are determined only by the photon monopole (energy density),
dipole (peculiar velocity), and quadrupole (more specifically, $\Pi$). In the
Boltzmann hierarchy, these are the result of some complicated transfer of
power between these lower moments of the photon distribution function and an
infinite tower of higher moments. The IE formalism shows, however, that the
lower moments, and in particular the quadrupole moment, at the surface of last
scatter (i.e., those that enter into the line-of-sight integration) are simply
described by the exact same equations that describe the lower moments that we
see.

I have shown that simple iterative solution of the combined system of integral
and differential equations does a pretty good job at reproducing the results
of the Boltzmann calculation in a fraction of the time. This exercise also
shows that the IE formalism can be implemented numerically without
(apparently) any significant numerical instabilities—this was not a foregone
conclusion, given the occurrence of instabilities in some IE solvers
Volterrabook , as well as those that may arise from finite $l_{\rm max}$ in
the Boltzmann hierarchy.

There is, however, far more work that needs to be done before we know whether
this approach can implemented to speed up a code like CLASS or CAMB. These
codes benefit from a number of insights and clever algorithms, whereas what I
have presented here is fairly naive. Those codes also have controlled errors,
whereas the grid spacing in my calculation was guessed to provide an
$O(0.1\%)$ precision in $\Delta_{2}^{\rm T}(\tau)$.

The spacing of the conformal-time grid points in the integral-equation solver
is an obvious thing to explore. In this calculation I simply estimated the
number of grid points that would be required for $O(0.1\%)$ precision.
However, the distribution of grid points can certainly be optimized to provide
the desired observables (e.g., CMB and matter power spectra) to the required
precision. Good results can probably also be obtained for smaller $k$ with
fewer grid points, given the smoother integrands at lower $k$. The current
code also sums over all prior grid points. However, given the high opacity at
early times, the sum can be restricted only to grid points that are at an
optical depth $\Delta\kappa\lesssim 5$ earlier. There are also algorithms,
more sophisticated than the trapezoidal-rule algorithm used here, on numerical
solution to Volterra equations (e.g., Ref. Baker:2000 ) in the literature that
may be worth exploring. Finally, there may be alternative implementations of
the integral/differential equations that may be better suited for numerics.
For example, it should be possible to eliminate the differential equations for
the photon monopole and dipole and replace the integral equation for the
quadrupole $\Delta_{2}^{\rm T}(\tau)$ with that for the monopole
$\Delta_{0}^{\rm T}(\tau)$. Or perhaps the differential equation for
$\Delta_{2}^{\rm T}(\tau)$ can be included and the integral equation replaced
by one for $\Delta_{3}^{\rm T}(\tau)$.

###### Acknowledgements.

I thank L. Ji, R. Caldwell, D. Grin, J. Bernal, and E. Kovetz for useful
discussions and comments on an earlier draft. This work was supported by NSF
Grant No. 1818899 and the Simons Foundation.

##  Appendix A Boltzmann hierarchy

For completeness and comparison with prior work, I provide the Boltzmann
equations for the photon moments in the notation used here. These equations
are derived by differentiating Eqs. (LABEL:eqn:Tfluctuation) and (5) with
respect to $\tau$. The independent variable $\tau$ appears in the limit of
integration, the visibility function, and in the radial eigenfunctions, and
all of the radial eigenfunctions satisfy the spherical-Bessel-function
relation, $(2l+1)j_{l}^{\prime}(x)=lj_{l-1}(x)-(l+1)j_{l+1}(x)$. The monopole
and dipole equations are already provided in Eq. (7). The equations for $l\geq
2$ are

$\displaystyle\dot{\Delta}^{\rm T}_{l}$ $\displaystyle=$
$\displaystyle-\dot{\kappa}\Delta^{\rm T}_{l}+\frac{kl}{2l+1}\Delta^{\rm
T}_{l-1}-\frac{k(l+1)}{2l+1}\Delta^{\rm T}_{l+1}$ $\displaystyle\ \ \ \ \ \
+\frac{1}{5}\left(\frac{\dot{\alpha}}{3}+\frac{\dot{\kappa}\Pi}{2}\right)\delta_{l2},$
$\displaystyle\dot{E}_{l}$ $\displaystyle=$
$\displaystyle-\dot{\kappa}E_{l}+\frac{k(l-2)}{2l+1}E_{l-1}-\frac{k(l+3)}{2l+1}E_{l+1}+\frac{1}{15}\dot{\kappa}\Pi\delta_{l2},$

with $\Pi=\Delta^{\rm T}_{l}+9E_{l}$.

##  Appendix B Details of the IE solver

We first define functions $I^{\rm T}(\tau,\tau^{\prime})$ and
$I^{\Pi}(\tau,\tau^{\prime})$ by writing

$\displaystyle\Delta^{\rm T}_{2}(\tau)$ $\displaystyle=$
$\displaystyle\int^{\tau}d\tau^{\prime}g(\tau,\tau^{\prime})I^{\rm
T}(\tau,\tau^{\prime}).$ $\displaystyle\Pi(\tau)$ $\displaystyle=$
$\displaystyle\int^{\tau}d\tau^{\prime}g(\tau,\tau^{\prime})I^{\Pi}(\tau,\tau^{\prime}).$
(15)

The integrals are then discretized, taking into account the fact that
$\Pi(\tau)$ appears in $I^{\rm T}(\tau,\tau^{\prime})$ and
$I^{\Pi}(\tau,\tau^{\prime})$, in the following way. We define two sums,

$\displaystyle\Delta^{0}_{2,i+1}$ $\displaystyle=$ $\displaystyle\sum_{j\leq
i}\left(I^{\rm T}_{j+1}W^{+}_{j}+I^{\rm
T}_{j}W_{j}\right)-\frac{1}{10}\Pi_{j+1}W^{+}_{i}$
$\displaystyle\Pi^{0}_{i+1}$ $\displaystyle=$ $\displaystyle\sum_{j\leq
i}\left(I^{\Pi}_{j+1}W^{+}_{j}+I^{\Pi}_{j}W_{j}\right)-\frac{3}{5}\Pi_{i+1}W^{+}_{j},$

where $\Pi_{i}=\Pi(\tau_{i})$, $I^{\rm T}_{j}=I^{\rm T}(\tau_{i+1},\tau_{j})$,
and $I^{\Pi}_{j}=I^{\Pi}(\tau_{i+1},\tau_{j})$. Here the weight functions are

$\displaystyle W^{+}_{j}$ $\displaystyle=$ $\displaystyle
e^{-\kappa(\tau_{i+1},\tau_{j+1})}\left(1-e^{-\Delta\kappa_{j}}-\frac{1-(1+\Delta\kappa_{j})e^{-\Delta\kappa_{j}}}{\Delta\kappa_{j}}\right),$
$\displaystyle W_{j}$ $\displaystyle=$
$\displaystyle\frac{e^{-\kappa(\tau_{i+1},\tau_{j+1})}}{\Delta\kappa_{j}}\left[1-(1+\Delta\kappa_{j})e^{-\Delta\kappa_{j}}\right]$
(16)

where $\Delta\kappa_{j}=\kappa(\tau_{j+1})-\kappa(\tau_{j})$. These weight
functions approach $W^{+}_{j}\to\Delta\kappa_{j}/2$ and
$W_{j}\to\Delta\kappa_{j}/2$ at late times, thus recovering Eq. (12) (written
as an integral over $\kappa$, rather than $\tau$). At early times,
$W_{j}^{+}\to 1-(\Delta\kappa)^{-1}$ and $W^{j}\to(\Delta\kappa)^{-1}$; this
then recovers the first-order tight-coupling approximation,
$\Delta_{2}=(2/5)\Pi=(4/45)(\dot{\alpha}+2\theta_{b})/\dot{\kappa}$, even from
one time step in the evaluation of the integral—the second-order TCA is
reproduced by two time steps.

The discretized quadrupoles are then,

$\displaystyle\Pi_{i+1}$ $\displaystyle=$
$\displaystyle\frac{\Pi^{0}_{i+1}+\Delta^{0}_{2,i+1}}{1-\frac{7}{10}W^{+}_{i}},$
$\displaystyle\Delta^{\rm T}_{2,i+1}$ $\displaystyle=$
$\displaystyle\Delta^{0}_{2,i+1}+\frac{1}{10}\Pi^{0}_{i+1}W^{+}_{i}.$ (17)

## References

  * (1) R. A. Sunyaev and Y. .B. Zeldovich, “Small scale fluctuations of relic radiation,” Astrophys. Space Sci. 7, 3 (1970); P. J. E. Peebles and J. T. Yu, “Primeval adiabatic perturbation in an expanding universe,” Astrophys. J. 162, 815 (1970); J. Silk, “Fluctuations in the Primordial Fireball,” Nature 215, no.5106, 1155-1156 (1967)
  * (2) J. R. Bond and G. Efstathiou, “Cosmic background radiation anisotropies in universes dominated by nonbaryonic dark matter,” Astrophys. J. 285, L45 (1984); J. R. Bond and G. Efstathiou, “The statistics of cosmic background radiation fluctuations,” Mon. Not. Roy. Astron. Soc. 226, 655 (1987); M. L. Wilson and J. Silk, “On the Anisotropy of the cosomological background matter and radiation distribution. 1. The Radiation anisotropy in a spatially flat universe,” Astrophys. J. 243, 14 (1981); N. Vittorio and J. Silk, “Fine-scale anisotropy of the cosmic microwave background in a universe dominated by cold dark matter,” Astrophys. J. 285, L39 (1984);
  * (3) U. Seljak and M. Zaldarriaga, “A Line of sight integration approach to cosmic microwave background anisotropies,” Astrophys. J. 469, 437-444 (1996) [arXiv:astro-ph/9603033 [astro-ph]].
  * (4) F. Y. Cyr-Racine and K. Sigurdson, “Photons and Baryons before Atoms: Improving the Tight-Coupling Approximation,” Phys. Rev. D 83, 103521 (2011) [arXiv:1012.0569 [astro-ph.CO]].
  * (5) D. Blas, J. Lesgourgues and T. Tram, “The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes,” JCAP 07, 034 (2011) doi:10.1088/1475-7516/2011/07/034 [arXiv:1104.2933 [astro-ph.CO]].
  * (6) A. Lewis, A. Challinor and A. Lasenby, “Efficient computation of CMB anisotropies in closed FRW models,” Astrophys. J. 538, 473-476 (2000) [arXiv:astro-ph/9911177 [astro-ph]].
  * (7) S. Weinberg, “Damping of tensor modes in cosmology,” Phys. Rev. D 69, 023503 (2004) [arXiv:astro-ph/0306304 [astro-ph]].
  * (8) D. Baskaran, L. P. Grishchuk and A. G. Polnarev, “Imprints of Relic Gravitational Waves in Cosmic Microwave Background Radiation,” Phys. Rev. D 74, 083008 (2006) [arXiv:gr-qc/0605100 [gr-qc]].
  * (9) R. Flauger and S. Weinberg, “Tensor Microwave Background Fluctuations for Large Multipole Order,” Phys. Rev. D 75, 123505 (2007) [arXiv:astro-ph/0703179 [astro-ph]].
  * (10) J. R. Pritchard and M. Kamionkowski, “Cosmic microwave background fluctuations from gravitational waves: An Analytic approach,” Annals Phys. 318, 2-36 (2005) [arXiv:astro-ph/0412581 [astro-ph]].
  * (11) S. Weinberg, “A No-Truncation Approach to Cosmic Microwave Background Anisotropies,” Phys. Rev. D 74, 063517 (2006) [arXiv:astro-ph/0607076 [astro-ph]].
  * (12) W. Hu and M. J. White, “CMB anisotropies: Total angular momentum method,” Phys. Rev. D 56, 596-615 (1997) [arXiv:astro-ph/9702170 [astro-ph]].
  * (13) L. Dai, M. Kamionkowski and D. Jeong, “Total Angular Momentum Waves for Scalar, Vector, and Tensor Fields,” Phys. Rev. D 86, 125013 (2012) [arXiv:1209.0761 [astro-ph.CO]].
  * (14) L. Ji and M. Kamionkowski, in preparation.
  * (15) M. Zaldarriaga and U. Seljak, “An all sky analysis of polarization in the microwave background,” Phys. Rev. D 55, 1830-1840 (1997) [arXiv:astro-ph/9609170 [astro-ph]].
  * (16) C. P. Ma and E. Bertschinger, “Cosmological perturbation theory in the synchronous and conformal Newtonian gauges,” Astrophys. J. 455, 7-25 (1995) [arXiv:astro-ph/9506072 [astro-ph]].
  * (17) W. Hu and N. Sugiyama, “Anisotropies in the cosmic microwave background: An Analytic approach,” Astrophys. J. 444, 489 (1995) [astro-ph/9407093];
  * (18) P. Linz, “Analytical and Numerical Methods for Volterra Equations,” (Society for Applied and Numerical Mathematics, 1985).
  * (19) W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, “Numerical Recipes in FORTRAN: The Art of Scientific Computing,” (Cambridge University Press, Cambridge, 1992).
  * (20) W. Hu, “Structure formation with generalized dark matter,” Astrophys. J. 506, 485-494 (1998) [arXiv:astro-ph/9801234 [astro-ph]].
  * (21) Y. Ali-Haimoud and C. M. Hirata, “HyRec: A fast and highly accurate primordial hydrogen and helium recombination code,” Phys. Rev. D 83, 043513 (2011) [arXiv:1011.3758 [astro-ph.CO]]; N. Lee and Y. Ali-Haïmoud, “HYREC-2: a highly accurate sub-millisecond recombination code,” Phys. Rev. D 102, no.8, 083517 (2020) [arXiv:2007.14114 [astro-ph.CO]].
  * (22) L. Shampine and M. Gordon, “Computer Solution of Ordinary Differential Equations: The Initial Value Problem,” (Freeman, 1975).
  * (23) C. T. H. Baker, “A perspective on the numerical treatment of the Volterra equations,” J. Comp. App. Math. 1250, 217 (2000).