Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.mrao.cam.ac.uk/~graca/NGsims/ME630rv.ps.gz Äàòà èçìåíåíèÿ: Mon Nov 1 21:09:06 2004 Äàòà èíäåêñèðîâàíèÿ: Sat Dec 22 03:50:44 2007 Êîäèðîâêà: IBM-866 Ïîèñêîâûå ñëîâà: http astrokuban.info astrokuban

Mon. Not. R. Astron. Soc. 000, 1--11 (2004) Printed 28 September 2004 (MN L A T E X style file v2.2)
Simulation of noníGaussian CMB maps
Graca Rocha 1,2,3 , M.P. Hobson 1 , Sarah Smith 1 , Pedro Ferreira 2 and Anthony Challinor 1
1 Astrophysics Group, Cavendish Laboratory, Magingley Road, Cambridge CB3 0HE, UK
2 Astrophysics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK
3 Centro de AstrofÒÐsica da Universidade do Porto, R. das Estrelas s/n, 4150í762 Porto, Portugal
Accepted ---. Received ---; in original form 28 September 2004
ABSTRACT
A simple method is presented for the rapid simulation of statisticallyíisotropic noníGaussian
maps of CMB temperature fluctuations with a given power spectrum and analyticallyí
calculable bispectrum and higheríorder polyspectra. The nthíorder correlators of the pixel
values may also be calculated analytically. The cumulants of the simulated map may be used
to obtain an expression for the probability density function of the pixel temperatures. The
statistical properties of the simulated map are determined by the univariate noníGaussian
distribution from which pixel values are drawn independently in the first stage of the simuí
lation process. We illustrate the method using a noníGaussian distribution derived from the
wavefunctions of the harmonic oscillator. The basic simulation method is easily extended to
produce noníGaussian maps with a given power spectrum and diagonal bispectrum.
Key words: cosmology: cosmic microwave background -- cosmology: theory -- statistics --
numerical simulations
1 INTRODUCTION
The study of noníGaussianity of Cosmic Microwave Background
(CMB) fluctuations is of major importance in understanding the
processes responsible for generating the fluctuations and in assessí
ing the contribution of foreground astrophysical processes and iní
strumental effects to observations of the CMB.
Singleífield inflationary scenarios predict, in general, that
CMB fluctuations are very nearly Gaussian (see e.g. Bartolo et al.
2004 for a recent review), if one assumes that the subíHubbleí
scale quantum fluctuations start off in the ground state (Contaldi,
Bean & Magueijo 1999; Martin, Riazuelo & Sakellarioudou 2000;
Gangui, Martin & Sakellarioudou 2002). In such models, noní
Gaussianity produced during inflation arises predominantly from
the nonílinear nature of gravitational interactions rather than from
selfíinteraction of the fluctuations of the inflaton field. Typically the
level of noníGaussianity is suppressed by the firstíorder slowíroll
parameters (Acquaviva et al. 2003; Maldacena 2003). Subsequent
nonílinear processing of the primordial fluctuations to secondíorder
in perturbation theory has been shown to amplify the tiny primorí
dial noníGaussianity to a level on the lastíscattering surface that
may be detectable with future CMB surveys (Bartolo et al. 2004).
Furthermore, secondíorder radiative transfer effects, such as graví
itational lensing of the CMB, should produce a detectable level of
noníGaussianity in the CMB (Zaldarriaga 2000). Larger levels of
noníGaussianity can be produced in inflation models with multiple
scalar fields. Examples include the curvaton (e.g. Lyth & Wands
2002) and the inhomogeneous reheating (e.g. Dvali, Gruzinov &
Zaldarriaga 2004) scenarios. Finally models which include toí
pological defects also produce significantly noníGaussian fluctuí
ations (Gangui, Pogosian & Winitzki 2002). It is also the case that
inevitable contaminants, such as discrete radio sources, Galactic
emission and systematic instrumental effects leave noníGaussian
signatures on CMB maps. Thus noníGaussianity tests are of fundaí
mental importance both for probing inflation physics and for isolí
ating systematic effects.
In order to investigate one's ability to detect and recover noní
Gaussian signals, it is useful to generate noníGaussian maps with
known statistical properties, to which putative analysis methods
may be applied. In particular, in many applications, it is desirable
that the simulated noníGaussian map is statistically isotropic, with
a prescribed power spectrum (or 2ípoint correlation function). The
generation of such noníGaussian CMB maps is, however, a surí
prisingly difficult task (see, for example, Vio et al. 2001, 2002).
Moreover, it is often useful for the noníGaussian map also to have a
known (or prescribed) bispectrum and oneípoint marginal probabí
ility density function. Several different techniques have been proí
posed to address various subsets of these requirements (Contaldi
& Magueijo 2001; MartÒÐnezíGonzÒalez et al. 2002; Komatsu et al.
2003; Liguori et al. 2003), but each carries a considerable comí
putational cost. The aim of this paper is to present a simple, fast
technique for simulating statisticallyíisotropic noníGaussian CMB
maps with a prescribed power spectrum, for which one can calcuí
late analytically the bispectrum, higheríorder polyspectra, nth orí
der pixel correlators and the oneídimensional marginalised distrií
bution (in terms of its cumulants). Moreover, our basic simulation
method is easily extended to produce noníGaussian maps with a
given power spectrum and diagonal bispectrum.
The problem of simulating noníGaussian CMB maps can be
formalised as follows. A real, random scalar field T (x) can be
c
# 2004 RAS

