How Does Powerbox Work?¶

It may be useful to understand the workings of powerbox to some extent – either to diagnose performance issues or to understand its behaviour in certain contexts.

The basic algorithm (for a Gaussian field) is the following:

Given a box length \(L\) (parameter boxlength) and number of cells along a side, \(N\) (parameter N), as well as Fourier convention parameters \((a,b)\), determine wavenumbers along a side of the box: \(k = 2\pi j/(bL)\), for \(j\in (-N/2,..., N/2)\).
From these wavenumbers along each side, determine the magnitude of the wavenumbers at every point of the \(d\)-dimensional box, \(k_j= \sqrt{\sum_{i=1}^d k_{i,j}^2}\).
Create an array, \(G_j\), which assigns a complex number to each grid point. The complex number will have magnitude drawn from a standard normal, and phase distributed evenly on \((0,2\pi)\).
Determine \(\delta_{k,j} = G_j \sqrt{P(k_j)}\).
Determine \(\delta_x = V \mathcal{F}^{-1}(\delta_k)\), with \(V = \prod_{i=1}^{d} L_i\).

For a Log-Normal field, the steps are slightly more complex, and involve determining the power spectrum that would be required on a Gaussian field to yield the same power spectrum for a log-normal field. The details of this approach can be found in Coles and Jones (1991) or Beutler et al. (2011).

One characteristic of this algorithm is that it contains no information below the resolution scale \(L/N\). Thus, a good rule-of-thumb is to choose \(N\) large enough to ensure that the smallest scale of interest is covered by a factor of 1.5, i.e., if the smallest length-scale of interest is \(s\), then use \(N = 1.5 L/s\).

The range of \(k\) used with this choice of \(N\) also depends on the Fourier Convention used. For the default convention of \(b=1\), the smallest scales are equivalent to \(k = \pi N/L\).