**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
^{(17)} 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*_{0}. 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*_{0} 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*_{0}. To ensure this
independence,
we shall keep *V*_{0} 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**_{1},
**k**_{2}, **k**_{3}) is defined only
for values of **k**_{i} that satisfy
Eq. (91). Given the scalar values *k*_{1},
*k*_{2}, *k*_{3}, 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**_{1},
**k**_{2}, **k**_{3}) 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*_{1}, *k*_{2}, *k*_{3}, while
preserving all angular
information. Hence, we can rewrite Eq. (92) as

(93) |

Then, using the simulation volume *V*_{0} explicitly, we have

(94) |

Given an arbitrary direction
_{1} and the
magnitudes
*k*_{1}, *k*_{2} and *k*_{3}, the
directions
_{2}
and _{3}
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 3*j*-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*_{0}^{3} contained in this
expression is correct: as shown in the next section, each term
_{l}
includes a factor
*V*_{0}^{-2/3}, while the average quantity
*P*_{l1l2l3}(*k*_{1},
*k*_{2}, *k*_{3})
*V*_{0}^{-3}
[cf. Eq. (94)], so that the arbitrary volume *V*_{0} 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*_{1}*l*_{2}*l*_{3}} ,
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*_{1},*k*_{2},*k*_{3}) .
It was not feasible for us to precompute the
values
*P*_{l1l2l3}
(*k*_{1}, *k*_{2}, *k*_{3}) on a
grid
before integration because of the large volume of data: for each set
{*l*_{1}*l*_{2}*l*_{3}} the grid
must contain ~ 10^{3}
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*_{1},*k*_{2},*k*_{3}) ,
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*_{1}, *k*_{2}, *k*_{3})
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*_{1} + *m*_{2} +
*m*_{3} = 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*_{1} + *l*_{2} +
*l*_{3} = 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*_{1} = *l*_{2} =
*l*_{3} = *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}^{m1m2m3}:
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*_{1} =
*l*_{2} = *l*_{3} = *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.

^{17}
http://physics.nyu.edu/matiasz/CMBFAST/cmbfast.html
Back.