2 Graca Rocha et al.
defined as a collection of random variables, one at each point
x = (x 1 , x 2 , . . . , x n ) in the nídimensional space (clearly n = 2 for
CMB maps on the sphere or flatípatches of sky). Thus, for each poí
sition x, T (x) # t, where t is a scalar random variable with a oneí
dimensional (marginal) probability density function (PDF) f T (t).
For a random field that is statistically homogeneous, f T (t) is the
same at all points in the space. The main difficulty in the numerical
simulation of a generic random field is that, in general, given two
arbitrary positions x 1 and x 2 , the quantities T (x 1 ) and T (x 2 ) are
not independent. In particular, it is often desirable for T (x) to have
a prescribed 2ípoint covariance function
x T (x 1 , x 2 ) = #T (x 1 )T (x 2 )#.
In the case of a statistically homogeneous random field in Eucí
lidean space, the covariance function depends only on # # x 1 -x 2 .
If the field is also isotropic then the dependence is only on t = |# |.
As discussed by Vio et al. 2001, 2002, most methods for simí
ulating noníGaussian maps with a prescribed 2ípoint covariance
function [and a prescribed marginal PDF f T (t)] are based on first
generating a zeroímean, unitívariance Gaussian random field G(x),
with an appropriate covariance structure x G (t). One then performs
the mapping transformation G(x) # T (x) according to
T (x) = h[G(x)],
where h represents an appropriate function. The usefulness of this
approach derives from the fact that there exist explicit analytical
(but complicated) formulae for the marginal PDF, f T (t), and the
covariance function, x T (t), of the transformed field in terms of
the covariance function, x G (t), of the original Gaussian field and
the mapping function h. In particular, we note that the formula for
x T (t) takes the form of a double integral of a function depending
on both x G (t) and h. There exist only a few mapping functions h
for which f T (t) and x T (t) may be calculated analytically. A still
smaller subset of these cases allows the resulting expressions to be
inverted analytically to obtain the required functions x G (t) and h
to be used in simulating the noníGaussian map. In general, one has
to resort to numerical methods to invert the general formulae for
f T (t) and x T (t) and this can be computationally very costly.
As mentioned above, it is often desirable for the simulated
noníGaussian map also to have a known (or prescribed) bispecí
trum. The method outlined above has not been extended to this
case, and any such generalisation is likely to be extremely comí
putationally demanding. An alternative method for generating noní
Gaussian maps with a prescribed power spectrum and bispectrum
has been suggested by Contaldi & Magueijo (2001), although, in
general, the marginal distribution f T (t) of the resulting map cannot
be obtained analytically. The method is based on choosing some
oneídimensional noníGaussian PDF, from which the real and imaí
ginary parts of (some subset of) the spherical harmonic coefficients
a #m of the map are drawn independently (the remaining coefficients
being drawn from a Gaussian PDF). However, statistical isotropy
imposes `selection rules' upon correlators of the a #m coefficients.
For example, the thirdíorder correlators must satisfy
#a # 1 m 1 a # 2 m 2 a # 3 m 3 # = # # 1 # 2 # 3
m 1 m 2 m 3
# B # 1 # 2 # 3 ,
where (§ § §) denotes the Wigner 3 j symbol and B # 1 # 2 # 3 are the bisí
pectrum coefficients. Hence the map corresponding to the drawn
a #m values is not only noníGaussian but also anisotropic since all
its thirdíorder correlators are zero, except for a subset of the form
#a 3
#m # for those a #m drawn from the noníGaussian PDF. It is thereí
fore necessary first to produce an ensemble of noníGaussian, aní
Figure 1. The noníGaussian PDF used in the simulations (solid line), which
has the form given in equation (A3) with a 3 = 0.2 and s 0 = 1, and a set of
samples drawn from the PDF (histogram).
isotropic maps and then create an isotropic ensemble by applying
a random rotation to each realisation. Contaldi & Magueijo (2001)
show that these random rotations produce the necessary correlaí
tions between the a #m coefficients to ensure isotropy of the ení
semble, but the method clearly requires significant computation.
In this paper, we discuss a very simple and computationally
fast method for simulating noníGaussian maps that are, by coní
struction, statistically isotropic, and for which numerous statistical
properties may be calculated analytically. In particular, it is posí
sible to produce a map with a prescribed power spectrum for which
one can obtain simple analytical expressions for the bispectrum and
higheríorder polyspectra, and nthíorder correlators of the pixel valí
ues. One may also calculate the cumulants of the map, which may
be used to obtain the oneídimensional marginalised distribution. A
simple extension of the method allows for the simulation of maps
with a prescribed power spectrum and diagonal bispectrum. In Secí
tion 2, we discuss our method for simulating noníGaussian maps,
the statistical properties of which are presented in Section 3. In
Section 4 we extend the method to allow simulation of maps with
prescribed power spectrum and diagonal bispectrum. Finally, our
conclusions are presented in Section 5.
2 SIMULATION METHOD
Since our goal is the simulation of CMB maps for use in later anaí
lysis, it is convenient at the outset to divide the celestial sphere
into pixels labelled by p = 1,2, . . . , N pix . For simplicity, we also
assume an equalíarea pixelisation, so that each pixel subtends the
same solid angle W pix = 4p/N pix . Examples of such pixelisation
schemes are HEALPix 1 (GÒorski et al. 1999) and IGLOO (Crittenden
2000). The distribution of pixel centres across the sphere is unimí
portant.
Our simulation method begins by drawing each pixel value
s p # S(x p ) = S(q p , f p ) independently from the same oneí
dimensional noníGaussian PDF f S (s). The precise PDF used is uní
important, but for the purposes of illustration, we adopt here a PDF
derived from the Hilbert space of a linear harmonic oscillator, as
developed by Rocha et al. (2001). This class of PDF is summarí
ised in Appendix A. In particular, we assume the PDF illustrated in
Fig. 1, which is chosen for convenience to have a mean of zero.
The resulting map S(x) will, by construction, be statistically
1 http://www.eso.org/science/healpix/
c
# 2004 RAS, MNRAS 000, 1--11

Simulation of noníGaussian CMB maps 3
isotropic to within the errors introduced by the pixelisation scheme.
In fact, it is a realisation of isotropic noníGaussian white noise,
with a variance given by the second (central) moment ² 2 of the
PDF f S (s). We will assume throughout that the mean ² 1 of the
generating noníGaussian PDF is zero, so that moments and centí
ral moments coincide. According to the `cumulant expansion theí
orem' (Ma 1985) the cumulants, or connected moments, of the mulí
tivariate distribution of pixel values are related to the logarithm of
the momentígenerating function,
M(k) # #e is§k
#, (1)
where s # (s 1 , . . . , s N pix ) is the vector of pixel values and k #
(k 1 , . . . , k N pix
). The connected moments are given in terms of the
derivatives of lnM(k) by
# s p 1 § § § s p n # c = (-i) n ´ n
´k p 1 § § § ´k p n
lnM(k) # # # k=0 . (2)
Given that the pixel values s p are independent random variables we
have ln M(k) = õ p ln M S (k p ), where
lnM S (k) # ln#e isk
# = ln Z å
-å
f S (s)e isk ds
=
å
õ
n=1
(ik) n
n!
k n (3)
is the logarithm of the momentígenerating function (the cumulantí
generating function) for the noníGaussian PDF f S (s) whose cumuí
lants are the k n . Substituting in equation (2), we find
# s p 1 § § § s p n # c = d p 1 §§§ p n
(-i) n d n
dk n ln M S (k) # # # k=0
= d p 1 §§§p n
k n . (4)
The symbol d p 1 §§§p n equals one if all the n pixels are the same and
vanishes otherwise. If required, the pixel correlations can always
be expanded in their connected parts (Ma 1985). In particular, since
#s p # = 0 for each pixel, one finds that, for example,
#s p 1 s p 2 # = #s p 1 s p 2 # c
#s p 1 s p 2 s p 3 # = #s p 1 s p 2 s p 3 # c
#s p 1 s p 2 s p 3 s p 4 # = #s p 1 s p 2 s p 3 s p 4 # c + #s p 1 s p 2 # c #s p 3 s p 4 # c
+#s p 1 s p 3 # c #s p 2 s p 4 # c + #s p 1 s p 4 # c #s p 2 s p 3 # c .
The next step in the simulation procedure is to transform the
map S(x) into sphericalíharmonic space (using, for example, the
map2alm routine from the HEALPix package) to obtain the coeffií
cients
a #m =
N pix
õ
p=1
Y #
#m (x p )s p W pix # Z 4p
Y #
#m (q,f)S(q,f)dW. (5)
Using equation (4) and the (approximate) orthogonality of the
(pixelised) spherical harmonics, we quickly find that the secondí
order correlator of the harmonic coefficients is given by
#a #m a # # # m # # = ² 2 W pix d ## d mm , (6)
where ² 2 = k 2 since f S (s) has vanishing mean. In order to obtain
a final noníGaussian map with a particular prescribed ensembleí
average power spectrum, C # , one then rescales the harmonic coefí
ficients to obtain
Õ
a #m = a #m
# C #
² 2 W pix
, (7)
such that # Õ
a #m Õ
a # # # m # # = C # d ## d mm . We note that the effect of a
Figure 2. A realisation of a noníGaussian allísky CMB map with a preí
scribed ensembleíaverage power spectrum C # , obtained using the noní
Gaussian PDF plotted in Fig. 1 with a 3 = 0.2, HEALPix resolution paraí
meter N side = 512 (WMAP resolution) and # max = 1500.
Figure 3. A realisation of a Gaussian allísky CMB map drawn from the
same ensembleíaverage power spectrum, C # , as the noníGaussian map
shown in Fig. 2.
spatiallyíinvariant, circularlyísymmetric observing beam on the fií
nal map is trivially included by letting C # #C # B 2
# , where B l are the
coefficients of the beam in a Legendre expansion. Finally, the harí
monic coefficients Õ
a #m are inverse spherical harmonic transformed
(using, for example, the alm2map routine from the HEALPix packí
age) to obtain the final map
t p # T (x p ) = õ
#,m
Õ
a #m Y #m (x p ), (8)
where the double summation extends from # = 0 to # max and m =
-# to #. The equivalent flatísky approximation for small patches is
discussed in Appendix C.
In Fig. 2 we plot a realisation of a noníGaussian allísky CMB
map generated as described above, using the noníGaussian PDF
plotted in Fig. 1 with a prescribed ensembleíaverage power specí
trum, C # . The map was produced using the HEALPix pixelisation
scheme with the N side parameter set to 512, which corresponds
to N pix = 12N 2
side # 3 ½ 10 6 equalíarea pixels. For comparison, in
Fig. 3, we plot a realisation of a Gaussian CMB map with the same
power spectrum C # , using the same pixelisation.
The source code to simulate the noníGaussian CMB maps for
c
# 2004 RAS, MNRAS 000, 1--11

