Gaussian Fields

Introduction

Note: To view full latex expression on mobile, you may need to visit the desktop site.

In Cosmology, we are often interested in understanding the spatial correlations amongst fluctuations of various kinds. Of particular importance are two-point correlation functions of the form $$\langle \phi(\bfx)\phi(\bfy)\rangle$$ where $$\phi(\bfx)$$ is some scalar quantity, say a temperature fluctuation.

Because of the homogeneous and isotropic nature of spatial slices, it’s convenient to discuss the Fourier transform of the correlator which takes on the form $$\langle\phi_{\bfk}\phi_{-\bfk}\rangle'=P(k)$$ where the power spectrum, $$P(k)$$, only depends on the magnitude of $$\bfk$$, as indicated (see below for more on the notation).

All of the physics of the two-point function is encoded in $$P(k)$$ and hence it’s typically the quantity we focus on. However, if handed a $$P(k)$$, it’s not very easy to visualize what a corresponding spatial map of fluctuations would look like. In order to gain intuition for the power spectrum and to understand how the fluctuations change as $$P(k)$$ is altered, it’s nice to be able to go from a given power spectrum to some concrete position space maps whose statistics realize a given power spectrum.

In this post we answer the following question: if we’re handed a $$P(k)$$ which describes a Gaussian random field, how do we generate a simple fluctuation map which has the correlations described by this power spectrum?

We start by reviewing some of the theory behind random Gaussian fields, which I found difficult to find in any one place in the literature, and then move on to the practical construction of our maps using Mathematica. A sample notebook is included at the end of the post. The methods used for constructing the realization are by no means optimal; our goal is just to understand a quick and dirty method of generating realizations.

Theory Background

We start by reviewing the necessary theory background.

Correlators

First, we set notation and review basic facts about equal time correlators of scalar quantities. We focus only on two-point functions.

A homogeneous, isotropic position space two-point corrrelator only depends on the magnitude of the distance between the points of interest: $$\langle \phi(\bfx)\phi(\bfy)\rangle=A(|\bfx-\bfy|)$$. Due to this fact, the Fourier transformed correlator is proportional to a delta function and only depends on the magnitude of the momentum (i.e. wavevector):

$\langle \phi_{\bfk}\phi_{\bfp}\rangle\equiv\int\rd^{d}x\rd^{d}y\, e^{-i\bfk\cdot\bfx-i\bfp\cdot\bfy}\langle \phi(\bfx)\phi(\bfy)\rangle=P(k)\tilde{\delta}^{d}(\bfk+\bfp)$

We then put a prime on the momentum space correlator to indicate that we’ve dropped the delta function and associated $$2\pi$$ factors:

$\langle \phi_{\bfk}\phi_{-\bfk}\rangle'=P(k)$

We define $$\tilde{\delta}(\bfk+\bfp)\equiv (2\pi)^{d}\delta^{d}(\bfk+\bfp)$$ and later use $$\tilde{\bfk}\equiv \bfk/(2\pi)$$ to minimize the number of explicit $$2\pi$$’s appearing.

Two functional forms of the power spectrum are particularly important: white noise and scale invariant power spectra.

White Noise: $$P(k)\propto C$$ with $$C$$ constant. This corresponds to a position space correlator of the form $$\langle \phi(\bfx)\phi(\bfy)\rangle\propto C\delta^{d}(\bfx-\bfy)$$, i.e. there are no correlations between fluctuations at separated points if the power spectrum is described by white noise.

Scale Invariant: $$P(k)\propto k^{-d}$$ in $$d$$ spatial dimensions. A scale invariant power spectrum corresponds to fluctuations which have the same correlations no matter how separated the two points are:

${\rm Scale \ Invariant}\Longleftrightarrow \langle \phi(\lambda\bfx)\phi(\lambda\bfy)\rangle=\langle \phi(\bfx)\phi(\bfy)\rangle\ {\rm for \ all} \ \lambda$

Probability Distributions: Functions and Functionals

We now review how probability distributions work for continuous, Gaussian random fields. We start with a brief review of the analogous theory for continuous random variables.

Probability Distribution Functions (PDF): Given some continuous random variable $$q$$, the probability that $$q$$ lies in some range $$\rd q$$ is given by $$\P(q)\rd q$$, where $$\P(q)$$ is the PDF. The PDF is normalized as $$\int\P(q)\rd q=1$$ and the expectation value of $$q^{n}$$ is given by

$\langle q^n\rangle= \int\rd q\, \P(q)q^n$

The Gaussian PDF is of primary importance. For $$n$$ random variables, $$q_i$$, the Gaussian PDF takes on the form:

