Panphasia: Reproducible Gaussian Phase Field

• Panphasia is a discrete Gaussian white-noise field defined to supply public and reproducible phase information for cosmological simulations.
• It features a hierarchical octree structure with 50 levels, allowing efficient on-demand evaluation of any simulation subregion without explicit storage.
• Its implementation using the MRGK5-93 pseudorandom generator and an S8 Legendre block basis ensures accurate phase reconstruction for multi-scale resimulations.

Searching arXiv for the cited Panphasia papers and closely related records. Panphasia is a very large predefined discrete realization of a Gaussian white-noise field, designed specifically to provide public, reproducible phase information for Gaussian initial conditions in cosmological simulations and resimulations of structure formation. Its defining features are a hierarchical octree structure with 50 octree levels fully populated, deterministic evaluation of arbitrary subregions on demand rather than explicit storage, and a compact descriptor format that identifies the phase volume used by a simulation. In this framework, Panphasia supplies the random phases, while the cosmological parameters, linear power spectrum, and any resimulation-region specification determine the amplitudes and physical interpretation of the initial conditions (Jenkins et al., 2013, Jenkins, 2013).

1. Scientific role and phase specification

Panphasia was introduced to make phase choice objective, public, and reproducible. In conventional practice, different simulation groups often generated distinct private realizations of Gaussian white noise, so reproducing a published result required more than matching the power spectrum: it required exactly the same realization of phases. Panphasia replaces that ad hoc situation with a single global white-noise realization from which any simulation or resimulation draws a specified subvolume (Jenkins, 2013).

The underlying rationale is that a white-noise field contains the phase information for any Gaussian field. Once a white-noise realization has been fixed, one may convolve it with a real kernel whose Fourier amplitude is $\sqrt{P(\mathbf{k})}$ to obtain a Gaussian field with the desired power spectrum while preserving phases. The construction therefore decouples phase specification from cosmology. A common misconception is that Panphasia itself defines a cosmology; it does not. It specifies phases only, and full reproducibility still requires the cosmological parameters, the linear power spectrum, and, for resimulations, the locations of the objects or regions of interest (Jenkins, 2013, Jenkins et al., 2013).

This objective phase specification is particularly important for multi-scale resimulations. Because fluctuations at all relevant physical scales are already present in one predefined hierarchical field, refinement does not invent additional short-wavelength structure ad hoc. Instead, deeper octree levels reveal smaller-scale structure already implicit in the same realization. This gives repeated zooms at increasing resolution a consistency with the parent simulation and with one another that is, in the paper’s formulation, “by construction” (Jenkins, 2013).

2. Hierarchical and mathematical construction

Panphasia is organized as an octree. At octree level $l$, each cell is labeled by integer Cartesian coordinates

$$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$

and, in the construction paper’s notation, cells at level $l$ have side length

$\Delta_l = \frac{L}{2^l}.$

A simulation phase volume is defined by selecting a cube of octree cells at some chosen level. In the released software, this selected cube is treated as a cubic 3-torus with periodic boundary conditions; although the general descriptor convention allows cuboids, the public code supports only the cubic case (Jenkins et al., 2013, Jenkins, 2013).

The mathematical basis is the invariance of Gaussian white noise under orthonormal change of basis. If $\{f_i(\mathbf{x})\}$ is a complete orthonormal set, then a white-noise field can be expanded as

$\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$

with coefficients $C_j$ that are again independent Gaussian variables. Panphasia exploits this by using an octree-localized orthogonal basis built from Legendre blocks and octree basis functions (Jenkins, 2013).

The 3D Legendre block functions are localized products of Legendre polynomials,

$p_{j_1j_2j_3}({\bf x}) = \begin{cases} \displaystyle \prod_{i=1}^3 (2j_i+1)^{1/2}P_{j_i}(2x_i) & \text{if %%%%5%%%%;} \ 0 & \text{otherwise,} \end{cases}$

with orthonormality over the cell volume. The construction paper considers three self-representing block sets: $S_1$, $l$0, and $l$1. “Self-representing” means that a block at one level can be represented exactly by blocks of the same set in the eight children one level deeper, which is essential for a hierarchical basis. The practical conclusion of the paper is that $l$2 is the best compromise: $l$3 reproduces large-scale power too poorly, with low-$l$4 errors that scale as $l$5, whereas $l$6 gives $l$7 convergence at low $l$8 and $l$9 would be much more expensive computationally (Jenkins, 2013).

In the chosen $$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$0 realization, each octree cell carries 64 localized functions. Eight are the Legendre blocks returned directly to the user-facing interface,

$$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$1

and the remaining 56 are the true octree basis functions used internally in the hierarchical construction. A useful property for resimulation work is that these $$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$2 octree basis functions have vanishing zeroth and first moments, which reduces boundary artifacts when they are assigned to nested Fourier grids (Jenkins, 2013, Jenkins et al., 2013).

3. Pseudorandom sequence, Gaussian mapping, and formal field definition

Panphasia assigns Gaussian coefficients to its basis functions using the MRGK5-93 pseudorandom number generator of L’Ecuyer et al. (1993). The generator’s state is a 5-component integer vector modulo