4 Graca Rocha et al.
both the full sky and for a small patch of the sky are available at the
NGSIMS webpage 2 .
As a guide to help the reader reproduce our method, we
give here a summary of the logical steps to be taken to generate
these noníGaussian maps (see also documentation at the NGSIMS
webpage):
(i) Draw independent identicallyídistributed pixel values from
a noníGaussian PDF to create a statisticallyíisotropic map of noní
Gaussian white noise;
(ii) Transform the map to harmonic space (with e.g. a fast spherí
ical transform);
(iii) Scale the harmonic coefficients to enforce the desired power
spectrum;
(iv) Transform back to real space to obtain a noníGaussian map
which has the desired 2ípoint correlation function.
We have implemented this method using two classes of noní
Gaussian PDF (see documentation in NGsims webpage), but for
the purposes of illustration, we adopt here PDF (i):
(i) A general PDF based on the energy eigenstates of a linear
harmonic oscillator, which takes the form of a Gaussian multiplied
by a series of Hermite polynomials (see Appendix A). In principle
this expansion can be used to generate any PDF, although if the exí
pansion is truncated then the available values of the relative skewí
ness are constrained.
(ii) The pixel values are drawn from a Gaussian distribution and
then raised to some (even) integer power. This method is very fast
to implement and easily generates distributions with large skewness
and so is useful for checking statistical tools.
As we shall see in the next section, the simulation method
described above enables one to generate noníGaussian maps for
which many of the statistical properties can be calculated analytí
ically. The range of possible correlators and polyspectra are, howí
ever, rather restricted, with the scale dependence of the latter coní
trolled solely by the angular power spectrum (see Section 3.2). It
is, however, straightforward to extend our basic method to allow
the simulation of noníGaussian maps with a much wider range of
statistical properties, as will be discussed in Section 4.
3 STATISTICAL PROPERTIES OF THE MAP
It is clear that, by construction, the noníGaussian map plotted in
Fig. 2 is statistically isotropic, with an ensembleíaverage power
spectrum, C # , corresponding to a generic inflationary cosmological
model. As a check, the power spectrum of the map was calculated
and found to agree with the input spectrum, as shown in Fig. 4.
Owing to the simple manner in which the simulated noníGaussian
map is produced, many more of its statistical properties may be
calculated analytically, such as the correlators of the pixel values,
the bispectrum and higheríorder polyspectra, and the marginalised
probability distribution of the map. In this section, we again coní
centrate on the allísky case; the corresponding discussion for the
flatísky approximation is given in Appendix C.
To facilitate the calculation of the statistical properties of the
map, it is convenient first to express the temperature t p in each pixel
in terms of the initial pixel values {s p } drawn from the original
2 http://www.mrao.cam.ac.uk/#graca/NGsims/
0
1000
2000
3000
4000
5000
6000
7000
8000
0 200 400 600 800 1000 1200 1400
l(l+1)C l
/2p
/
²K
2
l
Power spectrum
Input
Output
Figure 4. The power spectrum of the map shown in Fig. 2 (points) as comí
pared with the input ensembleíaverage spectrum C # (solid line).
noníGaussian PDF f S (s). Using equations (7) and (8), we begin by
writing
t p = 1
# ² 2 W pix
õ
#,m
# C # a #m Y #m (x p )
= # W pix
² 2
õ
p
s p õ
#
# C # õ
m
Y #m (x p )Y #
#m (x p ), (9)
where in the second line we have substituted for a #m from equation
(5). Using the spherical harmonic addition formula
õ m
Y #m (x p )Y #
#m (x p ) = 2#+1
4p P # (x p § x p ), (10)
where P # (z) is a Legendre polynomial, we may thus write equation
(9) as
t p = õ
p
W pp s p , (11)
where the elements of the weight matrix are given by
W pp = # W pix
² 2
õ
#
2#+1
4p # C # P # (x p § x p # ). (12)
It is also useful to express the spherical harmonic coefficients,
Õ
a #m , of the map as linear superposition of the original pixel values
{s p }. From equations (5) and (7), we immediately obtain
Õ
a #m = õ
p #
W #m,p s p , (13)
where the elements of this second weight matrix are simply scaled
spherical harmonics evaluated at the pth pixel and are given by
#
W #m,p = # W pix C #
² 2
Y #
#m (x p ). (14)
Equations (11)--(14) provide the basic expressions from which
the statistical properties of the noníGaussian map may be obtained.
3.1 Correlators of the pixel values
Using equations (4) and (11), a general expression for the connecí
ted nípoint correlator of the pixel values in the noníGaussian map
is given by
#t p 1 § § §t p n # c = õ
q 1 §§§q n
W p 1 q 1 § § §W p n q n #s q 1 § § § s q n # c
c
# 2004 RAS, MNRAS 000, 1--11