$\P^{\rm Gauss.}(q_1,\ldots,q_n)=\sqrt{\frac{\det (A)}{(2\pi)^n}}\exp\left[-\frac{1}{2}q_iA_{ij}q_j\right]\ ,$

where $$A_{ij}$$ is some positive definite matrix (and sums over repeated indices are implicit). The corresponding Gaussian two-point statistics are then given by:

$\langle q_i q_j\rangle =\left(\prod_{k=1}^{n}\int \rd q_k\right)\P^{\rm Gauss.}(\vec{q})q_iq_j=A^{-1}_{ij}$

where $$A^{-1}_{ij}$$ is the inverse of the matrix $$A_{ij}$$.

Probability Distribution Functionals are perhaps less familiar, but they are exactly the same idea as PDFs simply applied to the case of random fields such as $$\phi(\bfx)$$, rather than random variables such as $$q_i$$. We will only focus on the case of Gaussian functionals whose form and properties are defined in exact analogy to the PDF case. Starting with a position space field $$\phi(\bfx)$$, the Gaussian functional $$\P(\phi(\bfx))$$ is written as

${\cal P}[\phi(\bfx)]\propto \exp\left[-\frac{1}{2}\int\rd^dx\rd^dy\, \phi(\bfy)A(\bfy,\bfz)\phi(\bfz) \right]$

where we have omitted the normalization factor. Correlators are then computed by a path integral over all field values $$\phi(\bfx)$$. The outcome of the integral is determined in analogy to the PDF case. For instance, the two-point correlator is:

$\langle \phi(\bfx)\phi(\bfy)\rangle\equiv \int\D\phi\, \P[\phi(\bfx)]\phi(\bfy)\phi(\bfz)=A^{-1}(\bfy,\bfz)\ ,$

where $$A^{-1}(\bfy,\bfz)$$ obeys

$\int\rd^dy\, A(\bfx,\bfy)A^{-1}(\bfy,\bfz)=\delta^{d}(\bfx-\bfz)\ .$

White noise corresponds to the simple case where $$A(\bfx,\bfz)=C\delta^{d}(\bfx-\bfz)$$ and so $$A^{-1}(\bfx,\bfz)=\frac{1}{C}\delta^{d}(\bfx-\bfz)$$, for instance. Assuming spatial homogeneity and isotropy, it is straightforward to Fourier transform the above formulas in order to derive their momentum space analogues:

\begin{align} \P[\phi_{\bfk}]&\propto \exp\left[-\frac{1}{2}\int\rd^{d}\tilde{q}\rd^{d}\tilde{p}\, \phi_{\bfq}A(q)\tilde{\delta}^{d}(\bfq+\bfp)\phi_{\bfp}\right]\nn \langle\phi_{\bfk}\phi_{\bfk'}\rangle&=\frac{\tilde{\delta}^{d}(\bfk+\bfk')}{A(k)} \end{align}

Clearly, the function $$A(k)$$ above is simply the inverse of the power spectrum: $$P(k)=A^{-1}(k)$$. By adding terms cubic and higher order in $$\phi_\bfk$$ to the exponential in $$\P[\phi_\bfk]$$, the distribution is made non-Gaussian. The search for non-Gaussian signatures in Cosmological statistics is part of cutting-edge experimental efforts, but further discussion of these issues falls outside of the scope of this post.

Implementation

We now demonstrate how to generate a realization, given $$P(k)$$. We first present the steps to the process and explain the logic. Then we show how to practically implement the steps in Mathematica.

The Steps

A realization of a Gaussian random field with power spectrum can be created via the following steps:

1. Start with a white noise field with unit amplitude $$\varphi_{\bfk}$$ obeying $$\langle \varphi_{\bfk}\varphi_{-\bfk}\rangle'=1$$.
2. Generate a position space realization of the white noise, denoted by $$R_{\rm white}(\bfx)$$. That is, $$R_{\rm white}(\bfx)$$ is a particular map showing the values of $$\varphi$$ at various positions $$\bfx$$ and for which $$\langle\varphi(\bfx)\varphi(\bfy)\rangle'=\delta^{d}(\bfx-\bfy)$$.
3. Fourier transform the realization: $$R_{\rm white}(\bfx)\longrightarrow R_{\rm white}(\bfk)$$.
4. Multiply $$R_{\rm white}(\bfk)$$ by the square root of the power spectrum to create $$R_P(\bfk)=P^{1/2}(k)R_{\rm white}(\bfk)$$.
5. Fourier transform $$R_P(\bfk)$$ back to positio space to get the desired realization: $$R_P(\bfx)=\int\rd^d \tilde{k}\, e^{i\bfk\cdot\bfx} R_P(\bfx)$$.

