Power Spectra

Though SWIFT includes functionality for calculating power spectra, this tool runs on-the-fly, and as such after a run has completed you may wish to create a number of more non-standard power spectra.

These tools are available as part of the swiftsimio.visualisation.power_spectrum module. Making a power spectrum consists of two major steps: depositing the particles on grid(s), and then binning their fourier transform to get the one-dimensional power.

Depositing on a grid

Depositing your particles on a grid is performed using swiftsimio.visualisation.power_spectrum.render_to_deposit(). This function performs a nearest-grid-point (NGP) of all particles in the provided particle dataset. For example:

from swiftsimio import load
from swiftsimio.visualisation.power_spectrum import render_to_deposit

data = load("cosmo_volume_example.hdf5")

gas_mass_deposit = render_to_deposit(
    data.gas,
    resolution=512,
    project="masses",
    parallel=True,
)

The specific field being depositied can be controlled with the project keyword argument. The resolution` argument gives the one-dimensional resolution of the 3D grid, so in this case you would recieve a 512x512x512 grid. Note that the gas_mass_deposit is a cosmo_array, and as such includes cosmological and unit information that is used later in the process.

Generating a power spectrum

Once you have your grid deposited, you can easily generate a power spectrum using the deposition_to_power_spectrum() function. For example, using the above deposit:

from swiftsimio.visualisation.power_spectrum import deposition_to_power_spectrum

wavenumbers, power_spectrum, _ = deposition_to_power_spectrum(
    deposition=gas_mass_deposit,
    boxsize=data.metadata.boxsize,
)

This power spectrum can then be plotted. Units are included on both the wavenumbers and the power spectrum. Cross-spectra are also supported through the cross_deposition keyword, but by default this generates the auto power. Wavenumbers are calculated to be at the weighted mean of the k-values in each bin, rather than representing the center of the bin.

More complex scenarios

In a realistic simualted power spectrum, you will need to perform ‘folding’ to achieve a viable dynamic range within an achievable memory footprint. Consider, for instance, a 1 Gpc simulation volume with a 1 kpc resolution limit. In that case, you would need a deposition grid with \(10^{18}\) cells, amounting to a 4 EB (yes, exabyte!) memory footprint. That is, of course, not realistic!

As in a power spectrum we are only interested in the periodicity of the system, we can fold it back in on itself during the rendering process. The position of the particle in the box is set to be:

\(x_i' := \left( x_i \% \frac{L}{2^{n}} \right) \frac{2^{n}}{L}\)

where \(L\) is the box-size and \(n\) is some integer greater than or equal to zero. This allows you to probe modes in the reduced box-length \(L / 2^{n}\) with the same fixed resolution deposition buffer.

The folding parameter is available for both render_to_deposit and deposition_to_power_spectrum, but it may be easier to use the utility functions provided for automatically stitching together the folded spectra. The function folded_depositions_to_power_spectrum() allows you to do this easily:

import unyt as u
from swiftsimio.visualisation.power_spectrum import folded_depositions_to_power_spectrum
from swiftsimio.objects import cosmo_array

folded_depositions = {}

for folding in [x * 2 for x in range(5)]:
    folded_depositions[folding] = render_to_deposit(
        data.gas,
        resolution=512,
        project="masses",
        parallel=True,
        folding=folding,
    )

bins, centers, power_spectrum, foldings = folded_depositions_to_power_spectrum(
    depositions=folded_depositions,
    boxsize=data.metadata.boxsize,
    number_of_wavenumber_bins=128,
    wavenumber_range=cosmo_array(
        [1e-2, 1e2],
        u.Mpc**-1,
        comoving=True,
        scale_factor=data.metadata.a,
        scale_exponent=-1,
    ),
    log_wavenumber_bins=True,
    workers=4,
    minimal_sample_modes=8192,
    cutoff_above_wavenumber_fraction=0.75,
    shot_noise_norm=len(gas_mass_deposit),
)

The ‘used’ foldings of the power spectrum are shown in the foldings return vaule, which is an array containing the folding that was used for each given bin. This is useful for debugging and visualisation.

There are a few crucial parameters to this function:

1. workers is the number of threads to use for the calculation of the Fourier transforms.

2. minimal_sample_modes is the minimum number of modes that must be present in a bin for it to be included in the final power spectrum. Generally for a big simulation you want to set this to around 10,000, and this number is ignored for the lowest wavenumber bin.

3. cutoff_above_wavenumber_fraction is the fraction of the individual fold’s (as represented by the FFT itself) maximally sampled wavenumber. Ignored for the last fold, and we always cap the maximal wavenumber to the Nyquist frequency.

4. shot_noise_norm is the number of particles in the simulation that contribute to the power spectrum. This is used to normalise the power spectrum to the shot noise level. This is very important in this case because of the use of NGP deposition.

Foldings are stitched using a simple method where the ‘better sampled’ foldings are used preferentially, up to the cutoff value.