Simulation of noníGaussian CMB maps 5
= k n õ
q
W p 1 q § § §W p n q , (15)
where k n is the nth cumulant of f S (s). Inserting the expression (12)
into this result, we obtain
#t p 1 § § §t p n # c = k n
W pix
# W pix
² 2
# n/2
õ
# 1 §§§# n
a # 1 §§§# n I # 1 §§§# n
(x p 1 , . . . , x p n
),(16)
where we have defined the quantities
a # 1 §§§# n #
(2# 1 +1) § § § (2# n +1)
(4p) n # C # 1 § § §C # n , (17)
I # 1 §§§# n
(x p 1 , . . . , x p n
) # õ
q
P # 1
(x p 1 § x q ) § § § P # n
(x p n § x q )W pix . (18)
Taking the continuum limit of the sum over q, we see that an anaí
lytic expression for the nípoint correlator may be obtained by evalí
uating
I # 1 §§§# n
(x p 1 , . . . , x p n
) # Z 4p
P # 1
(x p 1 § x) § § § P # n
(x p n § x)dW. (19)
Clearly, I # 1 §§§# n is invariant under rigid rotations of its vector arguí
ments, and this property is inherited by the nípoint correlator as
required by statistical isotropy. The integral may be related to the ní
polar harmonics with zero total angular momentum (Varshalovich
et al. 1988); this reduction for the first few values of n is described
further in Appendix B.
Using the results derived in Appendix B, one recovers, for exí
ample, the wellíknown result for the 2ípoint correlator
#t p 1 t p 2 # = õ
#
2#+1
4p C # P # (x p 1 § x p 2
), (20)
and an explicit form for the 3ípoint correlator given by
#t p 1 t p 2 t p 3 # c
= ² 3
² 3/2
2
W 1/2
pix õ
# 1 ,# 2 ,# 3
# D # 1 D # 2 D # 3
4p
# # 1 # 2 # 3
0 0 0 #
½ õ
m 1 m 2 m 3
# # 1 # 2 # 3
m 1 m 2 m 3
# Y #
# 1 m 1
(x p 1
)Y #
# 2 m 2
(x p 2
)Y #
# 3 m 3
(x p 3
),
in which we have defined D # # (2# +1)C # .
3.2 Bispectrum and higheríorder polyspectra
To compute the polyspectra of the noníGaussian map we form the
nth connected correlator of its multipoles using equations (4) and
(13):
# Õ
a # 1 m 1 § § § Õ
a # n m n # c = õ
p 1 §§§p n #
W # 1 m 1 ,p 1 § § § #
W # n m n ,p n #s p 1 § § § s p n # c
= k n õ
p #
W # 1 m 1 ,p § § § #
W # n m n ,p . (21)
Inserting the expression (14) into this result, we obtain
# Õ
a # 1 m 1 § § § Õ
a # n m n # c = k n
W pix
# W pix
² 2
# n/2
# C # 1 § § §C # n J # 1 m 1 ,...,# n m n , (22)
where
J # 1 m 1 ,...,# n m n
= õ
p
Y #
# 1 m 1
(x p ) § § §Y #
# n m n
(x p )W pix
# Z 4p
Y #
# 1 m 1
(x) § § §Y #
# n m n
(x)dW. (23)
Integrals of this form may be reduced to products of 3 j symbols as
discussed in Appendix B.
For the case n = 2, one immediately recovers # Õ
a #m Õ
a # # # m # # =
C # d ## # d mm # . For n = 3, we find that
# Õ
a # 1 m 1
Õ
a # 2 m 2
Õ
a # 3 m 3 # = # # 1 # 2 # 3
m 1 m 2 m 3
# B # 1 # 2 # 3 , (24)
where the bispectrum coefficients of the simulated noníGaussian
map are given by
B # 1 # 2 # 3
= k 3
W pix
# W pix
² 2
# 3/2
# D # 1
D # 2
D # 3
4p
# # 1 # 2 # 3
0 0 0 # . (25)
It is convenient also to introduce the `normalised' (by the power
spectrum) reduced bispectrum Ó
b # 1 # 2 # 3 by the relation
B # 1 # 2 # 3
= # D # 1 D # 2 D # 3
4p # # 1 # 2 # 3
0 0 0 # Ó
b # 1 # 2 # 3 . (26)
Thus, for our simulated map, Ó
b # 1 # 2 # 3
has a constant value given by
Ó b # 1 # 2 # 3
= k 3
W pix
# W pix
² 2
# 3/2
. (27)
We see that the amplitude of the bispectrum is determined by
the variance ² 2 and skewness k 3 of the original noníGaussian PDF
f S (s). If f S (s) were Gaussian, for example, the bispectrum would
clearly vanish. Even if f S (s) has a nonízero skewness, however,
we note from equation (22) that, in the limit that the number of
pixels tends to infinity, all polyspectra with n # 3 vanish. This is
a consequence of each pixel being the weighted sum of independí
ent identicallyídistributed variates, so that in the limit of an infinite
number of pixels the processed map T (x) tends to Gaussian by the
centralílimit theorem. Fortunately, in practice, we can bypass this
property as N pix increases by scaling the cumulants of the initial
distribution f S (s). For example, to get a bispectrum independent of
N pix , we must scale k 3 such that
k 3
W pix
# W pix
² 2
# 3/2
= constant. (28)
As mentioned in Appendix A, for the particular noníGaussian PDF
used here, generating higher values of the relative skewness reí
quires one to have a larger range of the generalized cumulants paraí
meters a n non zero.
As a check on our calculations, the normalised reduced bisí
pectra of an ensemble of 300000 noníGaussian maps was calcuí
lated using an estimator defined in Spergel & Goldberg (1999) as
Ó b # 1 # 2 # 3
= 1
4p Z e # 1
(x)e # 2
(x)e # 3
(x)dW
½ # 4p
D # 1
D # 2
D # 3
# # 1 # 2 # 3
0 0 0 # -2
, (29)
where, e # (x) = # 4p
2#+1 õm a #m Y #m (x). The resulting mean values of
Ó b ### are plotted in Fig. 5, together with the associated uncertainties.
Note that it is important to divide the measured bispectrum by the
ensembleíaveraged values of the power spectrum; dividing by the
measured values for each simulation produces a biased estimator.
The predicted value of Ó b # 1 # 2 # 3 was calculated using (27) and is plotí
ted as the solid line in the figure. We see that the measured and
predicted values are fully consistent.
Returning to equation (22), if one considers n = 4 and follows
the notation of Hu (2001), one finds that the connected part of the
c
# 2004 RAS, MNRAS 000, 1--11

6 Graca Rocha et al.
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0 10 20 30 40 50
lll
>
/
C l
3/2
l
Fullísky noníGaussian simulation
Measured values
Theoretical values
Figure 5. Nonízero diagonal components of the bispectrum estimated from
300000 noníGaussian simulations with the same parameters as those used
in Fig. 2. The mean over the simulations and its standard error are plotted.
The dashed line is the theoretical ensembleíaverage value. Note that the
diagonal bispectrum necessarily vanishes for odd # by parity.
trispectrum of the noníGaussian simulation is given by
T # 1 # 2
# 3 # 4
(L) = W pix k 4
² 2
2 # D # 1
D # 2
D # 3
D # 4
2L+1
4p
½ # # 1 # 2 L
0 0 0 ## # 3 # 4 L
0 0 0 # . (30)
Analytic expressions for higheríorder polyspectra may be obtained
in an analogous manner.
3.3 Cumulants and marginal distribution
Finally, we consider the marginal PDF f T (t) of the processed map.
This is mostly easily carried out by considering the cumulants of
f T (t), and relating them to those of the original noníGaussian PDF
f S (s).
The cumulants of the marginal PDF are equal to the connected
parts of the correlators of the pixel temperatures with all pixels in
the correlator the same. From equation (15), we find that
k n = # N pix
õ
p #
W n
pp #
# k n
= k n
W pix
# W pix
² 2
# n/2
õ
# 1 §§§# n
a # 1 §§§# n I # 1 §§§# n
(x p , . . . , x p ), (31)
which is independent of the choice of pixel p. Evaluating the
rotationallyíinvariant function I # 1 §§§# n
(x p , . . . , x p ) along the polar
axis, where Y #m (Óz) = # (2l +1)/4pd m0 , we find that the first few
cumulants (beyond k # 1 = 0) are
k 2 = õ
#
2#+1
4p C #
k 3 = õ
# 1 # 2 # 3
B # 1 # 2 # 3
# 2# 1 +1
4p § § §
2# 3 +1
4p # # 1 # 2 # 3
0 0 0 #
k 4 = õ
# 1 §§§# 4
# 2# 1 +1
4p § § §
2# 4 +1
4p
½ õ
L
T # 1 # 2
# 3 # 4
(L) # # 1 # 2 L
0 0 0 ## # 3 # 4 L
0 0 0 # ,
Figure 6. Histogram of pixel temperatures for the final processed noní
Gaussian map for N side = 64 and # max = 128, created using an initial PDF
with parameters a 3 = 0.27 and s 0 = 1. Also shown are the Edgeworth exí
pansion for f T (t) (dashed line), given by equation (33) but truncated at
n = 4, and a Gaussian with the same variance (solid line).
where we have written the results in such as way that they are true
generally, for any statisticallyíisotropic map.
In principle, knowledge of the complete set of cumulants k # n
may be used to obtain an explicit expression for the marginal PDF
f T (t). This could be carried out, for example, by first obtaining its
momentígenerating function
M T (k) = #e ikt
# = exp
å
õ
n=1
(ik) n
n! k # n
# , (32)
and then performing an inverse Fourier transform to yield f T (t).
Alternatively, for a weakly noníGaussian distribution, one can emí
ploy the Edgeworth expansion. In this approach, the PDF is exí
pressed as an asymptotic expansion around a Gaussian with mean
zero and variance s 2 = k # 2 to yield
f T (t) #
e - Ó t 2
# 2ps 2 # 1+
k # 3 /s 3
12 # 2
H 3 ( Ó
t) + k # 4 /s 4
96 H 4 ( Ó
t)
+ k # 5 /s 5
480 # 2
H 5 ( Ó t) + (k 6 +10k # 3
2 )/s 6
5760 H 6 ( Ó t) + § § § # , (33)
where Ó t = t/ # 2s. When k # n /s n is small for n larger than some iní
teger one can use a finite number of cumulants as an acceptable
approximation to the distribution. As pointed out by Rocha et al.
(2001), however, this might no longer be a PDF and, in particular,
might deviate from the original distribution in the tails. In Fig.6
we plot the histogram of the pixel temperatures for a noníGaussian
map with N side = 64 and # max = 128 (somewhat lower resolution
than that plotted in Fig. 2), generated from an initial PDF with paraí
meters a 3 = 0.27 and s 0 = 1. Overplotted is the Edgeworth expaní
sion of f T (t), equation (33), truncated at n = 4. The cumulants of
f T (t) can be computed efficiently from the first expression in equaí
tion (31); we find k # 1 = 0, k # 2 # 13809²K 2 , k # 3 # 177234²K 3 and
k # 4 # 6307963²K 4 . The Edgeworth expansion agrees with the simí
ulation results better than a Gaussian with the same variance.
4 EXTENDED SIMULATION METHOD
We see from the previous section that our basic simulation method
generates maps with a rather restricted range of possible correlators
and polyspectra, with the scale dependence of the latter controlled
c
# 2004 RAS, MNRAS 000, 1--11

Simulation of noníGaussian CMB maps 7
Figure 7. A realisation of a noníGaussian allísky map simulation using
three bins: B 1 = [0, 500],a 3 = 0.1; B 2 = [501,1000],a 3 = 0.2; and B 3 =
[1001,1500],a 3 = 0.25.
solely by the angular power spectrum. It is straightforward, howí
ever, to extend our basic method to allow the simulation of noní
Gaussian maps with a much wider range of statistical properties. In
particular, there is no fundamental requirement for the method to
be restricted to the same set of cumulants over the whole range of
scales.
The procedure is as follows. One divides the range of mulí
tipoles, #, into noníoverlapping bins. A given bin B will be a set
B = [# B min , # B max
]. For each bin B we simulate a map t Bp , with noní
zero C # only for # # B and a noníGaussian PDF that can differ
between bins. The final map is a superposition of these band maps,
i.e. with pixel values
t p = õ
B
t Bp . (34)
Since each map t Bp is individually statistically isotropic, then so
too is their sum. Moreover, from equation (26), we see that the bisí
pectrum for each band map can only be nonízero within the corresí
ponding bin (this is also true for the connected parts of higheríorder
polyspectra). The nth cumulant of the final map is also simply the
sum of the nth cumulants of the individual band maps. With this
method we are thus able to generate maps with more general statí
istical properties.
For example, by choosing appropriate values of k 3 for the noní
Gaussian PDF used to simulate each band map, one can arrange for
the summed map (34) to have a given power spectrum C # and an
arbitrary prescribed constant value of the reduced normalised bisí
pectrum Ó b # 1 # 2 # 3
in each bin. As an illustration, in Fig. 7 we plot a
noníGaussian map generated using three bins. In each bin, the noní
Gaussian PDF used was of the form given in equation (A3) with
s 0 = 1. However, the values of a 3 used in each bin were a 3 = 0.1,
0.2 and 0.25 respectively. As a check on our calculations, the norí
malised reduced bispectrum of an ensemble of 300000 such noní
Gaussian maps was calculated using the estimator (29). The resí
ulting mean values of Ó b ### , for individual values of #, are plotted
in Fig. 8, together with the associated uncertainties. The predicted
value of Ó
b ### in each of the three broad bins was calculated using
equation (27) and are plotted as the dashed lines in the figure. We
see that once again the measured and predicted values are fully
consistent.
How general is our extended method? We can in principle efí
ficiently generate maps with arbitrary diagonal bispectra, i.e. with
any given C # and B ### . This is of some importance because curí
0
0.005
0.01
0.015
0.02
0.025
0.03
0 10 20 30 40 50
lll
>
/
C l
3/2
l
Fullísky noníGaussian simulation using 3 bins
Measured values
Theoretical values
Figure 8. Nonízero diagonal components of the bispectrum estimated from
300000 noníGaussian simulations using three bins: B 1 = [0, 21],a 3 = 0.1;
B 2 = [22,43],a 3 = 0.2; and B 3 = [44,56],a 3 = 0.25. The mean over the
simulations and its standard error are plotted. The dashed lines shows the
theoretical ensembleíaverage value in each bin.
rently known primordial theories of noníGaussianity lead to more
general combinations of angular power spectra and bispectra than
can be created from our method with a single univariate PDF f S (s).
It is not possible, however, to generate specific models of primí
ordial noníGaussianity exactly with our extended method. To do
that one would have to be able to choose arbitrary values for the
angular spectrum and for all components of the bispectrum (and
higher polyspectra). This would involve going beyond the simple
oneípoint PDF methods advocated here.
5 CONCLUSIONS
We presented a simple, fast method for simulating statisticallyí
isotropic noníGaussian CMB maps with a given power spectrum
and analytically calculable bispectrum. We showed that our techí
nique allows one to describe the statistical properties of the map
by computing analytically the nthíorder polyspectra, and the nthí
order correlators of the pixel values. We showed that these can be
expressed in terms of the nípolar harmonics with zero total angular
momentum, and we describe this reduction for the first few valí
ues of n. We also recovered analytically the oneídimensional marí
ginalised distribution function in terms of its cumulants. The unií
variate noníGaussian distribution, from which the pixel values are
drawn independently in the first stage of the simulation process,
fully determines the statistical properties of the final map. Here we
used a noníGaussian distribution derived from the wavefunctions
of the harmonic oscillator. Simulations of both the full sky and a
small patch of the sky were generated and corresponding statistical
analysis performed. As a check on our calculations we computed
both the power spectrum and bispectrum of the simulated maps and
found them to be fully consistent. The simulation method described
here clearly enables one to generate maps with wellídefined correlí
ators and polyspectra. We extended the method to encompass difí
ferent set of cumulants over the whole range of scales, generating
maps with arbitrary power spectra and diagonal bispectra for differí
ent scales. It is not possible, however, to generate specific models
of primordial noníGaussianity exactly with our extended method,
since these require offídiagonal bispectrum coefficients to be speí
cified arbitrarily. This would involve going beyond the simple oneí
point PDF methods advocated here.
c
# 2004 RAS, MNRAS 000, 1--11

8 Graca Rocha et al.
The source code to simulate the noníGaussian CMB maps for
both the full sky and for a small patch of the sky are available at the
NGSIMS webpage 3 .
A pertinent question is what other statistical properties can
be calculated analytically for the class of noníGaussian maps we
have investigated. Of particular interest are the phase associations
between different harmonic coefficients. In Matsubara (2003), a
general relationship between phase correlations and the hierarchy
of polyspectra in Fourier space is established. It is also stated that
the phase correlations are related to the polyspectra through the
noníuniform distribution of the phase sum q k 1
+ q k 2
+ § § § + q kN
with closed vectors k 1 +k 2 + § § § +k N = 0. We are currently investí
igating the form of the distribution function of this phase sum in our
maps. A study of the Minkowski functionals of our noníGaussian
maps is also underway.
6 ACKOWLEDGEMENTS
We thank Carlo Contaldi and Neil Turok for useful discussions, and
Martin Kunz and Grazia De Troia for providing us with the bispecí
trum code for the fullísky case. Some of the results in this paper
have been derived using the HEALPix (Gorski, Hivon, and Wandelt
1999) package. GR acknowledges a Leverhulme Fellowship at the
University of Cambridge. PF and AC acknowledge Royal Society
University Research Fellowships. SS acknowledges support by a
PPARC studentship.
REFERENCES
Acquaviva V., Bartolo N., Matarrese S., Riotto A., 2003, Nuclear
Physics B, 667, 119
Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004, Physics
Reports
Bartolo N., Matarrese S., Riotto A., 2004, J. High Energy Phys.,
04, 006
Contaldi C., Bean R., Magueijo J., 1999, Phys. Lett., B468, 189
Contaldi C., Magueijo J., 2001, Phys. Rev., D63, 103512
Crittenden R. G., 2000, Astrophysical Letters and Communicaí
tions, 37, 377
Dvali G., Gruzinov A., Zaldarriaga M., 2004, Phys. Rev. D69,
023505
Edmonds A. R., 1974, Angular Momentum in Quantum Mechaní
ics. Princeton University Press, Princeton, New Jersey
Gangui A., Martin J., Sakellarioudou M., 2002, Phys. Rev. D, 66,
083502
Gangui A., Pogosian L., Winitzki S., 2002, New Astronomy Reí
views, 46, 681
GÒorski K. M., Hivon E., Wandelt B. D., 1999, in Banday A. J.,
Sheth R. S., Costa L. D., eds, Proceedings of the MPA/ESO Cosí
mology Conference `Evolution of LargeíScale Structure' Printí
Partners Ipskamp, NL, pp 37--42
Hu W., 2001, Phys. Rev., D64, 083005
Komatsu E., et al., 2003, ApJS, 148, 135H
Liguori M., Matarrese S., Moscardini L., 2003, ApJ., 597, 57
Lyth D. H., Wands D., 2002, Phys. Lett. B, 524, 5
Ma S. K., 1985, Statistical Mechanics. World Scientific, Philí
adelphia
Maldacena J., 2003, J. High Energy Phys., 05, 013
3 http://www.mrao.cam.ac.uk/#graca/NGsims/
Martin J., Riazuelo A., Sakellarioudou M., 2000, Phys. Rev. D,
61, 083518
MartÒÐnezíGonzÒalez E., Gallegos J., Argueso F., CayÒon L., Sanz J.,
2002, MNRAS, 336, 22
Matsubara T., 2003, ApJ., 591, L79
Rocha G., Magueijo J., Hobson M., Lasenby A., 2001, Phys. Rev.,
D64, 063512
Smith S., et al., 2004, MNRAS, 352, 887
Spergel D. N., Goldberg D. M., 1999, Phys. Rev., D59, 103001
Varshalovich D. A., Moskalev A. N., Khersonskii V. K., 1988,
Quantum Theory of Angular Momentum. World Scientific,
Singapore
Vio R., Andeani P., Tenorio L., Wamsteker W., 2001, PASP, 113,
1009
Vio R., Andeani P., Tenorio L., Wamsteker W., 2002, PASP, 114,
1281
Zaldarriaga M., 2000, Phys. Rev. D62, 063510
This paper has been typeset from a T E X/ L A T E X file prepared by the
author.
c
# 2004 RAS, MNRAS 000, 1--11