The logic behind the above steps is pretty simple. If we’re given a white noise field $$\phi_\bfk$$ which has $$\langle \varphi_{\bfk}\varphi_{-\bfk}\rangle'=1$$ , then we can create a field with the desired power spectrum simply by defining $$\phi_\bfk\equiv P^{1/2}(k)\varphi_\bfk$$. The new field has the desired statistics, by construction.

Therefore, given a realization of the white noise field , we just have to multiply its Fourier space realization by and inverse Fourier transform to get a position space realization of $$\phi(\bfx)$$. The above is practical, because making a white noise realization is easy.

These steps are idealized, as we’ll have to make the realization on a discrete grid; we won’t be able to work with continuous variables. Practical complications which arise from working on a lattice are covered in the following sections.

From Continuous to Discrete

For concreteness, we focus on creating two-dimensional realizations on $$N\times N$$ grids, with $$N$$ an even number. The generalization to higher dimensions is straightforward. In passing to a lattice, our continuous fields $$\phi(\bfx)$$ and $$\phi_\bfk$$ are replaced by the discrete quantities $$\phi^{\bfx}_{ab}$$ and $$\phi^{\bfk}_{ab}$$, respectively, where $$a,b\in\{0,\ldots,N-1\}$$. $$\phi^{\bfx}_{ab}$$ is the value of the field at the point $$(a,b)$$ on the grid. The continuous Fourier transform which related $$\phi^{\bfx}_{ab}$$ and $$\phi^{\bfk}_{ab}$$ is replaced by a discrete one:

\begin{align} \phi^{\bfk}_{ab}&=\sum_{c,d=0}^{N-1}\exp\left(-ix_ck_a-ix_dk_b\right)\phi^{\bfx}_{cd}\nn \phi^{\bfx}_{ab}&=\frac{1}{N^2}\sum_{c,d=0}^{N-1}\exp\left(ix_ck_a+ix_dk_b\right)\phi^{\bfk}_{cd}\nn \end{align}

where $$x_a\equiv a$$ and $$k_a\equiv \frac{2\pi a}{N}$$. Given real $$\phi^{\bfx}_{ab}$$ and even $$N$$, $$\phi^{\bfk}_{ab}$$ has the following properties

• $$\phi^{\bfk}_{ab}=\phi^{\bfk}_{(a+\alpha N)b}=\phi^{\bfk}_{a(b+\alpha N)}$$ for any integer $$\alpha$$
• $$\phi^{*\bfk}_{ab}=\phi^{\bfk}_{-a-b}$$

The above further imply the following important property: $$\phi^{*\bfk}_{a(N/2+b)}=\phi^{\bfk}_{a(N/2-b)}$$ and similar with the other index. and similar with the other index.

Building a Realization

Now we walk through the steps involved in actually building a realization in Mathematica.

White Noise

First, we construct a grid of white noise. When passing to the grid, the unit white noise probability distribution functional becomes the following ordinary PDF:

$\P[\varphi(\bfx)]\longrightarrow \P[\varphi^{\bfx}_{ab}]=\prod_{c,d=0}^{N-1}\frac{\exp\left[-\frac{1}{2}\left(\phi^{\bfx}_{cd}\right)^2\right]}{\sqrt{2\pi}}$

Building a white noise realization is therefore simple: the value of at every grid point is just randomly chosen from a Gaussian distribution. An example of such a realization is shown in Fig. 1 below.

In Mathematica, the white noise array can be generated via:

 WhiteNoise=RandomVariate[NormalDistribution[],{Nsize,Nsize}] 

where Nsize is the size of the grid (Nsize=1024 in Fig. 1). WhiteNoise is our realization $$R_{\rm white}(\bfx)$$ and building $$R_{\rm white}(\bfk)$$ is very simple: we just pass WhiteNoise through the built-in Fourier function:

 WhiteNoiseFourier=Fourier[WhiteNoise] 

Constructing the Fourier Realization

The Fourier space realization then needs to be multiplied by the square root of the power spectrum. This is subtle, as the resulting field $$\phi^{\bfk}_{ab}$$ needs to obey the conditions we derived in the previous section, particularly $$\phi^{*\bfk}_{a(N/2+b)}=\phi^{\bfk}_{a(N/2-b)}$$. .