$$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$3

and the crucial implementation property is efficient jumping: one can advance by $$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$4 steps through matrix powering at cost $$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$5. This makes effectively random access into an astronomically long sequence practical and is central to mapping a deterministic pseudorandom stream reproducibly into octree cells (Jenkins, 2013).

The generator period is

$$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$6

with a subperiod

$$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$7

The paper discusses possible long-range correlations associated with the subperiod and concludes that they are irrelevant for realistic uses of Panphasia. Uniform deviates are transformed into Gaussian variables using the Box–Muller transform. Because the generator is discrete, the far Gaussian tail would otherwise be artificially truncated when the uniform variate is extremely small; the implementation therefore replaces $$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$8 by $$(i_x,i_y,i_z), \qquad 0 \le i_x,i_y,i_z \le 2^l - 1,$$9, with $l$0 taken from a branch reached by a very large jump of $l$1 (Jenkins, 2013).

The sequence is mapped to the octree level-by-level, cell-by-cell, and basis-function-by-basis-function. The first eight Gaussian numbers are reserved for the root-cell Legendre blocks. Thereafter, each octree cell receives 64 Gaussian numbers: 56 for the octree basis functions and 8 extra numbers for an auxiliary independent field used in initial-condition construction. The index mapping is

$l$2

This yields the formal field definition given in the construction paper,

$l$3

With 50 octree levels, the field spans a factor of more than $l$4 in linear scale. In the public-code convention, the root cell has physical size $l$5, and the smallest level-50 cells correspond to masses around $l$6 (Jenkins, 2013).

4. Descriptors, coordinates, and user-facing software

Panphasia’s operational innovation is the phase descriptor: a text string that compactly specifies the simulation volume and therefore its phases. The canonical cubic form is

$l$7

where the fields are the descriptor/version tag, octree level, absolute coordinates of the corner cell nearest the origin, cube side length in level-$l$8 cells, a check number, and a human-readable name. A typical example is

$l$9

The construction paper also defines a cuboid syntax using D(...), but the released software assumes that the descriptor defines a cube rather than a general cuboid (Jenkins, 2013, Jenkins et al., 2013).

The check integer is intended to catch transcription errors, because even a tiny change in a descriptor completely alters the phases. In the public code, utilities include parse_descriptor, compose_descriptor, validate_descriptor, and generate_random_descriptor. A notable behavior of validate_descriptor is that if the check-digit field is -999 and the rest is correct, it returns the proper check digit. The final name string is human-readable, may contain no spaces, and can be up to 20 characters long (Jenkins et al., 2013).

The released package is PANPHASIA_code_V1.000.tar.gz, unpacking to PANPHASIA_public_code with subdirectories src, demo, and cps. The src directory contains panphasia_routines.f and generic_lecuyer.f90; demo contains main.f and an example makefile; and cps contains a modified version of the public serial initial-conditions code of Crocce et al. (2006), referred to as the CPS code (Jenkins et al., 2013).

For standard periodic cosmological simulations, the user guide presents Panphasia as a black-box function accessed chiefly through two subroutines: $C_j$9 The constraint

$\Delta_l = \frac{L}{2^l}.$0

must hold, where $\Delta_l = \frac{L}{2^l}.$1 is the descriptor side length, so that there is a one-to-one correspondence between octree cells and mesh cells. For each requested cell, panphasia_cell_properties returns cell_prop(9): the eight Legendre block coefficients and a ninth independent Gaussian pseudorandom number. All nine returned quantities are independent Gaussian pseudorandom numbers with

$\Delta_l = \frac{L}{2^l}.$2

For resimulations, the more general interface $p_{j_1j_2j_3}({\bf x}) = \begin{cases} \displaystyle \prod_{i=1}^3 (2j_i+1)^{1/2}P_{j_i}(2x_i) & \text{if %%%%5%%%%;} \ 0 & \text{otherwise,} \end{cases}$0 allows explicit choice of octree levels and of the relative origin within the periodic 3-torus (Jenkins et al., 2013).

5. Use in cosmological boxes and nested resimulations

For a uniform cosmological simulation, one selects a cubic region made of $\Delta_l = \frac{L}{2^l}.$3 octree cells at some level $\Delta_l = \frac{L}{2^l}.$4, then chooses a Fourier mesh commensurate with a finer octree level so that grid points correspond one-to-one with octree cells. The Panphasia interface returns, for each finest-level cell, the 8 $\Delta_l = \frac{L}{2^l}.$5 Legendre-block coefficients plus the ninth independent Gaussian number. The initial-condition generator then places the eight coefficient fields on eight real-space grids, Fourier-transforms them, multiplies by the Fourier transform of the corresponding Legendre block, applies the half-grid phase correction required by the grid-centre versus grid-corner convention, and sums them to obtain the Panphasia-derived white-noise contribution in $\Delta_l = \frac{L}{2^l}.$6-space (Jenkins, 2013).

