{'fraction_non_alphanumeric': 0.03588705141185735, 'fraction_numerical': 0.0075374057824277195, 'mean_word_length': 4.363197586726998, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A No-Truncation Approach to Cosmic Microwave Background Anisotropies

We offer a method of calculating the source term in the line-of-sight integral for cosmic microwave background anisotropies without using a truncated partial-wave expansion in the Boltzmann hierarchy.


I. Introduction
Originally the Boltzmann equation for the photon distribution in cosmology [1] was solved numerically by expanding the components of the photon density matrix in a series of Legendre polynomials [2], but to get results that could be compared with observation of the cosmic microwave background this method requires the inclusion of hundreds or even thousands of partial waves, requiring hours or even days of computer time for each theoretical model. A great improvement was introduced with the suggestion to use instead a formal solution of the Boltzmann equation, in the form of a "line of sight" integral [3]. But this is still only a formal solution, in the sense that we still need to calculate source terms appearing in the integrand. These terms involve partial waves for the photon distribution up to ℓ = 2 for the scalar modes and ℓ = 4 for the tensor modes, and these of course are coupled to higher partial waves. In the original proposal of the "line of sight" method, and in the computer programs CMBfast and CAMB based on this method, these source terms are calculated numerically, by first finding an approximate solution of the Boltzmann equations for partial wave amplitudes. An accurate solution for the partial waves appearing in the source terms can be found by truncating the partial wave expansion at a sufficiently high value of ℓ. In the latest version of CMBfast, the integrand for scalar modes is calculated using partial waves up to ℓ = 12, in which case one has to solve at least 26 coupled ordinary differential equations for the evolution of the partial wave amplitudes, not counting the equations needed to follow the evolution of the baryonic plasma, cold dark matter, neutrinos, and gravitational field components. For tensor modes, the source terms are calculated by solving the 22 differential equations for partial waves with up to ℓ = 10 for photons. The results for ℓ ≤ 2 or ℓ ≤ 4 are used in this method to calculate the integrand of the line-of-sight integral, which then is used to calculate all the higher partial wave amplitudes measured in observations of the cosmic microwave background, up to values of ℓ over 1000.
This article will present an alternative approach, which does not use partial wave expansions to calculate the source terms, and hence obviates the need for any truncation of this expansion. Instead of a large number of coupled differential equations for the partial waves, we have a single integral equation for the tensor modes, and a trio of coupled integral equations for the scalar modes (including one for the plasma velocity). Of course, integral equations are generally harder to solve numerically than differential equations (no routine for solving them is supplied by Mathematica), but in the case at hand the integral equations can be solved numerically by simple iteration. In this method, the calculation itself provides an immediate way of judging its own reliability -if the nth iteration agrees with the n−1th iteration to a satisfactory degree of accuracy, one has a solution. In a sample calculation of the source term for the tensor modes, the results converge rapidly in just a few iterations. This paper concentrates in the next section on the calculation of the photon distribution, but a truncated partial wave expansion is also unnecessary for neutrinos. Indeed, it is already known [4] that the momentum distribution of massless neutrinos for a given metric perturbation can be calculated in terms of a simple line-of-sight integral, with no need to solve integral equations. CMBfast does not use this line-of-sight method for neutrinos, presumably because no one is interested in very high partial waves in the neutrino distribution, but to get good accuracy for the neutrino contribution to the energy-momentum tensor it carries the partial wave expansion to ℓ max = 25. In the last section of this paper the approach of reference [4], which dispenses with partial wave expansions, is extended to the case of massive as well as massless neutrinos, and to scalar as well as tensor modes.

II. Photons
First, some reminders about the Boltzmann equation for photons in cosmology. We will adopt a coordinate system in which the metric takes the form where h ij is a first-order perturbation. Weakly perturbed metrics will automatically be of this form in tensor modes, and can be put in this form for scalar modes by adopting a synchronous gauge.
For our purposes, it is important to write the Boltzmann equation for the photon distribution in a matrix form, rather than in the partial wave formalism in which it is usually presented. The photon distribution is described by a polarization density matrix n ij (x, p, t), defined so that if we measure whether photons have polarization in a direction e i rather than in an orthogonal direction, then the number of photons with polarization e i in a volume i dp i dx i of phase space at time t will be found to be g ik g jl e k e l n ij (x, p, t) m dp m dx m , with p i n ij = 0. (The polarization of a photon with 3-momentum p i is described by a polarization vector e i , satisfying p i e i = 0 and g ij e i e j = 1.) For small perturbations, this matrix can be put in the form Heren γ (p) is the equilibrium phase space number densitȳ (which is a time-independent function of its argument because in the era of interest T (t) ∝ 1/a(t)), and δn ij is a small perturbation. This perturbation satisfies a linearized Boltzmann equation: where p ≡ √ p i p i ,p k ≡ p k /p, δu k is the streaming velocity of the baryonic plasma, and ω c is the frequency with which a photon collides with electrons in the plasma. Instead of δn ij , it is sufficient to consider the intensity matrix perturbation J ij (x,p, t), defined by whereρ γ (t) ≡ a −4 (t) 4πp 3n γ (p) dp is the mean photon energy density. (This is all we need to calculate the photon contributions to the perturbations in the energy-momentum tensor.) To derive the Boltzmann equation for J ij (x,p, t) we multiply Eq. (4) with 4πp 3 and integrate over p ≡ √ p i p i , and find We can also calculate the equation of motion of the plasma, using the perturbed momentumconservation equation: whereρ B is the unperturbed density of the baryonic plasma. These partial differential equations can be converted to ordinary differential equations by writing the perturbed metric as a Fourier integral and looking for solutions in the form Then Eqs. (6) and (7) become ordinary differential equations, with ∂/∂x i replaced with iq i .
They have a formal solution in the form of a "line of sight" integral [3], which in the general case may be written We have here introduced the convenient abbreviations and, as usual, If we take the initial time t 1 early enough so that photons are in local thermal equilibrium with the plasma at that time, then It is the calculation of the source terms J ij (q, t) and I i (q, t) that concerns us in this paper. For this purpose, we now need to distinguish between tensor and scalar modes. We first consider tensor modes, which are computationally simpler.

Tensor Modes
For tensor modes the metric perturbation takes the form where β(q, λ) is a stochastic parameter; e ij (q, λ) is a polarization tensor for helicity λ, witĥ q i e ij (q, λ) = 0 and e kk (q, λ) = 0; and D q (t) is the solution of the wave equation that does not decay while outside the horizon. Here π T q (t) is the coefficient of λ e ij (q, λ)β(q, λ) in the Fourier transform of the tensor part of the anisotropic inertia tensor. (We are considering times sufficiently late so that the decaying solution makes a negligible contribution to δg ij .) The quantity J ij (q, t) will then take the form of a corresponding sum over graviton with J ij (q, t, λ) ordinary c-number functions, not stochastic fields, satisfying the Boltzmann where (The velocity perturbation δu i is absent in tensor modes.) Furthermore, because J ij (q, t, λ) for a given helicity λ must be a linear combination of the polarization tensor components e kl (q, λ) with the same λ, whileq k e kl (q, λ) and e kk (q, λ) both vanish, the only possible form of J ij (q, t, λ) allowed by rotational invariance is just e ij (q, λ) times some function of q ≡ |q| and t. This relation is conventionally written To make contact with the notation used in the usual calculation of the source term Ψ(q, t), we note that the intensity matrix perturbation may be written in the form (Here the superscript T stands for 'tensor,' while the subscript T stands for 'temperature.' The coefficients are chosen so that J ii is proportional to ∆ T , and the polarization is pro- To calculate the source term Ψ(q, t) in terms of partial waves, one first integrates Eq. (21) overp, and finds The functions ∆ and then Eq. (24) reads [5]: As already mentioned, experience shows that to accurately calculate the partial wave amplitudes up to ℓ = 4, which appear in Eq. (27), one needs to solve the Boltzmann equations for the partial wave amplitudes up to larger values of ℓ, up to ℓ = 10. Once Ψ is calculated in this way, we can calculate J ij (q,p, t, λ) for very much higher values of ℓ by using the "line of sight" integral (10), which for the tensor modes gives (In tensor modes there is no perturbation to either the temperature or the velocity of the baryonic plasma, so Eq. (14) gives the initial value J ij (q,p, t 1 ) = 0 for an initial time t 1 taken sufficiently early so that photons are in local thermal equilibrium with the plasma.) Here we suggest the alternative, of deriving an integral equation for Ψ(q, t) by simply analytically integrating Eq. (28) overp. Equating the coefficients of e ij on both sides gives the integral equation Here K(v) and F (v) are the functions Integral equations of this sort are harder to solve numerically than differential equations, because in a step-by-step calculation it is necessary to keep track of Ψ(t ′ ) for all t ′ < t in order to calculate Ψ(t). On the other hand, such equations can also be solved by simple iteration, in which the numerical work is reduced to doing a few integrals. We calculate the nth approximation Ψ (n) to Ψ by using the previous approximation Ψ (n−1) in the second term in square brackets in Eq. (29), and start this calculation by taking the lowest approximation Ψ (0) as just Eq. (29) with the second term in square brackets dropped: Using Ψ (0) for Ψ in Eq. (28) would give precisely the same photon distribution function that we would find if we assumed that photons remained unpolarized until the last scattering; subsequent iterations take account of the polarization produced in multiple scattering. When after n iterations we find that Ψ (n) = Ψ (n−1) to an adequate degree of accuracy, we have a solution of Eq. (29). This is in contrast with the usual truncation method, in which the accuracy of the calculation for a given maximum multipole order can only be judged by comparing with the results for a higher maximum multipole order.
To test the convergence of this iteration procedure, we adopt the usual ΛCDM model, with cosmological parameters We calculate the collision rate ω c by solving the kinetic equations for hydrogen recombination, with an un-ionized helium fraction Y = 0.26. The numerical calculations are done with q chosen so that the physical wave number q/a comes within the horizon at just the time of matter-radiation equality. The gravitational wave amplitude D q (t) is calculated ignoring the effect of anisotropic inertia. Over the interesting time interval when the matter/radiation density ratio y increases from 2 to 4, which includes the time of recombination, the first iteration Ψ (1) q (t) has roughly the same time-dependence as Ψ (0) q (t), but is about 13% to 33 % larger. On the other hand, after five iterations we get a result Ψ (5) (t) that at worst (at y ≃ 3) is within 0.3% of the previous iteration Ψ (4) , and is much closer at other values of y.
As a further illustration of the convergence of the iteration procedure, Figure 1 shows the first few iterations for several wave numbers and a broader range of values of y.
We will not go on here to use this source function in Eq. (28) to calculate the intensity matrix perturbation, because we are not proposing anything new in that part of the calculation of the microwave background anisotropies, but only in the calculation of the source terms appearing in the integrand of the line of sight integral. The sample calculation described here shows that this is a practical method of calculating the source terms, as well as having the conceptual advantage of avoiding a more-or-less arbitrary truncation of a partial wave expansion.

Scalar Modes
For scalar modes, the metric perturbation in synchronous gauge takes the form Assuming the perturbation to be dominated by a single mode (presumably the adiabatic mode that does not decay while outside the horizon), the dependence of A(q, t) and B(q, t) on the direction of q is entirely contained in a stochastic factor β(q) for this mode: In this case, the source terms (12) and (13) appearing in the integrand of the line-of-sight integral must take the form Then, using Eq. (11) to evaluate the plasma velocity, the line-of-sight formula (10) becomes Conventionally, the intensity matrix perturbation for scalar modes is written with ∆ By expanding ∆ and integrating Eq. (37) over the directions ofp, one obtains well known expressions for the source functions in terms of the partial wave amplitudes with ℓ ≤ 2: Instead, we suggest the alternative of deriving coupled integral equations for Φ, Π, and Ω by analytically integrating overp in Eq. (36). This gives where (We have dropped the terms involving the plasma temperature and velocity perturbations at the initial time t 1 , because as long as the collision rate is much larger than the expansion rate the photon distribution is given in terms of δT /T and δu by an equilibrium formula like Eq. (14), and as long as the perturbation is outside the horizon δT /T and δu grow in the adiabatic mode, so that if t 1 is taken sufficiently early in the era when the collision rate is large and the perturbation is outside the horizon, the initial-data terms become negligible.) One may again expect that these coupled integral equations may be solved by iteration, as done for the tensor modes.

III. Neutrinos
In the absence of collisions, the Boltzmann equation for either massless or massive neutrinos for the metric (1) reads Here n(x, p, t) is the phase space density of neutrinos, regarded as a function of the components x i , p i and the time t, while p i and p 0 are functions of the x i and t as well as the p i , given by p i = g ij p j and p 0 = g ij p i p j + m 2 . (The derivation of Eq. (58) is given in Ref. [4] for massless neutrinos, but precisely the same derivation applies also for massive neutrinos.) For weak perturbations, we write wheren ν is a time-independent function of its argument (60) Here δn is a small correction, representing the dynamical rather than purely geometric effect of metric perturbations on the neutrinos, and t 1 is the time the neutrinos went out of equilibrium with the baryonic plasma. (Probably all neutrinos have masses much less than k B T (t 1 ), in which case the term a 2 (t 1 )m 2 in the square root can be neglected.) Expanding to first order in δn and δg ij = a 2 h ij , after many cancellations this gives We again write h ij as a Fourier integral (8), and seek a solution in the form δn(x, p, t) = d 3 q e iq·x δn(q, p, t) .
Then Eq. (61) becomes an ordinary differential equation, with ∂/∂x i replaced with iq i .
Instead of solving this equation numerically with a truncated partial wave expansion, it can be solved analytically, as a line of sight integral where t 1 is any time taken early enough so that δn(x, p, t 1 ) is negligible.
For the foreseeable future, the neutrino distribution will be important only in calculating the neutrino contribution to the energy-momentum tensor. The first-order perturbation to the contribution of each species of neutrino or antineutrino to the mixed components of the energy-momentum tensor takes the simple form with all other first-order contributions nicely cancelling. In the tensor mode, for which q i h ij = h ii = 0, the only non-vanishing perturbations are to the space-space components of the energy-momentum tensor, which are needed in calculating the viscous damping of gravitational waves [4]. These take the form where K(v) ≡ j 2 (v)/v 2 . For scalar modes, we get contributions to all components, of the form in which δp ν , π ν , δu ν and δρ ν may be identified as the pressure perturbation, scalar anisotropic inertia, velocity potential, and energy density perturbation, respectively, for a given neutrino gives It is only in the case m = 0 that the integrals over neutrino energies can be done separately from the integrals over time, and give results simply proportional toρ ν .
The contribution of photons to the perturbations in the energy momentum tensor can be calculated in a similar way, by using Eqs.  Here the source function is calculated for a gravitational perturbation D q (t) whose value D o q before horizon entry is unity; for any other initial condition, Ψ should be multiplied by the value of D o q . The quantity κ is the wave number q in units of the wave number that just comes into the horizon at matter-radiation equality, and y is the Robertson-Walker scale factor, in units of the scale factor at matter-radiation equality.