This is important because if we did the naive construction and took the white noise realization values $$\varphi^{\bfk}_{ab}$$ and simply multiplied by $$P^{1/2}(k)$$ with $$k=\frac{2\pi}{ N}\sqrt{a^2+b^2}$$ to build $$\phi^{\bfk}_{ab}=P^{1/2}(k)\varphi^{\bfk}_{ab}$$, then we'd find $$\phi^{*\bfk}_{a(N/2+b)}\neq\phi^{\bfk}_{a(N/2-b)}$$, due to the power spectrum factor. If we then transformed the naive $$\phi^{\bfk}_{ab}$$, we would be led to an imaginary $$\phi^{\bfx}_{ab}$$. We can get around this issue by building the array $$\Phi^{\bfk}_{ab}$$ whose components are given by:

$\Phi^{\bfk}_{ab}= \begin{cases} P^{1/2}(k)\varphi^{\bfk}_{ab} & a,b\le N/2\\ P^{1/2}(k)\varphi^{\bfk}_{a(b-N)} & a\le N/2, b>N/2\\ P^{1/2}(k)\varphi^{\bfk}_{(a-N)b} & a> N/2, b\le N/2\\ P^{1/2}(k)\varphi^{\bfk}_{(a-N)(b-N)} & a,b> N/2 \end{cases}$

In the second line above, the $$k$$ in $$P^{1/2}(k)\varphi^{\bfk}_{a(b-N)}$$ is evaluated at $$k=\frac{2\pi}{N}\sqrt{a^2+(b-N)^2}$$ and similar for other lines. The factors of $$N$$ can also be removed in the indices of the $$\varphi^{\bfk}$$'s, using the previously mentioned properties of $$\varphi^{\bfk}$$.

Note that these manipulations effectively result in the the $$a,b$$ indices running over the values $$\{0,\ldots, N/2,-N/2+1,\ldots\}$$, rather than $$\{0,\ldots, N-1\}$$. In the limit $$P(k)\to 1$$, Fourier transforming $$\Phi^{\bfk}_{ab}$$ yields the original position space white noise realization we started with, $$\varphi^{\bfk}_{ab}$$, , despite the above manipulations. The position space field is then built via

$\phi^{\bfx}_{ab}=\frac{1}{N^2}\sum_{c,d=0}^{N-1}\exp\left(ix_ck_a+ix_dk_b\right)\Phi^{\bfk}_{cd}\ .$

This field realizes the desired spectrum. The construction of can be carried out in Mathematica as follows:

1. We first build the vector $$k_a=\frac{2\pi}{N}(0,\ldots, N/2,-N/2+1,\ldots, -1)$$ as:
 kVector=N[(2π)/Nsize*Join[Range[0,Nsize/2],Range[-Nsize/2+1,-1]]] 
2. Next, build the power spectrum function. Take it to be of the power law form: $$P(k)\propto k^{-n}$$. Because this $$P(k)$$ is divergent at $$k=0$$ (assuming $$n>0$$), we will set the value of the power spectrum to zero at this one point. This is implemented in Mathematica as:
 PowerLawPowerSpectrum[n_,k_]:= If[k==0,0,1/Abs[k^n]] 
3. We want to evaluate $$P^{1/2}(k)$$ at all points $$(k_a,k_b)$$ with $$k=\sqrt{k_a^2+k_b^2}$$. This array is constructed using the Outer command:
 sqrtPowerSpectrumArray=Outer[N[Sqrt[PowerLawPowerSpectrum[n,Norm[{##}]]]]&,kVector,kVector] 
where $$n$$ should be replaced by the desired number determining the power law.
4. Finally, we mutiply the $$a,b$$ entry of sqrtPowerSpectrumArray by the $$a,b$$ entry of WhiteNoiseFourier to get $$\Phi^{\bfk}_{ab}$$:
PowerSpectrumRealizationFourier=sqrtPowerSpectrumArray*WhiteNoiseFourier 
Note that the above multiplication is carried out with the Times function, not Dot.

Plotting the Realization

Finally, we transform from $$\Phi^{\bfk}_{ab}\longrightarrow \phi^{\bfx}_{ab}$$ and plot the result. This is done in a single command as:

MatrixPlot[Re[InverseFourier[PowerSpectrumRealizationFourier]]]  

Above, the real part of the transform was taken using Re because numerical errors induce tiny imaginary components in the final calculation of $$\phi^{\bfx}_{ab}$$. The results of the realization for various power spectra can be seen in Fig. 2.

An Example in $$d=3$$

Finally, we present the results of a higher dimensional realization: Fig. 3 is a GIF which shows successive slices of a $$d=3$$ realization for a power spectrum $$P(k)=k^{-4}$$.