5.2. CMB bispectrum from active models

Different cosmological models differ in their predictions for the statistical distribution of the anisotropies beyond the power spectrum. Future MAP and Planck satellite missions will provide high-precision data allowing definite estimates of non-Gaussian signals in the CMB. It is therefore important to know precisely which are the predictions of all candidate models for the statistical quantities that will be extracted from the new data and identify their specific signatures.

Of the available non-Gaussian statistics, the CMB bispectrum, or the three-point function of Fourier components of the temperature anisotropy, has been perhaps the one best studied in the literature [Gangui & Martin, 2000a]. There are a few cases where the bispectrum may be deduced analytically from the underlying model. The bispectrum can be estimated from simulated CMB sky maps; however, computing a large number of full-sky maps resulting from defects is a much more demanding task. Recently, a precise numerical code to compute it, not using CMB maps and similar to the CMBFAST code ⁽¹⁷⁾ for the power spectrum, was developed in [Gangui, Pogosian & Winitzki 2001b]. What follows below is an account of this work.

In a few words, given a suitable model, one can generate a statistical ensemble of realizations of defect matter perturbations. We used a modified Boltzmann code based on CMBFAST to compute the effect of these perturbations on the CMB and found the bispectrum estimator for a given realization of sources. We then performed statistical averaging over the ensemble of realizations to compute the expected CMB bispectrum. (The CMB power spectrum was also obtained as a byproduct.) As a first application, we then computed the expected CMB bispectrum from a model of simulated string networks first introduced by Albrecht et al. [1997] and further developed in [Pogosian & Vachaspati, 1999] and in [Gangui, Pogosian & Winitzki 2001].

We assume that, given a model of active perturbations, such as a string simulation, we can calculate the energy-momentum tensor T_µ(x, ) for a particular realization of the sources in a finite spatial volume V₀. Here, x is a 3-dimensional coordinate and is the cosmic time. Many simulations are run to obtain an ensemble of random realizations of sources with statistical properties appropriate for the given model. The spatial Fourier decomposition of T_µ can be written as

(85)

where k are discrete. If V₀ is sufficiently large we can approximate the summation by the integral

(86)

and the corresponding inverse Fourier transform will be

(87)

Of course, the final results, such as the CMB power spectrum or bispectrum, do not depend on the choice of V₀. To ensure this independence, we shall keep V₀ in all expressions where it appears below.

It is conventional to expand the temperature fluctuations over the basis of spherical harmonics,

(88)

where is a unit vector. The coefficients a_lm can be decomposed into Fourier modes,

(89)

Given the sources _µ(k, ), the quantities _l(k) are found by solving linearized Einstein-Boltzmann equations and integrating along the line of sight, using a code similar to CMBFAST [Seljak & Zaldarriaga, 1996]. This standard procedure can be written symbolically as the action of a linear operator _l^µ(k) on the source energy-momentum tensor, _l(k) = _l^µ(k) _µ(k, ), so the third moment of _l(k) is linearly related to the three-point correlator of _µ(k, ). Below we consider the quantities _l(k), corresponding to a set of realizations of active sources, as given. The numerical procedure for computing _l(k) was developed in [Albrecht et al. 1997] and in [Pogosian & Vachaspati, 1999].

The third moment of a_lm, namely <a_{l1 m1} a_{l2 m2} a_{l3 m3}> , can be expressed as

(90)

A straightforward numerical evaluation of Eq. (90) from given sources _l(k) is prohibitively difficult, because it involves too many integrations of oscillating functions. However, we shall be able to reduce the computation to integrations over scalars [a similar method was employed in Komatsu & Spergel, 2001 and in Wang & Kamionkowski, 2000]. Due to homogeneity, the 3-point function vanishes unless the triangle constraint is satisfied,

(91)

We may write

(92)

where the three-point function P_l1l2l3(k₁, k₂, k₃) is defined only for values of k_i that satisfy Eq. (91). Given the scalar values k₁, k₂, k₃, there is a unique (up to an overall rotation) triplet of directions _i for which the RHS of Eq. (92) does not vanish. The quantity P_l1l2l3(k₁, k₂, k₃) is invariant under an overall rotation of all three vectors k_i and therefore may be equivalently represented by a function of scalar values k₁, k₂, k₃, while preserving all angular information. Hence, we can rewrite Eq. (92) as

(93)

Then, using the simulation volume V₀ explicitly, we have

(94)

Given an arbitrary direction ₁ and the magnitudes k₁, k₂ and k₃, the directions ₂ and ₃ are specified up to overall rotations by the triangle constraint. Therefore, both sides of Eq. (94) are functions of scalar k_i only. The expression on the RHS of Eq. (94) is evaluated numerically by averaging over different realizations of the sources and over permissible directions _i; below we shall give more details of the procedure.

Substituting Eqs. (93) and (94) into (90), Fourier transforming the Dirac delta and using the Rayleigh identity, we can perform all angular integrations analytically and obtain a compact form for the third moment,

(95)

where, denoting the Wigner 3j-symbol by , we have

(96)

and where we have defined the auxiliary quantities b_l1l2l3 using spherical Bessel functions j_l,

(97)