Simulation of noníGaussian CMB maps 9
APPENDIX A: NONíGAUSSIAN PDFS BASED ON THE HARMONIC OSCILLATOR
In this appendix, we summarise the class of probability distribution functions (PDFs) derived from the Hilbert space of a linear harmonic
oscillator, which was developed by Rocha et al. (2001). The original noníGaussian distribution, f S (s), used in the main text to produce the
simulated noníGaussian maps is an example of such a PDF.
This general PDF is based on the coordinateíspace wavefunctions of the energy eigenstates of a linear harmonic oscillator, and takes
the form of a Gaussian multiplied by the square of a (possibly finite) series of Hermite polynomials whose coefficients a n are used as
noníGaussian qualifiers. In particular, if x is a general random variable, the most general PDF has the form
p(x) = |y| 2 = e -x 2 /(2s 2
0 ) # # # # õ n
a n C n H n # x
# 2s 0
## # # #
2
, (A1)
where H n (x) are the Hermite polynomials, and the quantity s 2
0 is the variance associated with the (Gaussian) probability distribution for the
ground state |y 0 | 2 . The constants C n are fixed by normalising the individual states. The only constraint upon the amplitudes a n is
õ |a n | 2 = 1. (A2)
This is a simple algebraic expression which can be eliminated explicitly by writing a 0 = # 1-õ å
1 |a n | 2 . Thus the coefficients a n can be
independently set to zero without mathematical inconsistency (Rocha et al. 2001). Moreover, these coefficients can be written as series of
cumulants (Contaldi, Bean & Magueijo 1999) and should indeed be regarded as noníperturbative generalisations of cumulants.
For the simulations in the main text, we use the noníGaussian PDF for which all a n are set to zero, except for the real part of a 3 (and
consequently a 0 ). The reason for this choice is that this quantity reduces to the skewness in the perturbative regime. The imaginary part of
a 3 is only meaningful in the noníperturbative regime (and can be set to zero independently without inconsistency). Hence we consider a PDF
of the form
p(x) = e -x 2 /(2s 2
0 )
# 2ps 0
# a 0 + a 3
# 48
H 3 # x
# 2s 0
## 2
, (A3)
with a 0 = # 1-a 2
3 . It is straightforward to show that the first, second and third moments of our PDF are related to a 3 and s 0 by (Contaldi
& Magueijo 2001)
² 1 = 0
² 2 = s 2
0 # 1+6a 2
3 #
² 3 = # 2s 2
0 # 3
2
# 3
# a 2
3 # 1-a 2
3 ## . (A4)
The PDF therefore has zero mean and a fixed variance and skewness. In the simulations discussed in the main text, we choose a 3 = 0.2 and
s 0 = 1. This resulting PDF is plotted in Fig. 1.
We note that the space of possible PDFs is constrained as a result of restricting the set of coefficients a n to two nonízero values. This
implies that we cannot generate distributions with arbitrarily large relative skewness. Indeed, ² 3 /² 3/2
2 is bounded above by 0.74, and takes
this maximum value for a 2
3 = (7 - # 43)/6 = 0.27 2 . However, in general our method can generate higher values of the relative skewness
(since it can generate any distribution) but for that purpose one needs more nonízero coefficients a n (Contaldi & Magueijo 2001).
APPENDIX B: SOME USEFUL INTEGRALS
We give here useful results concerning integrals involving products of Legendre polynomials:
I # 1 §§§# n
(x 1 , . . . , x n ) # Z 4p
P # 1
(x 1 § x) § § § P # n
(x n § x)dW. (B1)
Using the addition theorem (equation 10), this reduces to evaluating integrals of products of spherical harmonics,
J # 1 m 1 ,...,# n m n # Z 4p
Y # 1 m 1
(x) § § §Y # n m n
(x)dW (B2)
since
I # 1 ,§§§,# n
(x 1 , . . . , x n ) = 4p
2# 1 +1 § § §
4p
2# n +1 õ
m 1 §§§m n
Y #
# 1 m 1
(x 1 ) § § §Y #
# n m n
(x n )J # 1 m 1 ,...,# n m n
. (B3)
First, consider the case n = 2. Using the orthonormality of the spherical harmonics, and the relation Y #
#m (x) = (-1) m Y #-m (x), to evaluate
J # 1 m 1 ,# 2 m 2
= (-1) m 1 d # 1 # 2
dm 1 -m 2 , and then applying the addition theorem we find the wellíknown result
I # 1 # 2
(x 1 , x 2 ) = 4p
2# 1 +1 P # 1
(x 1 § x 2 )d # 1 # 2 . (B4)
For integrals involving products of three or more spherical harmonics, the general strategy is to combine pairs of harmonics using the
ClebschíGordan series (e.g. Varshalovich et al. 1988; Edmonds 1974)
c
# 2004 RAS, MNRAS 000, 1--11

10 Graca Rocha et al.
Y # 1 m 1
(x)Y # 2 m 2
(x) = õ
#m
# (2# 1 +1)(2# 2 +1)(2#+1)
4p
# # 1 # 2 #
m 1 m 2 m ## # 1 # 2 #
0 0 0 # Y #
#m (x), (B5)
until we have only a single pair left which can then be integrated trivially using orthonormality. We illustrate this procedure for the case
of n = 3 and n = 4 since these are needed for the calculation of the bispectrum and trispectrum. For n = 3, from equation (B5) we find
immediately that
J # 1 m 1 ,...,# 3 m 3
= # (2# 1 +1)(2# 2 +1)(2# 3 +1)
4p
# # 1 # 2 # 3
0 0 0 ## # 1 # 2 # 3
m 1 m 2 m 3
# , (B6)
and therefore
I # 1 # 2 # 3
(x 1 , x 2 , x 3 ) = (4p) 2 # 4p
(2# 1 +1)(2# 2 +1)(2# 3 +1) # # 1 # 2 # 3
0 0 0 # õ
m 1 m 2 m 3
# # 1 # 2 # 3
m 1 m 2 m 3
# Y #
# 1 m 1
(x 1 )Y #
# 2 m 2
(x 2 )Y #
# 3 m 3
(x 3 ). (B7)
The final term in this equation (the summation over m 1 , m 2 and m 3 ) ensures that I # 1 # 2 # 3 is invariant under rigid rotations of its vector
arguments x 1 , x 2 and x 3 . As expected, the summation can be expressed in terms of the tripolar spherical harmonics with zero total angular
momentum (Varshalovich et al. 1988).
Consider now the case n = 4. There is now some freedom in the choice of spherical harmonics to combine. If we couple Y # 1 m 1 with Y # 2 m 2
and Y # 3 m 3 with Y # 4 m 4 , we find
J # 1 m 1 ,...,# 4 m 4
= õ
#m
(-1) m # (2# 1 +1)(2# 2 +1)(2#+1)
4p
# (2# 3 +1)(2# 4 +1)(2#+1)
4p (B8)
½ # # 1 # 2 #
m 1 m 2 m ## # 3 # 4 #
m 3 m 4 -m
## # 1 # 2 #
0 0 0 ## # 3 # 4 #
0 0 0 # . (B9)
The expression on the right is not manifestly symmetric with respect to interchange of e.g. (# 1 m 1 ) and (# 3 m 3 ) since the latter involves a differí
ent coupling scheme. However, the symmetry is easily verified by switching between the two schemes with the 6 j coefficients (Varshalovich
et al. 1988). Finally, we find that
I # 1 §§§# 4
(x 1 , . . . , x 4 ) = (4p) 3/2 # (2# 1 +1) § § § (2# 4 +1)
4p õ
l
2l +1
4p # # 1 # 2 #
0 0 0 ## # 3 # 4 #
0 0 0 #
½ õ
mm 1 §§§m 4
(-1) m # # 1 # 2 #
m 1 m 2 m ## # 3 # 4 #
m 3 m 4 -m # Y #
# 1 m 1
(x 1 )Y #
# 2 m 2
(x 2 )Y #
# 3 m 3
(x 3 )Y #
# 4 m 4
(x 4 ). (B10)
It is straightforward to verify that the last term on the right (the summation over m, m 1 , . . . , m 4 ) is invariant under rigid rotations of x 1 , . . . , x 4 .
APPENDIX C: FLATíSKY APPROXIMATION
For analysis over a small patch of the sky we can use the flatísky approximation and replace spherical transforms by Fourier transforms.
Our starting point is again a pixelised map of noníGaussian white noise, with each pixel value drawn from the noníGaussian PDF f S (s).
Approximating the Fourier transform a(### # #) by a discrete Fourier transform we have
a(### # #) = Z d 2 x
2p
S(x)e -i### # #§x
#
W pix
2p
õ
p
s p e -i### # #§x p , (C1)
where W pix is the pixel area. We evaluate a(### # #) on a regular grid in Fourier space with a Fast Fourier Transform. For a square patch of sky
with N pix pixels, the cell size in Fourier space is (2p) 2 /(N pix W pix ). The secondíorder correlator of the discrete a(### # #) evaluates to
#a(### # #)a # (### # # # )# = # W pix
2p # 2
N pix ² 2 d ### # ## # # # # , (C2)
where ² 2 is the variance of the zeroímean f S (s). In the continuum limit, equation (C2) becomes
#a(### # #)a # (### # # # )# =W pix ² 2 d(### # # -### # # # ), (C3)
where we have used
d ### # ## # # # # # = 1
N pix
õ p
e i(### # #-# # # # # )§x p
#
1
N pix W pix
Z d 2 xe i(### # #-# # # # # # )§x = (2p) 2
N pix W pix
d(### # # -### # # # ). (C4)
We scale the a(### # #) defining Õ
a(### # #) # # C # /(W pix ² 2 )a(### # #) (as for the fullísky case described in the main text) such that the Õ
a(### # #) have the required
power spectrum:
# Õ
a(### # #) Õ
a # (### # # # )# =C # d(### # # -### # # # ). (C5)
Finally, we inverse Fourier transform to obtain our noníGaussian map, T (x), with the prescribed twoípoint statistics:
c
# 2004 RAS, MNRAS 000, 1--11

Simulation of noníGaussian CMB maps 11
T (x p ) = 2p
N pix W pix
õ
### # #
Õ
a(### # #)e i### # #§x p
# Z d 2 ### # #
2p
Õ
a(### # #)e i### # #§x p . (C6)
As in the fullísky case, we can express the pixel values, t p , in the final map as linear combinations of those in the original map, s p :
t p = õ
p
W pp s p , (C7)
where in the continuum approximation
W pp = 1
2p
# W pix
² 2
Z d 2 ### # #
2p # C # e i### # #§(x p -x p # ) (C8)
= # W pix
² 2
Z #d#
2p # C # J 0 (#|x p -x p |),
with J 0 (z) the Bessel function of order zero. Using the asymptotic result J 0 (#q) # P # (cosq), it is straightforward to see that W pp obtained
here in the flatísky limit is equivalent to the fullísky expression (equation 12). Forming the connected nípoint function, as in equation (15),
we find
#t p 1 § § §t p n # c = k n
W pix
# W pix
² 2
# n/2
Z d 2 ### # # 1
(2p) 2 § § §
d 2 ### # # n
(2p) 2 # C # 1 § § §C # n Z d 2 xe i### # # 1 §(x p 1 -x) § § § e i### # # n §(x pn -x)
= (2p) 2 k n
W pix
# W pix
² 2
# n/2
Z d 2 ### # # 1
(2p) 2 § § §
d 2 ### # # n
(2p) 2 # C # 1 § § §C # n d(### # # 1 + § § § +### # # n )e i(### # # 1 §x p 1 +§§§+### # # n §x pn ) . (C9)
For example, for n = 2 we have
#t p 1 t p 2 # c = Z d 2 ### # #
(2p) 2 C # e i### # #§(x p 1 -x p 2 )
# Z #d#
2p C # J 0 (#|x p 1 -x p 2 |). (C10)
In the asymptotic limit this result reduces to the fullísky expression (equation 20).
Finally, we consider the polyspectra of the processed noníGaussian maps. Expressing the Fourier transform of the final map, Õ
a(### # #), in
terms of the original map, i.e.
Õ
a(### # #) = 1
2p
# C # W pix
² 2
õ
p
s p e -i### # #§x p , (C11)
we have
# Õ
a(### # # 1 ) § § § Õ
a(### # # n )# c = k n
(2p) n # W pix
² 2
# n/2
# C # 1 § § §C # n
õ
p
e -i(### # # 1 +§§§+### # # n )§x p
#
k n
W pix
# W pix
² 2
# n/2
(2p) 2-n
# C # 1 § § §C # n
d(### # # 1 + § § § +### # # n ), (C12)
where in the last line we have taken the continuum approximation. This form for the correlator is clearly consistent with rotational, translaí
tional and parity invariance.
As in the fullísky case, we have produced simulated noníGaussian maps and calculated their power spectra and bispectra; the latter were
estimated using the code described in Smith et al. (2004). We find that the flatísky simulations behave as expected with no discernible bias.
This paper has been typeset from a T E X/ L A T E X file prepared by the author.
c
# 2004 RAS, MNRAS 000, 1--11