Because the 8 $\Delta_l = \frac{L}{2^l}.$7 Legendre blocks at a single depth are not a complete basis, finite-depth reconstruction is not exactly white noise. The missing power is restored by adding the independent ninth field with spectrum

$\Delta_l = \frac{L}{2^l}.$8

This yields a true white-noise spectrum with isotropic high-$\Delta_l = \frac{L}{2^l}.$9 cutoff, although it also means that the smallest-scale phases are not perfectly those of Panphasia. The construction paper characterizes the associated phase error through a smoothed fractional RMS measure $\{f_i(\mathbf{x})\}$0 and concludes that $\{f_i(\mathbf{x})\}$1 reconstructs the true phases far more accurately than $\{f_i(\mathbf{x})\}$2 (Jenkins, 2013).

For zoom or resimulation initial conditions, the central principle is that each octree basis function must be assigned to exactly one nested Fourier grid. The paper gives four heuristic placement rules: any octree basis function that can be included should be included; except for the outermost grid, only whole octree basis functions can be assigned to a grid; subject to a boundary rule, each basis function should be placed on the innermost allowed grid; and, on inner grids, except for the smallest basis functions, a basis function should not be placed within a perpendicular distance less than $\{f_i(\mathbf{x})\}$3 times its own size from a grid boundary. These rules use the vanishing zeroth and first moments of the basis functions to suppress edge artifacts and improve sampling on the innermost feasible grid (Jenkins, 2013).

The user guide also documents a concrete integration path through the modified CPS code. The principal structural change was to separate phase generation from the computation of Zel’dovich displacements in $\{f_i(\mathbf{x})\}$4-space, allowing a conditional choice between the original phase generator and Panphasia-derived phases. The new routine initialise_wnf, in set_up_wn_field.f, computes the Panphasia-based white-noise field. Because the CPS code stores six values per grid cell for 2LPT, those arrays are reused as temporary storage, so the Panphasia modification does not increase memory requirements. Since Panphasia returns nine values per cell, however, the property routine must be called twice per grid cell, and the nine values are combined into a single $\{f_i(\mathbf{x})\}$5-space field in two stages (Jenkins et al., 2013).

6. Validation, performance characteristics, and limitations

The user guide reports basic statistical tests on a $\{f_i(\mathbf{x})\}$6 grid for a randomly chosen descriptor. For each of the 9 returned fields, the expected properties are a Gaussian one-point distribution with $\{f_i(\mathbf{x})\}$7 and $\{f_i(\mathbf{x})\}$8, vanishing two-point autocorrelation at non-zero lag, and vanishing cross-correlations between distinct outputs. For the $\{f_i(\mathbf{x})\}$9 field, using 61 bins covering $\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$0 to $\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$1, the reported result is

$\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$2

against an expected value near 61, which the guide interprets as satisfactorily Gaussian. For the full set of cross-correlations and autocorrelations, the values followed the expected Gaussian with variance $\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$3, and the absence of outliers beyond $\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$4 implied an upper bound on the maximum correlation of roughly

$\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$5

The guide is explicit that these tests are necessary but not sufficient and encourages additional randomness tests; their primary purpose was to catch software errors, and deliberately introduced mistakes were readily detected (Jenkins et al., 2013).

The construction paper extends validation to full N-body comparisons. Its strongest test compares a very expensive reference cosmological run based on a $\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$6 Fourier transform with cheaper cosmological runs using $\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$7 or $\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$8 transforms and with zoom simulations built from nested grids but using Panphasia. The main conclusions are that $\delta(\mathbf{x}) = \sum_j C_j f_j(\mathbf{x}),$9 is dramatically better than $C_j$0; that cosmological initial conditions made from Panphasia with $C_j$1 and sufficiently large Fourier grids reproduce halo bulk properties essentially at sub-percent level; and that a zoom simulation with 4 nested Fourier meshes, each $C_j$2, using $C_j$3 and $C_j$4, reproduces the reference halo extremely well, with $C_j$5, $C_j$6, and $C_j$7 differing by less than a percent, often much less (Jenkins, 2013).

Several practical limitations are emphasized. The released evaluator is serial and uses only a small amount of memory; for parallel use, the guide recommends multiple private instances rather than a native distributed-memory implementation. Access pattern matters strongly because the code caches up to eight octree cells at each tree level near the current spatial location: random access is expensive, while raster or locality-preserving access is much faster. The construction paper reports a typical evaluation cost of about $C_j$8 per cell, with white-noise generation around 20% of total initial-condition generation time. Floating-point outputs can vary slightly with compiler, optimization flags, and order of operations, though the guide states that such differences are very small and not significant for cosmological initial-condition generation (Jenkins et al., 2013, Jenkins, 2013).

A final limitation is conceptual rather than algorithmic. Publishing a descriptor makes the phases citable and reconstructible, but it does not by itself specify a complete simulation. Reproduction additionally requires the cosmological parameters, the linear power spectrum, and, for resimulations, the region or object definition. Within those bounds, Panphasia’s main contribution is to turn phase information from a hidden implementation detail into a compact, public, deterministic component of cosmological numerical methodology (Jenkins et al., 2013, Jenkins, 2013).