The volume factor V₀³ contained in this expression is correct: as shown in the next section, each term _l includes a factor V₀^-2/3, while the average quantity P_l1l2l3(k₁, k₂, k₃) V₀^-3 [cf. Eq. (94)], so that the arbitrary volume V₀ of the simulation cancels.

Our proposed numerical procedure therefore consists of computing the RHS of Eq. (95) by evaluating the necessary integrals. For fixed {l₁l₂l₃} , computation of the quantities b_l1l2l3(r) is a triple integral over scalar k_i defined by Eq. (97); it is followed by a fourth scalar integral over r [Eq. (95)]. We also need to average over many realizations of sources to obtain P_l1l2l3 (k₁,k₂,k₃) . It was not feasible for us to precompute the values P_l1l2l3 (k₁, k₂, k₃) on a grid before integration because of the large volume of data: for each set {l₁l₂l₃} the grid must contain ~ 10³ points for each k_i. Instead, we precompute _l (k) from one realization of sources and evaluate the RHS of Eq. (94) on that data as an estimator of P_l1l2l3 (k₁,k₂,k₃) , averaging over allowed directions of _i. The result is used for integration in Eq. (97).

Because of isotropy and since the allowed sets of directions _i are planar, it is enough to restrict the numerical calculation to directions _i within a fixed two-dimensional plane. This significantly reduces the amount of computations and data storage, since _l (k) only needs to be stored on a two-dimensional grid of k.

In estimating P_l1l2l3 (k₁, k₂, k₃) from Eq. (94), averaging over directions of _i plays a similar role to ensemble averaging over source realizations. Therefore if the number of directions is large enough (we used 720 for cosmic strings), only a moderate number of different source realizations is needed. The main numerical difficulty is the highly oscillating nature of the function b_l1l2l3(r). The calculation of the bispectrum for cosmic strings presented in the next Section requires about 20 days of a single-CPU workstation time per realization.

We note that this method is specific for the bispectrum and cannot be applied to compute higher-order correlations. The reason is that higher-order correlations involve configurations of vectors k_i that are not described by scalar values k_i and not restricted to a plane. For instance, a computation of a 4-point function would involve integration of highly oscillating functions over four vectors k_i which is computationally infeasible.

From Eq. (95) we derive the CMB angular bispectrum _l1l2l3, defined as [Gangui & Martin, 2000b]

(98)

The presence of the 3 j-symbol guarantees that the third moment vanishes unless m₁ + m₂ + m₃ = 0 and the l_i indices satisfy the triangle rule |l_i - l_j| l_k l_i + l_j. Invariance under spatial inversions of the three-point correlation function implies the additional `selection rule' l₁ + l₂ + l₃ = even, in order for the third moment not to vanish. Finally, from this last relation and using standard properties of the 3 j-symbols, it follows that the angular bispectrum _l1l2l3 is left unchanged under any arbitrary permutation of the indices l_i.

In what follows we will restrict our calculations to the angular bispectrum C_l1l2l3 in the `diagonal' case, i.e. l₁ = l₂ = l₃ = l. This is a representative case and, in fact, the one most frequently considered in the literature. Plots of the power spectrum are usually done in terms of l(l + 1)C_l which, apart from constant factors, is the contribution to the mean squared anisotropy of temperature fluctuations per unit logarithmic interval of l. In full analogy with this, the relevant quantity to work with in the case of the bispectrum is

(99)

For large values of the multipole index l, G_lll l^3/2 C_lll. Note also what happens with the 3 j-symbols appearing in the definition of the coefficients _l1l2l3^m₁m₂m₃: the symbol is absent from the definition of C_l1l2l3, while in Eq. (99) the symbol is squared. Hence, there are no remnant oscillations due to the alternating sign of .

However, even more important than the value of C_lll itself is the relation between the bispectrum and the cosmic variance associated with it. In fact, it is their comparison that tells us about the observability `in principle' of the non-Gaussian signal. The cosmic variance constitutes a theoretical uncertainty for all observable quantities and comes about due to the fact of having just one realization of the stochastic process, in our case, the CMB sky [Scaramella & Vittorio, 1991].

The way to proceed is to employ an estimator _l1l2l3 for the bispectrum and compute the variance from it. By choosing an unbiased estimator we ensure it satisfies C_l1l2l3 = <_l1l2l3>. However, this condition does not isolate a unique estimator. The proper way to select the best unbiased estimator is to compute the variances of all candidates and choose the one with the smallest value. The estimator with this property was computed in [Gangui & Martin, 2000b] and is

(100)

The variance of this estimator, assuming a mildly non-Gaussian distribution, can be expressed in terms of the angular power spectrum C_l as follows

(101)

The theoretical signal-to-noise ratio for the bispectrum is then given by

(102)

In turn, for the diagonal case l₁ = l₂ = l₃ = l we have

(103)

Incorporating all the specifics of the particular experiment, such as sky coverage, angular resolution, etc., will allow us to give an estimate of the particular non-Gaussian signature associated with a given active source and, if observable, indicate the appropriate range of multipole l's where it is best to look for it.

¹⁷ http://physics.nyu.edu/matiasz/CMBFAST/cmbfast.html Back.